大地测量学小脚本分享

大地测量学中的一些简单变换

经过两个多月的学习，对大地测量学的知识有了一定的了解，课堂上老师也布置了一些简单的任务，在这里把自己的一些成果和经验和大家分享一下，希望对将要接触或正在学习本学科的同学有所帮助。
在这节课的第一个任务就是子午线于卯酉线曲率半径的描绘。MATLAB代码如下：

a=6378137;
b=6356752.3142;
e=sqrt((a^2-b^2)/a^2);
a1=0:0.01:90;
a1=pi/180*a1;
W=sqrt(1-e^2*sin(a1).^2);
M=a*(1-e^2)./W.^3;
N=a./W;
R=sqrt(M.*N);
plot(a1,M,a1,N,a1,R);
xlabel('纬度');
ylabel('长度');
title('曲率半径');
legend('子午线曲率半径','卯酉线曲率半径','平均曲率半径');

得到如下图象：

之后又学习了高斯变换，通过大地坐标（B，L）->高斯投影(X，Y)的变换，MATLAB代码如下：

pi=3.1415926535898;
B=-pi/2:0.01:pi/2;
l=-pi/60:0.01:pi/60;
[B,l]=meshgrid(B,l);
a=6378137;
b=6356752.3142;
e1=sqrt((a^2-b^2)/a^2);
e2=sqrt((a^2-b^2)/b^2);
t=tan(B);
n=sqrt(e2^2.*cos(B).^2);
W=sqrt(1-e1^2*sin(B).^2);
M=a*(1-e1^2)./W.^3;
N=a./W;
m0=a*(1-e1^2);
m2=3/2*e1^2*m0;
m4=5/4*e1^2*m2;
m6=7/6*e1^2*m4;
m8=9/8*e1^2*m6;
a0=m0+m2/2+3/8*m4+5/16*m6+35/128*m8;
a2=m2/2+m4/2+15/32*m6+7/16*m8;
a4=m4/8+3/16*m6+7/32*m8;
a6=m6/32+m8/16;
a8=m8/128;
%克拉索夫斯基
% X=111134.861.*B/pi*180-16036.480.*sin(2*B)+16.828.*sin(4*B)-0.022.*sin(6*B);
X=a0*B-a2/2*sin(2*B)+a4/4*sin(4*B)-a6/6*sin(6*B)+a8/8.*sin(8*B);
x=X+N./2.*sin(B).*cos(B).*l.^2+N./24.*sin(B).*cos(B).^3.*(5-t.^2+9*n.^2+4*n.^4).*l.^4+N/720.*sin(B).*cos(B).^5.*(61-58*t.^2+t.^4).*l.^6;
y=N.*cos(B).*l+N./6.*cos(B).^3.*(1-t.^2+n.^2).*l.^3+N./120.*cos(B).^5.*(5-18*t.^2+t.^4+14*n.^2-58*n.^2.*t.^2).*l.^5;
subplot(3,1,1);mesh(l,B,x)
xlabel('经差');ylabel('纬度');zlabel('x坐标');
subplot(3,1,2);mesh(l,B,y)
xlabel('经差');ylabel('纬度');zlabel('y坐标');
subplot(3,1,3);
plot(y,x);
xlabel('y坐标');ylabel('x坐标');

得到如下的图像：


可以发现高斯投影的部分规律。
再后来又学习了高斯反算，MATLAB代码如下：

clc;
clear;
x=6.655021690543837e+06;
y=8.369522497201736e+04;
pi=3.1415926535898;
a=6378137;
b=6356752.3142;
e1=sqrt((a^2-b^2)/a^2);
e2=sqrt((a^2-b^2)/b^2);
m0=a*(1-e1^2);
m2=3/2*e1^2*m0;
m4=5/4*e1^2*m2;
m6=7/6*e1^2*m4;
m8=9/8*e1^2*m6;
a0=m0+m2/2+3/8*m4+5/16*m6+35/128*m8;
a2=m2/2+m4/2+15/32*m6+7/16*m8;
a4=m4/8+3/16*m6+7/32*m8;
a6=m6/32+m8/16;
a8=m8/128;
B0=x/a0;
%B1=(x-(-a2/2*sin(2*B0)+a4/4*sin(4*B0)-a6/6*sin(6*B0)+a8/8.*sin(8*B0)))/a0;
for i=1:2000
    B1=(x-(-a2/2*sin(2*B0)+a4/4*sin(4*B0)-a6/6*sin(6*B0)+a8/8.*sin(8*B0)))/a0;
    if(abs(B1-B0)<10e-10)
        break;
    end
    B0=B1;
