MATLAB实现高斯-克吕格投影正算

MATLAB实现高斯-克吕格投影正算-即经纬度转为x和y

高斯-克吕格投影简介

更新2020-06，重新整理一下脚本函数

高斯-克吕格投影，是由德国数学家、物理学家、天文学家高斯于1822年代首次提出，后经德国大地测量学家克吕格于1912年对投影公式加以补充，故称高斯-克吕格投影，又名"等角横切椭圆柱投影”，是地球椭球面和平面间正形投影的一种。该投影是用一个设想的圆柱筒横置于地球表面，与地球相切于某一经线（中央经线），圆柱的中心轴位于赤道面内，按等角条件将地球椭球面投影于椭球圆柱面上。
为控制投影变形，先按一定的经度差（6°或者3°）将地球表面划分为若干投影带，再使圆柱面依次和每一带的中央经线相切，并把各带中央经线东西两侧一定经度差范围内的经纬线网投影到圆柱上，然后从两级将该圆柱面切开展平，构成地球各带经纬线网在平面上的图形。

高斯-克吕格投影正算公式

x = X + N t c o s 2 B l 2 ρ 2 [ 0.5 + 1 24 ( 5 − t 2 + 9 η 2 + 4 η 4 ) c o s 2 B l 2 ρ 2 + 1 720 ( 61 − 58 t 2 + t 4 ) c o s 4 B l 4 ρ 4 ] ( 1 ) y = N c o s B l ρ [ 1 + 1 6 ( 1 − t 2 + η 2 ) c o s 2 B l 2 ρ 2 + 1 120 ( 5 − 18 t 2 + t 4 + 14 η 2 − 58 t 2 η 2 ) c o s 4 B l 4 ρ 4 ] ( 2 ) \begin{alignedat}{4} x=X+Ntcos^2{B}\frac{l^2}{\rho^2} \lbrack 0.5+ \frac{1}{24}(5-t^2+9\eta^2+4\eta^4)cos^2{B}\frac{l^2}{\rho^2} +\frac{1}{720} (61-58t^2+t^4) cos^4{B}\frac{l^4}{\rho^4}\rbrack \qquad (1) \\ y=Ncos{B}\frac{l}{\rho}[1+\frac{1}{6}(1-t^2+\eta^2)cos^2{B}\frac{l^2}{\rho^2} + \frac{1}{120}(5-18t^2+t^4+14\eta^2-58t^2\eta^2)cos^4{B}\frac{l^4}{\rho^4}] \qquad (2) \end{alignedat} x=X+Ntcos2Bρ2l2[0.5+241(5−t2+9η2+4η4)cos2Bρ2l2+7201(61−58t2+t4)cos4Bρ4l4](1)y=NcosBρl[1+61(1−t2+η2)cos2Bρ2l2+1201(5−18t2+t4+14η2−58t2η2)cos4Bρ4l4](2)其中， X X X为中央子午线弧长， N N N为卯酉圈曲率半径， t = t a n B t=tan{B} t=tanB， ρ = 180 × 3600 / π \rho=180×3600/\pi ρ=180×3600/π为弧度秒， η 2 = e ′ 2 c o s 2 B \eta^2 = e'^2cos^2{B} η2=e′2cos2B， e ′ e' e′为地球椭球第二偏心率， B B B为当地纬度。卯酉圈曲率半径及中央子午线弧长公式如下： N = a 1 − e 2 s i n 2 B ( 3 ) N=\frac{a}{\sqrt{1-e^2sin^2{B}}} \qquad (3) N=1−e2sin2B a(3) X = a ( 1 − e 2 ) ( A ′ a r c B − B ′ s i n 2 B + C ′ s i n 4 B − D ′ s i n 6 B + E ′ s i n 8 B − F ′ s i n 10 B + G ′ s i n 12 B ) ( 4 ) X =a(1-e^2)(A'arcB-B'sin{2B}+C'sin{4B}-D'sin{6B}\\ +E'sin{8B}-F'sin{10B}+G'sin{12B}) \qquad (4) X=a(1−e2)(A′arcB−B′sin2B+C′sin4B−D′sin6B+E′sin8B−F′sin10B+G′sin12B)(4)子午线弧长计算公式中的各项符号具体公式如下： A ′ = 1 + 3 4 e 2 + 45 64 e 4 + 175 256 e 6 + 11025 16384 e 8 + 43659 65536 e 10 + 693693 1048576 e 12 A'=1+\frac{3}{4}e^2+\frac{45}{64}e^4+\frac{175}{256}e^6+\frac{11025}{16384}e^8+\frac{43659}{65536}e^{10}+\frac{693693}{1048576}e^{12} A′=1+43e2+6445e4+256175e6+1638411025e8+6553643659e10+1048576693693e12 B ′ = 3 8 e 2 + 15 32 e 4 + 525 1024 e 6 + 2205 4096 e 8 + 72765 131072 e 10 + 297297 524288 e 12 B'=\frac{3}{8}e^2+\frac{15}{32}e^4+\frac{525}{1024}e^6+\frac{2205}{4096}e^8+\frac{72765}{131072}e^{10}+\frac{297297}{524288}e^{12} B′=83e2+3215e4+1024525e6+40962205e8+13107272765e10+524288297297e12 C ′ = 15 256 e 4 + 105 1024 e 6 + 2205 16384 e 8 + 10396 65536 e 10 + 1486485 8388608 e 12 C'=\frac{15}{256}e^4+\frac{105}{1024}e^6+\frac{2205}{16384}e^8+\frac{10396}{65536}e^{10}+\frac{1486485}{8388608}e^{12} C′=25615e4+1024105e6+163842205e8+6553610396e10+83886081486485e12 D ′ = 35 3072 e 6 + 105 4096 e 8 + 10395 262144 e 10 + 55055 1048576 e 12 D'=\frac{35}{3072}e^6+\frac{105}{4096}e^8+\frac{10395}{262144}e^{10}+\frac{55055}{1048576}e^{12} D′=307235e6+4096105e8+26214410395e10+104857655055e12 E ′ = 315 131072 e 8 + 3465 524288 e 10 + 99099 8388608 e 12 E'=\frac{315}{131072}e^8+\frac{3465}{524288}e^{10}+\frac{99099}{8388608}e^{12} E′=131072315e8+5242883465e10+838860899099e12 F ′ = 693 1310720 e 10 + 9909 5242880 e 12 F'=\frac{693}{1310720}e^{10}+\frac{9909}{5242880}e^{12} F′=1310720693e10+52428809909e12 G ′ = 1001 8388608 e 12 G'=\frac{1001}{8388608}e^{12} G′=83886081001e12其中， e e e为地球椭球第一偏心率， a a a为地球椭球长半轴。
注意 1 ：当纬度已经为弧度 ( π = 180 ° ) 时 ρ = 1 ; \color{blue}注意1：当纬度已经为弧度(\pi=180°)时\rho=1; 注意1：当纬度已经为弧度(π=180°)时ρ=1;
注意 2 ：中国地区内，为了避免出现负数， y 需要加上 500000 m . \color{blue}注意2：中国地区内，为了避免出现负数，y需要加上500000m. 注意2：中国地区内，为了避免出现负数，y需要加上500000m.
在WGS-84坐标系中， a = 6378137 m ， f = 1 / 298.257223563 ， e 2 = 2 f − f 2 ， e ′ 2 = ( 2 f − f 2 ) / ( 1 − f 2 ) a=6378137m，f = 1/298.257223563，e^2=2f-f^2，e'^2={(2f-f^2)}/{(1-f^2)} a=6378137m，f=1/298.257223563，e2=2f−f2，e′2=(2f−f2)/(1−f2)

高斯-克吕格投影的MATLAB函数

function Pos = GaussProDirect(Coord)
% 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;
% Extract the longitude and latitude
Pos = Coord;
Eth.D2R = 0.0174532925199433;  % pi/180
Lat = Coord(1)*Eth.D2R;
Lon = Coord(2)*Eth.D2R - MerLon*Eth.D2R;
%% 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));
%% 高斯投影正算公式
RN = Eth.R0/sqrt(1 - Eth.e12*sin(Lat)*sin(Lat));
Lon2 = Lon*Lon;   Lon4 = Lon2*Lon2;
tnLat = tan(Lat); tn2Lat = tnLat*tnLat; tn4Lat = tn2Lat*tn2Lat;
csLat = cos(Lat); cs2Lat = csLat*csLat; cs4Lat = cs2Lat*cs2Lat;
Eta2 = Eth.e22*cs2Lat;   Eta4 = Eta2*Eta2;
NTBLP = RN*tnLat*cs2Lat*Lon2;
COE1 = (5 - tn2Lat + 9*Eta2 + 4*Eta4)*cs2Lat*Lon2/24;
COE2 = (61 - 58*tn2Lat + tn4Lat)*cs4Lat*Lon4/720;
x = Merdian(Eth,Lat) + NTBLP*(0.5 + COE1 + COE2);
NBLP = RN*csLat*Lon;
COE3 = (1 - tn2Lat + Eta2)*cs2Lat*Lon2/6;
COE4 = (5 - 18*tn2Lat + tn4Lat + 14*Eta2 - 58*tn2Lat*Eta2)*cs4Lat*Lon4/120;
y = NBLP*(1 + COE3 + COE4) + 500000;
Pos(1) = x; Pos(2) = y;
endfunction X0 = Merdian(Eth,Lat)
% REF//过家春.子午线弧长公式的简化及其泰勒级数解释[J].测绘学报,2014,43(2):125-130.
S0 = Eth.R0*(1 - Eth.e12);
e2 = Eth.e12;
e4 = e2*e2;
e6 = e4*e2;
e8 = e6*e2;
e10 = e8*e2;
e12 = e10*e2;
A1 = 1 + 3*e2/4 + 45*e4/64 + 175*e6/256 + 11025*e8/16384 + 43659*e10/65536 + 693693*e12/1048576;
B1 = 3*e2/8 + 15*e4/32 + 525*e6/1024 + 2205*e8/4096 + 72765*e10/131072 + 297297*e12/524288;
C1 = 15*e4/256 + 105*e6/1024 + 2205*e8/16384 + 10395*e10/65536 + 1486485*e12/8388608;
D1 = 35*e6/3072 + 105*e8/4096 + 10395*e10/262144 + 55055*e12/1048576;
E1 = 315*e8/131072 + 3465*e10/524288 + 99099*e12/8388608;
F1 = 693*e10/1310720 + 9009*e12/5242880;
G1 = 1001*e12/8388608;
X0 = S0*(A1*Lat - B1*sin(2*Lat) + C1*sin(4*Lat) - D1*sin(6*Lat) +...E1*sin(8*Lat) - F1*sin(10*Lat) + G1*sin(12*Lat));
end

函数代码使用示例

示例中纬度为30.4691868227°，经度为114.3510760836°，中央子午线经度为114°(武汉)，参考真值为：x(北向) = 3372178.140m，y(东向) = 533713.649m。计算结果如下。可以看出：计算结果正确。

参考资料

[1] 程鹏飞,成英燕,文汉江,等.2000国家大地坐标系实用宝典[M]. 北京:测绘出版社,2008.
[2] 孔祥元, 郭际明. 大地测量学基础[M]. 武汉: 武汉大学出版社, 2005.
[3] 过家春.子午线弧长公式的简化及其泰勒级数解释[J].测绘学报,2014,43(2):125-130.

完毕