高斯投影(高斯-克吕格投影)的反算
更新2020-06，将坐标系统统一换为WGS-84坐标系，整理一下脚本函数
高斯投影的反算是指由当地的局部坐标系(x,y)转换为当地的地理坐标系(B: 纬度, L: 经度)。由于之前的博文MATLAB实现高斯-克吕格投影正算已经对高斯投影进行过简要的说明，故本博文不再对高斯-克吕格投影的原理进行介绍，只给出高斯投影反算的算法流程和实现的MATLAB脚本。本博文参考文献资料如下：
[1] 孔祥元, 郭际明, 刘宗泉, 等. 大地测量学基础(第二版)[M]. 武汉: 武汉大学出版社, 2009.
[2] 程鹏飞,成英燕,文汉江,等.2000国家大地坐标系实用宝典[M]. 北京:测绘出版社, 2008.
[3] 牛丽娟. 测量坐标转换模型研究与转换系统实现[D]. 长安大学, 2010.
[4] 李厚朴,边少锋.高斯投影的复变函数表示[J].测绘学报,2008(01):5-9
高斯投影反算公式
当地纬度 B B B的计算公式如下：
B = B f − ρ t f 2 M f y ( y N f ) [ 1 − 1 12 ( 5 + 3 t f 2 + η f 2 − 9 η f 2 t f 2 ) ( y N f ) 2 + 1 360 ( 61 + 90 t f 2 + 45 t f 4 ) ( y N f ) 4 ] \begin{aligned} B=& B_{f}-\frac{\rho t_{f}}{2 M_{f}} y\left(\frac{y}{N_{f}}\right)\left[1-\frac{1}{12}\left(5+3 t_{f}^{2}+\eta_{f}^{2}-9 \eta_{f}^{2} t_{f}^{2}\right)\left(\frac{y}{N_{f}}\right)^{2}\right.\\ &+\left.\frac{1}{360}\left(61+90 t_{f}^{2}+45 t_{f}^{4}\right)\left(\frac{y}{N_{f}}\right)^{4}\right] \end{aligned} B=Bf−2Mfρtfy(Nfy)[1−121(5+3tf2+ηf2−9ηf2tf2)(Nfy)2+3601(61+90tf2+45tf4)(Nfy)4]
当 y = 0 y=0 y=0时的纬度称为底点纬度 B f B_f Bf，底点纬度的迭代计算公式为：
X
=
a
0
B
−
a
2
2
sin
2
B
+
a
4
4
sin
4
B
−
a
6
6
sin
6
B
+
a
8
8
sin
8
B
X=a_{0} B-\frac{a_{2}}{2} \sin 2 B+\frac{a_{4}}{4} \sin 4 B-\frac{a_{6}}{6} \sin 6 B+\frac{a_{8}}{8} \sin 8 B
X=a0B−2a2sin2B+4a4sin4B−6a6sin6B+8a8sin8B其中，
{
a
0
=
m
0
+
1
2
m
2
+
3
8
m
4
+
5
16
m
6
+
35
128
m
8
a
2
=
1
2
m
2
+
1
2
m
4
+
15
32
m
6
+
7
16
m
8
a
4
=
1
8
m
4
+
3
16
m
6
+
7
32
m
8
a
6
=
1
32
m
6
+
1
16
m
8
a
8
=
1
128
m
8
;
{
m
0
=
a
(
1
−
e
2
)
m
2
=
3
2
e
2
m
0
m
4
=
5
4
e
2
m
2
m
6
=
7
6
e
2
m
4
m
8
=
9
8
e
2
m
6
\left\{\begin{aligned} a_{0}=m_{0}+\frac{1}{2} m_{2}+\frac{3}{8} m_{4}+\frac{5}{16} m_{6}+\frac{35}{128} m_{8} \\ a_{2}=\frac{1}{2} m_{2}+\frac{1}{2} m_{4}+\frac{15}{32} m_{6}+\frac{7}{16} m_{8} \\ a_{4}=\frac{1}{8} m_{4}+\frac{3}{16} m_{6}+\frac{7}{32} m_{8} \\ a_{6}=\frac{1}{32} m_{6}+\frac{1}{16} m_{8} \\ a_{8}=\frac{1}{128} m_{8} \end{aligned};\quad \right. \left\{\begin{aligned} m_{0}=a(1-e^2) \\ m_{2}=\frac{3}{2} e^2 m_{0} \\ m_{4}=\frac{5}{4} e^2 m_{2}\\ m_{6}=\frac{7}{6} e^2 m_{4} \\ m_{8}=\frac{9}{8} e^2 m_{6} \end{aligned}\right.
⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧a0=m0+21m2+83m4+165m6+12835m8a2=21m2+21m4+3215m6+167m8a4=81m4+163m6+327m8a6=321m6+161m8a8=1281m8;⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧m0=a(1−e2)m2=23e2m0m4=45e2m2m6=67e2m4m8=89e2m6
B
f
B_f
Bf迭代计算过程为：
迭代初始: 令
B
f
0
=
x
a
0
B_f^0 = \frac{x}{a_0}
Bf0=a0x
迭代过程: 令
F
i
=
−
a
2
2
sin
2
B
f
i
+
a
4
4
sin
4
B
f
i
−
a
6
6
sin
6
B
f
i
+
a
8
8
sin
8
B
f
i
F^{i}=-\frac{a_{2}}{2} \sin 2 B_{f}^{i}+\frac{a_{4}}{4} \sin 4 B_{f}^{i}-\frac{a_{6}}{6} \sin 6 B_{f}^{i}+\frac{a_{8}}{8} \sin 8 B_{f}^{i}
Fi=−2a2sin2Bfi+4a4sin4Bfi−6a6sin6Bfi+8a8sin8Bfi; 令
B
f
i
+
1
=
X
−
F
i
a
0
B_f^{i+1} = \frac{X-F^i}{a_0}
Bfi+1=a0X−Fi
迭代停止: 当
∣
B
f
i
+
1
−
B
f
i
∣
<
ϵ
|B_f^{i+1} - B_f^i|<\epsilon
∣Bfi+1−Bfi∣<ϵ时则停止迭代.
注：
ϵ
\epsilon
ϵ为给定的阈值，阈值通常小于1E-8.
当地经度 l l l的计算公式如下:
l = ρ cos B f ( y N f ) [ 1 − 1 6 ( 1 + 2 t f 2 + η f 2 ) ( y N f ) 2 + 1 120 ( 5 + 28 t f 2 + 24 t f 4 + 6 η f 2 + 8 η f 2 t f 2 ) ( y N f ) 4 ] \begin{aligned} l=& \frac{\rho}{\cos B_{f}}\left(\frac{y}{N_{f}}\right)\left[1-\frac{1}{6}\left(1+2 t_{f}^{2}+\eta_{f}^{2}\right)\left(\frac{y}{N_{f}}\right)^{2}+\frac{1}{120}\left(5+28 t_{f}^{2}\right.\right.\\ &+\left.\left.24 t_{f}^{4}+6 \eta_{f}^{2}+8 \eta_{f}^{2} t_{f}^{2}\right)\left(\frac{y}{N_{f}}\right)^{4}\right] \end{aligned} l=cosBfρ(Nfy)[1−61(1+2tf2+ηf2)(Nfy)2+1201(5+28tf2+24tf4+6ηf2+8ηf2tf2)(Nfy)4]
上述式中:
- x x x为北方向的坐标，
- y y y为东方向的坐标(中国地区内需减去500000m)，
- ρ = 180 × 3600 / π \rho=180×3600/\pi ρ=180×3600/π为弧度秒，
- a a a为地球椭球长半径，
- e e e为地球椭球第一偏心率，
- e ′ e' e′为地球椭球第二偏心率，
- N f N_f Nf为底点纬度 B f B_f Bf处地球的卯酉圈曲率半径，
- M f M_f Mf为底点纬度 B f B_f Bf处地球的子午圈曲率半径，
- η f 2 = e ′ 2 c o s 2 B f \eta_f^2 = e'^2cos^2{B_f} ηf2=e′2cos2Bf，
- t f = t a n B f t_f=tan{B_f} tf=tanBf。
高斯投影反算的MATLAB函数
function Coord = GaussProInverse(Pos)
% INPUT // Units of longitude and latitude is RAD (°)
% REF[1]// 程鹏飞,等.2000国家大地坐标系实用宝典[M].北京:测绘出版社,2008.
% REF[2]// 孔祥元,郭际明.大地测量学基础(第二版)[M].武汉:武汉大学出版社,2010.
% Longitude of the Earth central meridian
MerLon = 114; % Wuhan is 114 deg
% Extract the local projected coordinate (x & y)
Coord = Pos;
Eth.D2R = 0.0174532925199433; % pi/180
x = Pos(1);
y = Pos(2) - 500000;
%% Earth orientation parameters of WGS 84 Coordinate System
Eth.R0 = 6378137.0;
Eth.f = 1/298.257223563;
Eth.Rp = 6356752.3142452; % R0*(1 - f);
% First Eccentricity and its Squared
Eth.e12 = 0.006694379990141; % (2f - f*f);
Eth.e11 = 0.081819190842622; % sqrt(2f - f*f);
% Second Eccentricity and its Squared
Eth.e22 = 0.00673949674227643; % (2f - f*f)/(1 + f*f - 2*f);
Eth.e21 = 0.08209443794969570; % sqrt((2f - f*f)/(1 + f*f - 2*f));
Bf = Meridian2Latitude(x, Eth.R0, Eth.e12);
tnBf = tan(Bf); tn2Bf = tnBf*tnBf; tn4Bf = tn2Bf*tn2Bf;
csBf = cos(Bf); cs2Bf = csBf*csBf; Eta2 = Eth.e22*cs2Bf;
COE = sqrt(1 - Eth.e12*sin(Bf)*sin(Bf));
Nf = Eth.R0/COE;
Mf = Eth.R0*(1 - Eth.e12)/COE^3;
%% Calculate Latitude(B)
YNf = y/Nf;
YNf2 = YNf*YNf;
YNf4 = YNf2*YNf2;
TYMN = 0.5*tnBf*y*YNf/Mf;
COE1 = (5 + 3*tn2Bf + Eta2 -9*Eta2*tn2Bf)*YNf2/12;
COE2 = (61 + 90*tn2Bf + 45*tn4Bf)*YNf4/360;
Lat = Bf - TYMN*(1 - COE1 + COE2);
%% Calculate Longitude(L)
YBNf = YNf/csBf;
COE3 = (1 + 2*tn2Bf + Eta2)*YNf2/6;
COE4 = (5 + 28*tn2Bf + 24*tn4Bf + 6*Eta2 + 8*Eta2*tn2Bf)*YNf4/120;
Lon = MerLon*Eth.D2R + YBNf*(1 - COE3 + COE4);
Coord(1) = Lat/Eth.D2R;
Coord(2) = Lon/Eth.D2R;
end
function Bf = Meridian2Latitude(x,a,e12)
m0 = a*(1 - e12); m2 = 3*e12*m0/2;
m4 = 5*e12*m2/4; m6 = 7*e12*m4/6;
m8 = 9*e12*m6/8; a8 = m8/128;
a6 = m6/32 + m8/16; a4 = m4/8 + 3*m6/16 + 7*m8/32;
a0 = m0 + m2/2 + 3*m4/8 + 5*m6/16 + 35*m8/128;
a2 = m2/2 + m4/2 + 15*m6/32 + 7*m8/16;
B0 = x/a0;
while 1
F = -a2*sin(2*B0)/2 + a4*sin(4*B0)/4 - a6*sin(6*B0)/6 + a8*sin(8*B0)/8;
Bf = (x - F)/a0;
if abs(B0 - Bf)<1e-10
break;
end
B0 = Bf;
end
end
扫一扫
3563
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