end
t0=tan(B0);
n0=sqrt(e2^2.*cos(B0).^2);
W0=sqrt(1-e1^2*sin(B0).^2);
M0=a*(1-e1^2)./W0.^3;
N0=a./W0;
B=B0-t0/(2*M0*N0)*y^2+t0/(24*M0*N0^3)*(5+3*t0^2+n0^2-9*n0^2*t0^2)*y^4-t0/(720*M0*N0^5)*(61+90*t0^2+45*t0^4)*y^6;
L=1/(N0*cos(B0))*y-1/(6*N0^3*cos(B0))*(1+2*t0^2+n0^2)*y^3+1/(120*N0^5*cos(B0))*(5+28*t0^2+24*t0^4+6*n0^2+8*n0^2*t0^2)*y^5;
B= radtoganle(B)
l=radtoganle(L)

其中包含radtoganle函数是自己设置的一个函数，方便观察结果，内容如下：

function angle = radtoganle(rad)
d=fix(rad*180/pi);
f=fix((rad*180/pi-d)*60);
m=((rad*180/pi-d)*60-f)*60;
%fprintf('%f',(((rad*180/pi-d)*60-f)*60)-m);
%if ((rad*180/pi-d)*60-f)*60-fix(m)>0.98
%    m=m+1;
%end
if m>=60
    m=m-60;
    f=f+1;
end
if f>=60
    f=f-60;
    d=d+1;
end
angle=[d,f,m]; 

最后还做了一个测试函数来检验函数的正确性。

pi=3.1415926535898;
B=pi/3;
l=pi/120;
radtoganle(B)
radtoganle(l)
a=6378137;
b=6356752.3142;
e1=sqrt((a^2-b^2)/a^2);
e2=sqrt((a^2-b^2)/b^2);
t=tan(B);
n=sqrt(e2^2.*cos(B).^2);
W=sqrt(1-e1^2*sin(B).^2);
M=a*(1-e1^2)./W.^3;
N=a./W;
m0=a*(1-e1^2);
m2=3/2*e1^2*m0;
m4=5/4*e1^2*m2;
m6=7/6*e1^2*m4;
m8=9/8*e1^2*m6;
a0=m0+m2/2+3/8*m4+5/16*m6+35/128*m8;
a2=m2/2+m4/2+15/32*m6+7/16*m8;
a4=m4/8+3/16*m6+7/32*m8;
a6=m6/32+m8/16;
a8=m8/128;
%克拉索夫斯基
% X=111134.861.*B/pi*180-16036.480.*sin(2*B)+16.828.*sin(4*B)-0.022.*sin(6*B);
X=a0*B-a2/2*sin(2*B)+a4/4*sin(4*B)-a6/6*sin(6*B)+a8/8.*sin(8*B);
x=X+N./2.*sin(B).*cos(B).*l.^2+N./24.*sin(B).*cos(B).^3.*(5-t.^2+9*n.^2+4*n.^4).*l.^4+N/720.*sin(B).*cos(B).^5.*(61-58*t.^2+t.^4).*l.^6;
y=N.*cos(B).*l+N./6.*cos(B).^3.*(1-t.^2+n.^2).*l.^3+N./120.*cos(B).^5.*(5-18*t.^2+t.^4+14*n.^2-58*n.^2.*t.^2).*l.^5;
B0=x/a0;
%B1=(x-(-a2/2*sin(2*B0)+a4/4*sin(4*B0)-a6/6*sin(6*B0)+a8/8.*sin(8*B0)))/a0;
for i=1:2000
    B1=(x-(-a2/2*sin(2*B0)+a4/4*sin(4*B0)-a6/6*sin(6*B0)+a8/8.*sin(8*B0)))/a0;
    if(abs(B1-B0)<10e-10)
        break;
    end
    B0=B1;
end
t0=tan(B0);
n0=sqrt(e2^2.*cos(B0).^2);
W0=sqrt(1-e1^2*sin(B0).^2);
M0=a*(1-e1^2)./W0.^3;
N0=a./W0;
B=B0-t0/(2*M0*N0)*y^2+t0/(24*M0*N0^3)*(5+3*t0^2+n0^2-9*n0^2*t0^2)*y^4-t0/(720*M0*N0^5)*(61+90*t0^2+45*t0^4)*y^6;
l=1/(N0*cos(B0))*y-1/(6*N0^3*cos(B0))*(1+2*t0^2+n0^2)*y^3+1/(120*N0^5*cos(B0))*(5+28*t0^2+24*t0^4+6*n0^2+8*n0^2*t0^2)*y^5;
%rad2deg(B)
%rad2deg(l)
B= radtoganle(B)
l=radtoganle(l)

不过在纬度到达九十度的结果会和应有结果相差较大，会有两秒左右的差值，这个问题目前还未解决，由于测试集过小，目前函数可能还有其他漏洞，希望有兴趣的小伙伴或是大佬能帮忙解决。