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

速度Vxyz转Ven平面速度"/>

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');
elsec = textscan(fid, '%s %f %f %f %f %f %f');                             % 读取文件fclose(fid);data_name=c{1};      %第一列为站点X=c{2};              %第二列为坐标XY=c{3};              %第三列为坐标YZ=c{4};              %第二列为坐标ZVx=c{5};             %第三列为速度VxVy=c{6};             %第二列为速度VyVz=c{7};             %第三列为速度VzN1=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_datafid2=fopen('endresult.txt','w');for i=1:N1data_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));            %Vnfprintf(fid2,' %8.6f',Venu(i,1));            %Veif i ~= N1fprintf(fid2,'\n');endendfclose(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 1N = 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-6break;endB0 = B;N = v ./ sqrt(1 - e2 .* sin(B) .* sin(B));
endend

3、VXYZ2VNEU.m

function [Venu] = VXYZ2NEU(Vx,Vy,Vz,B_end,L_end,N1)
%B_end纬度，L_end经度
for num=1:N1B = 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