MATLAB实现高斯-克吕格投影正算
MATLAB实现高斯-克吕格投影正算
高斯-克吕格投影简介
高斯-克吕格投影，是由德国数学家、物理学家、天文学家高斯于18代首次提出，后经德国大地测量学家克吕格于19对投影公式加以补充，故称高斯-克吕格投影，又名"等角横切椭圆柱投影”，是地球椭球面和平面间正形投影的一种。该投影是用一个设想的圆柱筒横置于地球表面，与地球相切于某一经线(中央经线)，圆柱的中心轴位于赤道面内，按等角条件将地球椭球面投影于椭球圆柱面上。
为控制投影变形，先按一定的经度差(6°或者3°)将地球表面划分为若干投影带，再使圆柱面依次和每一带的中央经线相切，并把各带中央经线东西两侧一定经度差范围内的经纬线网投影到圆柱上，然后从两级将该圆柱面切开展平，构成地球各带经纬线网在平面上的图形。
高斯-克吕格投影正算公式
x=X+Ntcos2Bl2ρ2[0.5+124(5?t2+9η2+4η4)cos2Bl2ρ2+1720(61?58t2+t4)cos4Bl4ρ4](1)y=NcosBlρ[1+16(1?t2+η2)cos2Bl2ρ2+1120(5?18t2+t4+14η2?58t2η2)cos4Bl4ρ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)?其中，XXX为中央子午线弧长，NNN为卯酉圈曲率半径，t=tanBt=tan{B}t=tanB，ρ=180×3600/π\rho=180×3600/\piρ=180×3600/π为弧度秒，η2=e′2cos2B\eta^2 = e'^2cos^2{B}η2=e′2cos2B，e′e'e′为地球椭球第二偏心率，BBB为当地纬度。卯酉圈曲率半径及中央子午线弧长公式如下：N=a1?e2sin2B(3) N=\frac{a}{\sqrt{1-e^2sin^2{B}}} \qquad (3)N=1?e2sin2B?a?(3) X=a(1?e2)(A′arcB?B′sin2B+C′sin4B?D′sin6B+E′sin8B?F′sin10B+G′sin12B)(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+34e2+4564e4+175256e6+1102516384e8+4365965536e10+6936931048576e12
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+43?e2+6445?e4+256175?e6+1638411025?e8+6553643659?e10+1048576693693?e12 B′=38e2+1532e4+5251024e6+22054096e8+72765131072e10+297297524288e12
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′=83?e2+3215?e4+1024525?e6+40962205?e8+13107272765?e10+524288297297?e12 C′=15256e4+1051024e6+220516384e8+1039665536e10+14864858388608e12
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′=25615?e4+1024105?e6+163842205?e8+6553610396?e10+83886081486485?e12 D′=353072e6+1054096e8+10395262144e10+550551048576e12
D'=\frac{35}{3072}e^6+\frac{105}{4096}e^8+\frac{10395}{262144}e^{10}+\frac{55055}{1048576}e^{12}D′=307235?e6+4096105?e8+26214410395?e10+104857655055?e12 E′=315131072e8+3465524288e10+990998388608e12
E'=\frac{315}{131072}e^8+\frac{3465}{524288}e^{10}+\frac{99099}{8388608}e^{12}E′=131072315?e8+5242883465?e10+838860899099?e12 F′=6931310720e10+99095242880e12
F'=\frac{693}{1310720}e^{10}+\frac{9909}{5242880}e^{12}F′=1310720693?e10+52428809909?e12 G′=10018388608e12
G'=\frac{1001}{8388608}e^{12}G′=83886081001?e12其中，eee为地球椭球第一偏心率，aaa为地球椭球长半轴。
注意1：当纬度已经为弧度(π=180°)时ρ=1;\color{blue}注意1：当纬度已经为弧度(\pi=180°)时\rho=1;注意1：当纬度已经为弧度(π=180°)时ρ=1;
注意2：中国地区内，为了避免出现负数，y需要加上500000m.\color{blue}注意2：中国地区内，为了避免出现负数，y需要加上500000m.注意2：中国地区内，为了避免出现负数，y需要加上500000m.
在WGS-84坐标系中，a=6378137m，f=1/298.257223563，e2=2f?f2，e′2=(2f?f2)/(1?f2)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 [x,y] = GaussProWGS84(Lat,Lon)
% Lat: Latitude(rad)
% Lon: longitude(rad)
% REF//程鹏飞,成英燕,文汉江,等.2000国家大地坐标系实用宝典[M].
% //北京:测绘出版社,:144-148.
Lat = Lat*pi/180;
Lon = Lon*pi/180;
MedLon = 114*pi/180; %武汉的中央子午线经度
Eth.R0 = 6378137.0;
Eth.f = 1/298.257223563;
Eth.e12 = 2*Eth.f - Eth.f*Eth.f; % 0.00669437999014132
Eth.e22 = Eth.e12/((1 - Eth.f)*(1 - Eth.f));
%% 高斯投影正算公式
RN = Eth.R0/sqrt(1 - Eth.e12*sin(Lat)*sin(Lat));
Lon = Lon - MedLon;
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;
NTBLP = RN*tnLat*cs2Lat*Lon2;
coe1 = (5 - tn2Lat + 9*Eta2 + 4*Eta2*Eta2)*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;
end
function X0 = Merdian(Eth,Lat)
% REF//过家春.子午线弧长公式的简化及其泰勒级数解释[J].测绘学报,,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]. 北京:测绘出版社,.
[2] 孔祥元, 郭际明. 大地测量学基础[M]. 武汉: 武汉大学出版社, .
[3] 过家春.子午线弧长公式的简化及其泰勒级数解释[J].测绘学报,,43(2):125-130.
完毕
MATLAB实现高斯-克吕格投影正算相关教程
