MATLAB实现高斯-克吕格投影反算
高斯投影(高斯-克呂格投影)的反算
更新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
高斯投影反算公式
當地緯度BBB的計算公式如下:
B=Bf?ρtf2Mfy(yNf)[1?112(5+3tf2+ηf2?9ηf2tf2)(yNf)2+1360(61+90tf2+45tf4)(yNf)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?ρtf??y(Nf?y?)[1?121?(5+3tf2?+ηf2??9ηf2?tf2?)(Nf?y?)2+3601?(61+90tf2?+45tf4?)(Nf?y?)4]?
當y=0y=0y=0時的緯度稱為底點緯度BfB_fBf?,底點緯度的迭代計算公式為:
X=a0B?a22sin?2B+a44sin?4B?a66sin?6B+a88sin?8BX=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 BX=a0?B?2a2??sin2B+4a4??sin4B?6a6??sin6B+8a8??sin8B其中,{a0=m0+12m2+38m4+516m6+35128m8a2=12m2+12m4+1532m6+716m8a4=18m4+316m6+732m8a6=132m6+116m8a8=1128m8;{m0=a(1?e2)m2=32e2m0m4=54e2m2m6=76e2m4m8=98e2m6\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?+21?m2?+83?m4?+165?m6?+12835?m8?a2?=21?m2?+21?m4?+3215?m6?+167?m8?a4?=81?m4?+163?m6?+327?m8?a6?=321?m6?+161?m8?a8?=1281?m8??;????????????????????????????m0?=a(1?e2)m2?=23?e2m0?m4?=45?e2m2?m6?=67?e2m4?m8?=89?e2m6??
BfB_fBf?迭代計算過程為:
迭代初始: 令Bf0=xa0B_f^0 = \frac{x}{a_0}Bf0?=a0?x?
迭代過程: 令Fi=?a22sin?2Bfi+a44sin?4Bfi?a66sin?6Bfi+a88sin?8BfiF^{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=?2a2??sin2Bfi?+4a4??sin4Bfi??6a6??sin6Bfi?+8a8??sin8Bfi?; 令Bfi+1=X?Fia0B_f^{i+1} = \frac{X-F^i}{a_0}Bfi+1?=a0?X?Fi?
迭代停止: 當∣Bfi+1?Bfi∣<?|B_f^{i+1} - B_f^i|<\epsilon∣Bfi+1??Bfi?∣<?時則停止迭代.
注:?\epsilon?為給定的閾值,閾值通常小于1E-8.
當地經度lll的計算公式如下:
l=ρcos?Bf(yNf)[1?16(1+2tf2+ηf2)(yNf)2+1120(5+28tf2+24tf4+6ηf2+8ηf2tf2)(yNf)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?ρ?(Nf?y?)[1?61?(1+2tf2?+ηf2?)(Nf?y?)2+1201?(5+28tf2?+24tf4?+6ηf2?+8ηf2?tf2?)(Nf?y?)4]?
上述式中:
高斯投影反算的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; endfunction 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 1F = -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-10break;endB0 = Bf; end end該方法為經典的高斯反算的算法,對于高斯反算的詳細推導可以參考大地測量的教材或者最新的文獻。
總結
以上是生活随笔為你收集整理的MATLAB实现高斯-克吕格投影反算的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 啥叫旁路电容?啥叫去耦?可以不再争论了吗
- 下一篇: 听说有人不了解柔性数组