Matlab中空间直角坐标系中三维速度Vxyz转Ven平面速度

关于空间直角坐标系中三维速度Vxyz转Ven平面速度，需要完成连个步骤的工作：①空间直角坐标XYZ转为大地坐标BLH；②
三维速度Vxyz转Ven平面速度；
一、空间直角坐标XYZ转为大地坐标BLH

二、三维速度Vxyz转Ven平面速度

三、正式Matlab代码
（1）工程文件下设置的文件目录，data下存放输入文件，result下存生成结果文件：

（2）data文件夹下存放的输入文件：

（3）输入文件的格式，第一列为站点名，第二列至第四列为对应XYZ坐标，第五列至第六列为对应Vx,Vy,Vz速度值：

（4）编写Matlab代码：
1、main_Vxyz2Ven.m

clc,clear

%///读入文件数据
file_address = ['data/','raw_end.txt'];                                    %加载文件路径
%[data_name,X,Y,Z,Vx,Vy,Vz]=textread(file_address,'%s%f%f%f%f%f%f');       %第二种打开文件方式

fid = fopen(file_address, 'r');                                            %以只读的方式打开文件
if (fid  == -1 )
    disp('Error opening file');
else
    c = textscan(fid, '%s %f %f %f %f %f %f');                             % 读取文件
    fclose(fid);
    data_name=c{1};      %第一列为站点
    X=c{2};              %第二列为坐标X
    Y=c{3};              %第三列为坐标Y
    Z=c{4};              %第二列为坐标Z
    Vx=c{5};             %第三列为速度Vx
    Vy=c{6};             %第二列为速度Vy
    Vz=c{7};             %第三列为速度Vz
    
    
    N1=0;
    N1=length(X);
    %///调用函数将原XYZ坐标转为BL大地坐标
    [B,L] = XYZtoBLH(X,Y,Z);
    
    
    %///XYZ方向速度转为NE平面速度速度
    [Venu] = VXYZ2NEU(Vx,Vy,Vz,B,L,N1);
    
    
    %///转为度
    B_END=B*180/pi;
    L_END=L*180/pi;
    
    
    %///结果写入文件
    cd Result_data
    fid2=fopen('endresult.txt','w');
    for i=1:N1
        data_NAME = data_name(i);
        data_NAME = cell2mat(data_NAME);
        fprintf(fid2,' %5s',data_NAME(1:4));
        fprintf(fid2,' %8.3f',L_END(i));             %经度
        fprintf(fid2,' %8.3f',B_END(i));             %纬度
        fprintf(fid2,' %8.6f',Venu(i,2));            %Vn
        fprintf(fid2,' %8.6f',Venu(i,1));            %Ve
        if i ~= N1
            fprintf(fid2,'\n');
        end
    end
    fclose(fid2);
    cd ..
    
    disp('*******转换程序运行结束*******');
end

2、XYZ2BLH.m

function [B,L] = XYZ2BLH(X,Y,Z)
%///WGS84坐标转换到大地经纬度
v=6378137;
e2=(1.0 / 298.257223563) * (2 - (1.0 / 298.257223563));

L = atan2(Y, X);
B0 = atan(Z ./sqrt(X .* X + Y .* Y));
while 1
    N = v ./ sqrt(1 - e2 .* sin(B0) .* sin(B0));
    H = Z ./ sin(B0) - N .* (1 - e2);
    B = atan(Z .* (N+H) ./ (sqrt(X .* X + Y .* Y) .* (N .* (1-e2)+H)));
    if abs(B - B0) < 1e-6
        break;
    end
	B0 = B;
	N = v ./ sqrt(1 - e2 .* sin(B) .* sin(B));
end

end

3、VXYZ2VNEU.m

function [Venu] = VXYZ2NEU(Vx,Vy,Vz,B_end,L_end,N1)
%B_end纬度，L_end经度
for num=1:N1
    B = B_end(num);
    L = L_end(num);
    V = [Vx(num);Vy(num);Vz(num)];
    T = [-sin(L) cos(L) 0;-sin(B)*cos(L) -sin(B)*sin(L) cos(B);cos(B)*cos(L) cos(B)*sin(L) sin(B)];
    Venu(num,:)=T*V;
end
end