1、直角坐标转换为地理坐标设计算法将地理坐标转换为直角坐标，通过迭代次数和欧氏误差评价算法性能。总体思路是：先将地理坐标（理论值）转换为直角坐标（理论值） ，由于该算法没有误差，两坐标是等价的。再用设计的算法将直角坐标系转换为地理坐标（计算值） ，在同样的纬度和高度情况下，取出经度误差（理论值与计算值之差的绝对值）最大的点，将该点的地理坐标转换为直角坐标，该坐标与直角坐标理论值即为评价算法的误差。当高度在-6300000m 到 30000000m 时，迭代次数和误差曲线如下图所示：迭代次数：误差：18000 个点中，发散点个数为：convergenum =0当高度在-50000m到1000000m
2、时，迭代次数和误差曲线如下图所示：迭代次数：误差：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); %纬度 il=linspace(-180,180,30); %经度 j%h=linspace(-50000,1000000,20); %高度 k%h=linspace(-630000
3、0,30000000,20); %高度 kh=linspace(-6378000,-6300000,20); %高度 khhj=pi/180; %角度换算成弧度需乘的量w=w*hhj; %纬度换算成弧度l=l*hhj; %经度换算成弧度ww=; %存放初始纬度数组ll=; %存放初始经度数组hh=; %存放初始高度数组for i=1:30 %纬度 ifor j=1:30 %经度 jfor k=1:20 %高度 k ww(i,j,k)=w(i);ll(i,j,k)=l(j);hh(i,j,k)=h(k);endendendxxe,yye,zze=wlh_change_xyz(ww,ll,hh);
4、 %经纬高 wlh 转换为直角坐标系 xyzNEWw,NEWl,NEWh,flag_converge,Nkk=xyz_change_wlh(xxe,yye,zze);%直角坐标系 xyz 转换为经纬高 wlhxxe1,yye1,zze1=wlh_change_xyz(NEWw,NEWl,NEWh);%经纬高 wlh 转换为直角坐标系 xyzconvergenum=length(xxe(:,1,1)*length(xxe(1,:,1)*length(xxe(1,1,:)-sum(sum(sum(flag_converge)wuchax=xxe-xxe1;%x 的误差wuchay=yye-yye1
5、;%y 的误差wuchaz=zze-zze1;%z 的误差wucha=sqrt(wuchax.2+wuchay.2+wuchaz.2);%欧氏的误差wuchamax=;%存放相同经度时的最大误差Nkkmax=;%存放相同经度时的最大迭代次数for i=1:30 %纬度 ifor k=1:20 %高度 k wuchamax(i,k)=max(abs(wucha(i,:,k);%求放相同经度时的最大误差Nkkmax(i,k)=max(Nkk(i,:,k);%求放相同经度时的最大迭代次数endend%画图plotw=w;ploth=h;pplotw=ones(size(plotw)*ploth;pp
6、loth=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)-%经纬高 wlh 转换为直角坐标系 xyzfunction xxe,yye
7、,zze=wlh_change_xyz(w,l,h)R=6378.173*1000;r=6356.7523142*1000;f=(R-r)/R;hhj=pi/180; %角度换算成弧度需乘的量w=w*hhj; %纬度换算成弧度l=l*hhj; %经度换算成弧度o=acot(cot(w)/(1-f); %通过纬度求角度 oxxe=; %存放直角坐标系 x 的数组yye=; %存放直角坐标系 y 的数组zze=; %存放直角坐标系 z 的数组for i=1:length(w(:,1,1) %纬度 ifor j=1:length(w(1,:,1) %经度 jfor k=1:length(w(1,1,
8、:) %高度 km=R*cos(o(i,j,k)+h(i,j,k)*cos(w(i,j,k);xe=m*cos(l(i,j,k);ye=m*sin(l(i,j,k);ze=r*sin(o(i,j,k)+h(i,j,k)*sin(w(i,j,k); %经纬高转换为直角坐标系xxe(i,j,k)=xe;yye(i,j,k)=ye;zze(i,j,k)=ze;endendend-%直角坐标系 xyz 转换为经纬高 wlhfunction NEWw,NEWl,NEWh,flag_converge,Nkk=xyz_change_wlh(xxe,yye,zze,kkmax,eps)if(nargin=3)
9、eps=1.0e-9; % %缺省误差参数时默认为 eps=1.0e-9kkmax=20; %最大迭代次数endR=6378.173*1000;r=6356.7523142*1000;f=(R-r)/R;NEWw=; %存放纬度的数组NEWl=; %存放经度的数组NEWh=; %存放高度的数组flag_converge=;%存放发散标志的数组Nkk=;for i=1:length(xxe(:,1,1) %纬度 ifor j=1:length(xxe(1,:,1) %经度 jfor k=1:length(xxe(1,1,:) %高度 kx0=sqrt(xxe(i,j,k)2+yye(i,j,k)
10、2);%根据直角坐标系求 x0.z0z0=zze(i,j,k); %step2A=r*z0; %系数 A，在求牛顿迭代算法时使用B=2*R*x0;C=2*(R2-r2);if z0=0 %step3 z0=0 的情况，neww=0; %求纬度newh=x0-R; %求高度newl=atan2(yye(i,j,k),xxe(i,j,k); %step10 ，求经度elsekk=0; %迭代次数初值t=-(1-f)*(x0/z0)+sign(z0)*sqrt(1+(1-f)*(x0/z0)2); %step4 ，t0 初值flag=1; %标志位，判断是否收敛while flag=1 kk=kk+
11、1; %step5 迭代次数加 1ff=A*t4+(B+C)*t3+(B-C)*t-A;df=4*A*t3+3*(B+C)*t2+(B-C);tnew=t-ff/df; %step6 t 值的更新oldt=t; %存储旧的 tt=tnew; %将新的 tnew 值赋给 tif abs(t-oldt)=kkmax %step8，若小于，则到 step5（通过标志位不改变来实现）flag=0; %若大于，则跳出循环endendendendNEWw(i,j,k)=neww;NEWl(i,j,k)=newl; NEWh(i,j,k)=newh; flag_converge(i,j,k)=flag;Nkk(i,j,k)=kk;endendendNEWw=NEWw*(180/pi); %纬度弧度转换为角度NEWl=NEWl*(180/pi);%经度弧度转换为角度
