直角坐标转换为地理坐标.DOC

直角坐标转换为地理坐标 设计算法将地理坐标转换为直角坐标，通过迭代次数和欧氏误差评价算法性能。总体思路是：先将地理坐标（理论值）转换为直角坐标（理论值），由于该算法没有误差，两坐标是等价的。再用设计的算法将直角坐标系转换为地理坐标（计算值），在同样的纬度和高度情况下，取出经度误差（理论值与计算值之差的绝对值）最大的点，将该点的地理坐标转换为直角坐标，该坐标与直角坐标理论值即为评价算法的误差。 当高度在-6300000m时，迭代次数和误差曲线如下图所示： 迭代次数： 误差： 18000个点中，发散点个数为：convergenum =0 当高度在-50000m到1000000m时，迭代次数和误差曲线如下图所示： 迭代次数： 误差： 18000个点中，发散点个数为：convergenum =0 当高度在-6378000m到-6300000m时，迭代次数和误差曲线如下图所示： 迭代次数： 误差： 18000个点中，发散点个数为：convergenum =0 主程序： R=6378.173*1000; r=6356.7523142*1000; f=(R-r)/R; w=linspace(-90,90,30); %纬度 i l=linspace(-180,180,30); %经度 j %h=linspace(-50000,1000000,20); %高度 k %h=linspace(-630000020); %高度 k h=linspace(-6378000,-6300000,20); %高度 k hhj=pi/180; %角度换算成弧度需乘的量 w=w*hhj; %纬度换算成弧度 l=l*hhj; %经度换算成弧度 ww=[]; %存放初始纬度数组 ll=[]; %存放初始经度数组 hh=[]; %存放初始高度数组 for i=1:30 %纬度 i for j=1:30 %经度 j for k=1:20 %高度 k ww(i,j,k)=w(i); ll(i,j,k)=l(j); hh(i,j,k)=h(k); end end end [xxe,yye,zze]=wlh_change_xyz(ww,ll,hh); %经纬高wlh转换为直角坐标系xyz [NEWw,NEWl,NEWh,flag_converge,Nkk]=xyz_change_wlh(xxe,yye,zze);%直角坐标系xyz转换为经纬高wlh [xxe1,yye1,zze1]=wlh_change_xyz(NEWw,NEWl,NEWh);%经纬高wlh转换为直角坐标系xyz convergenum=length(xxe(:,1,1))*length(xxe(1,:,1))*length(xxe(1,1,:))-sum(sum(sum(flag_converge))) wuchax=xxe-xxe1;%x的误差 wuchay=yye-yye1;%y的误差 wuchaz=zze-zze1;%z的误差 wucha=sqrt(wuchax.^2+wuchay.^2+wuchaz.^2);%欧氏的误差 wuchamax=[];%存放相同经度时的最大误差 Nkkmax=[];%存放相同经度时的最大迭代次数 for i=1:30 %纬度 i for k=1:20 %高度 k wuchamax(i,k)=max(abs(wucha(i,:,k)));%求放相同经度时的最大误差 Nkkmax(i,k)=max(Nkk(i,:,k));%求放相同经度时的最大迭代次数 end end %画图 plotw=w; ploth=h; pplotw=ones(size(plotw))*ploth; pploth=plotw*ones(size(ploth)); surf(pplotw,pploth,Nkkmax)%画迭代次数曲线 xlabel(Height(m)),ylabel(Geodetic latitude(degree)),zlabel(Iteration numbers) grid on ; surf(pplotw,pploth,wuchamax)%画误差曲线 xlabel(Height(m)),ylabel(Geodetic latitude(degree)),zlabel(Euclid distance errors(m)) --------------------------------------