关于空间直角坐标系中三维速度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