高斯投影正算与反算的理论方法与实

一、实习背景高斯投影是一种广泛应用的地图投影方法，它将地球表面的经纬度坐标转换为平面直角坐标。

高斯投影在测绘、地理信息系统、地图编制等领域有着重要的应用。

为了更好地掌握高斯投影的相关知识，提高自己的实践能力，我们进行了高斯投影反算的实习。

二、实习目的1. 理解高斯投影的基本原理和方法；2. 掌握高斯投影反算的计算步骤；3. 提高自己的实践操作能力；4. 培养团队协作精神。

三、实习内容1. 高斯投影原理高斯投影是一种等角投影，其基本原理是将地球椭球面上的经纬度坐标转换为平面直角坐标。

高斯投影具有以下特点：（1）等角投影：保持地球椭球面上任意两点间的夹角不变；（2）等积投影：保持地球椭球面上任意两块区域的面积比不变；（3）高斯-克吕格投影：以中央子午线和赤道为基准线，将地球椭球面投影到平面上。

2. 高斯投影反算步骤高斯投影反算是指将平面直角坐标转换为地球椭球面上的经纬度坐标。

其计算步骤如下：（1）计算投影面大地坐标（φ，λ）：根据给定的平面直角坐标（X，Y），利用高斯投影公式计算投影面大地坐标（φ，λ）；（2）计算大地坐标（φ，λ）：根据投影面大地坐标（φ，λ）和投影带参数，计算大地坐标（φ，λ）；（3）计算经纬度坐标（B，L）：根据大地坐标（φ，λ）和椭球参数，计算经纬度坐标（B，L）。

3. 实习过程在实习过程中，我们首先学习了高斯投影的基本原理和方法，了解了高斯投影在地图编制、地理信息系统等领域的应用。

然后，我们通过查阅相关资料，掌握了高斯投影反算的计算步骤。

在实践操作环节，我们使用高斯投影软件，对给定的平面直角坐标进行反算，得到对应的经纬度坐标。

在操作过程中，我们遇到了一些问题，如坐标转换误差、投影带参数设置等。

通过查阅资料、请教老师，我们解决了这些问题，最终完成了实习任务。

四、实习总结通过本次高斯投影反算实习，我们取得了以下成果：1. 掌握了高斯投影的基本原理和方法；2. 熟悉了高斯投影反算的计算步骤；3. 提高了实践操作能力；4. 培养了团队协作精神。

正形投影的一般条件基本出发点：在正形投影中，长度比与方向无关。

1、长度比的通用公式如图4-42，在微分直角三角形P1P2P3及P1′P2′P3′中有：其中l=L-L0，L0通常是中央子午线的经度，L是P点的经度令：()()222222d=d cos dd=d dS M B N B ls x y++(1)m平方可为：()()()22222222222d d d d dd d cos d dcos dcoss x y x ymS M B N B l M BN B lN B++⎛⎫===⎪⎡⎤⎝⎭+⎛⎫+⎢⎥⎪⎝⎭⎢⎥⎣⎦(2)为简化公式，令：ddcosM BqN B=dc o sB M BqN B=⎰(3) q称为等量纬度，因为它只与纬度B有关。

这样，式(2)可表示为：()()222222d dd dx ymr q l+=⎡⎤+⎣⎦(4)我们投影的目的是：建立平面坐标xy和大地坐标BL之间的函数关系，由式(3)可知，即建立xy和bl的函数关系。

令()(),,x x l q y y l q==(5) 对上式进行全微分可得：d d dd d dx xx q lq ly yy q lq l∂∂⎧=+⎪∂∂⎪⎨∂∂⎪=+⎪∂∂⎩(6)将上式代入式(1)中第二项，并令：2222x yEq qx x y yFq l q lx yGl l⎧⎛⎫⎛⎫∂∂=+⎪ ⎪ ⎪∂∂⎪⎝⎭⎝⎭⎪∂∂∂∂⎪=+⎨∂∂∂∂⎪⎪∂∂⎛⎫⎛⎫⎪=+⎪ ⎪∂∂⎪⎝⎭⎝⎭⎩(7)可得： ()()()()222d d 2d d d s E q F q l G l =++ (8) 则式(4)可写为： ()()()()()()222222d 2d d d d d E q F q l G l m r q l ++=⎡⎤+⎣⎦ (9)2 柯西-黎曼条件在上式引入方向，如图4-42所示：2313d d cot d d P P M B q A PP r l l === (10) 即： d tan d l A q = (11)将式(11)代入式(9)可得：注意sec 1cos A A =()()()()()222222222222222d 2tan d tan d d tan d 2tan tan sec cos 2sin cos sin E q F A q G A q m r q A q E F A G A r AE AF A AG A r ++=⎡⎤+⎣⎦++=++=(12)要想让m 和A 无关，必须使F=0，E=G ，即22220x x y y q l q l x y x y q q l l ∂∂∂∂⎧+=⎪∂∂∂∂⎪⎨⎛⎫⎛⎫∂∂∂∂⎛⎫⎛⎫⎪+=+⎪ ⎪ ⎪ ⎪⎪∂∂∂∂⎝⎭⎝⎭⎝⎭⎝⎭⎩ (13) 由上式第一式可得：y y x q l x lq∂∂∂∂∂=-∂∂∂(14)代入第二式可得： 222222y x y y x lq q q q x q ∂⎛⎫⎪⎡⎤⎛⎫⎛⎫⎛⎫⎛⎫∂∂∂∂∂⎝⎭+=+⎢⎥ ⎪ ⎪ ⎪ ⎪∂∂∂∂⎢⎥⎝⎭⎝⎭⎝⎭⎝⎭⎛⎫∂⎣⎦⎪∂⎝⎭(15) 消去公共项可得： 22x y q l ⎛⎫∂∂⎛⎫= ⎪ ⎪∂∂⎝⎭⎝⎭ (16)开方并代入式(13)的第一项：x y q l x y lq ∂∂⎧=⎪∂∂⎪⎨∂∂⎪=-⎪∂∂⎩ (17)高斯投影坐标正算高斯投影三条件：L0为直线；L0长度不变；正形投影 1、幂级数展开公式（x 偶y 奇）l /ρ微小量(ρ''=206265)，可进行级数展开，可得：2402435135x m m l m l y m l m l m l ⎧=+++⎪⎨=+++⎪⎩(18) 式中mi 为待定系数，是q 、B 的函数。

高斯投影正、反算及换带程序执行条件※数组投影选择T、换算点个数“Z=0 F≠0”、=0正算0、≠0反算※坐标系选择“54 ≠54”、=54换算为1954年北京坐标系输入54、≠54换算为1988年西安坐标系M、中央子午线经度（°′″）输入※大地坐标I、序列号B、L：大地纬度和经度（地理坐标）（°′″）※高斯平面坐标轴子午线I、序列号X、Y：高斯平面坐标（m） Z、轴子午线（°）输出※大地坐标子午收敛角N、序列号B、L：大地纬度和经度（地理坐标）（°′″） R、子午收敛角（°′″）※高斯平面坐标子午收敛角N、序列号X、Y：高斯平面坐标（m） R、子午收敛角（°′″）注：1、程序执行前必须进行数组定位。

如：Defm 10 T×2=5×2=102、Y坐标值要去掉带号及避免出现负值的500公里；4、本程序运算时，各已知数据、观测变量不会随之变化，可非常方便地进行各数据的核对;5、本程序在进行换带计算时采用的是间接换带计算法。

Prog GSXYDefm 10:TA“Z=0 F≠0”G“54 ≠54”Z:Fixm:I=0:「b」0:I=I+1◢J=2I-1:M=Z[J:L=Z[J+1:A=0=>Prog“3”:B=M:M=L+Z:Prog“3”:L=M:{BL}:M=B:Prog“2”: B=M:M=L:Prog“2”:L=M-Z:≠>X=M:Y=L:{XY}:B=X:L=Y⊿Z[J]=B:Z[J+1]=L:I<T=>Goto 0⊿G=54=>C=6399698.90178271:E=.006738525414684:≠>C=6399596.65198801:E=.006 739501819473⊿I=0:「b」0:I“N”=I+1◢J=2I-1:B=Z[J:L=Z[J+1:A≠0=>X=B:Y=L:Goto 2⊿S=sin B:G=54=>F=111134.8611B-(32 005.7799S+133.9238S∧3+.6973S∧5+.0039S∧7)cos B:≠>F=111133.0047B-(32009.857 S+133.9602S∧3+.6976S∧5+.0039S∧7)cos B⊿U=√Ecos B:V=√(1+U2:N=C÷V:W=tan B: M=cos B(Lπ÷180:X=F+NW(.5M2+1┛24(5-W2+9U2+4U∧4)M∧4+1┛720(61-58W2+W∧4)M∧6◢Y=N(M+1┛6(1-W 2+U 2)M ∧3+1┛120(5-18W 2+W ∧4+14U 2-58U 2W 2)M ∧5◢M=W ┛π(180M+60(1+3U 2+2U ∧4)M ∧3+12(2-W 2)M ∧5:Goto 3:「b 」2:W=E ﹣6X-3:G=54=>F=27.11115372595+9.024********W-.00579740442W 2-4.3532572E ﹣4W ∧3+4.857285E ﹣5W ∧4+2.15727E ﹣6W ∧5-1.9399E ﹣7W ∧6:≠>F=27.11162289465+9.024********W-.00579850656W2-4.3540029E ﹣4W ∧3+4.858357E ﹣5W ∧4+2.15769E ﹣6W ∧5-1.9404E ﹣7W ∧6⊿U=√Ecos F:V=√(1+U 2:Q=YV ÷C:W=tan F:M=F-(1+U 2)W ┛π(90Q 2-7.5(5+3W 2+U 2-9U 2W 2)Q ∧4+.25(61+90W 2+45W ∧4)Q ∧6:Prog “3”:B=M ◢M=Z+1┛(πcos F)(180Q-30(1+2W 2+U 2)Q ∧3+1.5(5+28W 2+24W ∧4)Q ∧5:Prog “3”:L=M ◢M=W ┛π(180Q-60(1+W 2-U 2)Q ∧3+12(2+5W 2+3W ∧4)Q ∧5:「b 」3:Prog “3”:R=M ◢ I<T=>Goto 1⊿“END ”概要说明：我国的经度范围西边自73°起，东边至135°，可分成6°带共11带或3°共22带。

高斯投影正反算原理高斯投影是一种常用于地图制图的投影方式，也被广泛应用于其他领域的空间数据处理。

高斯投影正反算是对于已知的地球坐标系上的位置（经纬度），通过计算得到该点的平面坐标（东、北坐标），或者对于已知的平面坐标（东、北坐标），通过计算得到该点的地球坐标系上的位置（经纬度）的过程。

本文将详细介绍高斯投影正反算的原理。

一、高斯投影简介高斯投影是一种圆锥投影，其投影面在地球表面的某个经线上，也就是说，投影面是以该经线为轴的圆锥面。

经过对圆锥体的调整后，使其切于地球椭球面，在该经线上进行投影，同时保持沿该经线方向的比例尺一致，从而达到地图上各点在包括该经线的垂直面上映射的目的。

这种投影方式在某一特定区域内得到高精度的结果，因此广泛应用于地图制图。

二、高斯投影数学模型对于高斯投影正反算，需要先建立高斯投影坐标系与地球坐标系的转换模型。

1.高斯投影坐标系的建立高斯投影坐标系的建立需要确定圆锥面的基本参数，首先需要确定其所处的中央子午线，再确定该子午线上的经度为零点，并利用该经线上某一点的经度和该点的高度来确定该点所在的圆锥体。

圆锥体的底面包括所有与地球椭球面相切的圆面，通过对这些圆面进行调整，使得圆锥体转动后能够在中央子午线上进行投影。

在此基础上，可建立高斯投影坐标系，其中投影面为圆锥面，且中央子午线与投影面的交点称为该投影坐标系的中心，投影面的上端点和下端点分别对应正北方向和正南方向。

2.地球坐标系的建立地球坐标系是以地球椭球体为基础建立的，其坐标系原点确定为地球椭球体上的一个特定点。

在已知该点经纬度和高度的前提下，可确定以该点为中心的地球椭球体，并可根据它与地球坐标系之间的转换关系得到平面坐标系。

3.高斯投影坐标系与地球坐标系之间的转换关系由于高斯投影坐标系与地球坐标系存在不同的坐标体系和基准面，因此需要通过数学关系式来建立它们之间的转换关系。

（1）高斯投影坐标系转地球坐标系：已知高斯投影坐标系中任意一点的东北坐标（N，E），以及所属的中央子午线经度λ0、椭球参数a和e，则可通过以下公式求出该点的地球坐标系经纬度（φ，λ）和高度H：A0为以地球椭球体中心为原点，高斯投影坐标系中心投影坐标为（0，0）的点到椭球面的距离。

高斯投影坐标正反算一、基本思想：高斯投影正算公式就是由大地坐标（L ，B ）求解高斯平面坐标（x ，y ），而高斯投影反算公式则是由高斯平面坐标（x ，y ）求解大地坐标（L ，B ）。

二、计算模型：基本椭球参数：椭球长半轴a椭球扁率f椭球短半轴：(1)b a f =-椭球第一偏心率：e a= 椭球第二偏心率：e b'=高斯投影正算公式：此公式换算的精度为0.001m6425644223422)5861(cos sin 720)495(cos 24cos sin 2l t t B B N l t B simB N l B B N X x ''+-''+''++-''+''⋅''+=ρηηρρ 5222425532233)5814185(cos 120)1(cos 6cos l t t t B N l t B N l B N y ''-++-''+''+-''+''⋅''=ηηρηρρ其中：角度都为弧度B 为点的纬度，0l L L ''=-，L 为点的经度，0L 为中央子午线经度； N 为子午圈曲率半径，1222(1sin )N a e B -=-；tan t B =； 222cos e B η'=1803600ρπ''=*其中X 为子午线弧长：2402464661616sin cos ()(2)sin sin 33X a B B B a a a a a B a B ⎡⎤=--++-+⎢⎥⎣⎦02468,,,,a a a a a 为基本常量，按如下公式计算：200468242684468686883535281612815722321637816323216128m a m m m m m m a m m m a m m m m a m a ⎧=++++⎪⎪⎪=+++⎪⎪⎪=++⎨⎪⎪=+⎪⎪⎪=⎪⎩02468,,,,m m m m m 为基本常量，按如下公式计算：22222020426486379(1);;5;;268m a e m e m m e m m e m m e m =-====；高斯投影反算公式：此公式换算的精度为0.0001’’.()()()()2222243246532235242225053922461904572012cos 6cos 5282468120cos f f f f f f f f f f f f f f f f f f f f f ff f f f f f ft t B B y t t yM N M N t y t t yM N y y l t N B N B y t t t N B L l L ηηηηη=-+++--++=-+++++++=+其中： 0L 为中央子午线经度。

高斯投影坐标正反算及程序设计目录1研究背景和意义 (1)2正形投影的特性及其公式 (2)2.1正形投影特性 (2)2.2正形投影的一般条件 (2)3高斯坐标正反算公式 (6)3.1高斯投影正算 (6)3.2高斯坐标反算 (7)4使用matlab编写高斯坐标正反算函数 (10)4.1高斯坐标正算函数编写 (10)4.2高斯坐标反算函数编写 (13)5高斯坐标正反算程序的界面设计 (17)5.1单点换算模块的编写 (17)5.2批量换算模块的编写 (23)总结与讨论 (33)致谢 (34)参考文献 (34)附录A (35)1研究背景和意义高斯投影即高斯-克吕格投影这个投影是在1920年左右被拟定的，拟定者是德国人高斯，一位著名的数学家、物理学家、天文学家[1]，之后有又有一名叫克吕格的测量学家在1912 年对投影公式加以补充，所以这个投影的全称是高斯-克吕格投影，又名"等角横切椭圆柱投影”，是地球椭球面和平面间正形投影的一种。

即高斯投影是正形投影的一种，所以它具有正形投影的特点，即角度不变性、图形相似性以及在某点方向上长度比的同一性[2]。

而为了使投影区域的长度变形不致过大，采用的措施就是分带。

这样既保证了长度变形不致过大，又可以在不同的投影带里用一样的数学公式来进行各种大地问题的计算，且带与带间的互相换算也能用相同的公式和方法进行[3]。

正是因为这种特点高斯投影被广泛应用到了测绘工作中。

虽然高斯投影作用非常广泛但是因为其复杂繁琐的计算过程往往需要通过程序设计来达到简化计算，提高效率目的。

本文所使用的MATLAB 软件编程技术广泛应用于测量数据处理领域，MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室）。

是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境[4]。

它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，能够开发出界面友好、使用方便的图形界面[5]，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。

高斯投影正算与反算的理论方法与实现代码高斯投影是正形投影的一种，同一坐标系中的高斯投影换带计算公式是根据正形投影原理推导出的两个高斯坐标系间的显函数式。

在同一大地坐标系中(例如1954北京坐标系或1980西安坐标系)，如果两个高斯坐标系只是主子午线的经度不同，那么显函数式前的系数可以根据坐标系使用的椭球元素和主子午线经度唯一确定。

但如果两个高斯坐标系除了主子午线的经度不同以外，还存在其他线性系，则将线性变换公式代入换带计算的显函数式中，仍然可以得到严密的坐标变换公式。

此时显函数式前的系数等价于使用两个坐标系主子午线的经度和线性变换参数联合求解得到的，可以唯一确定。

//6度带宽 54北京坐标系//高斯投影由大地坐标(Unit:Metres)反算经纬度(Unit:DD)void GaussProjInvCal(double X, double Y, double *longitude, double *latitude) {int ProjNo; int ZoneWide; ////带宽double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;double e1,e2,f,a, ee, NN, T,C, M, D,R,u,fai, iPI;iPI = 0.0174532925199433; ////3.1415926535898/180.0;a = 6378245.0; f = 1.0/298.3; //54年北京坐标系参数////a=6378140.0; f=1/298.257; //80年西安坐标系参数ZoneWide = 6; ////6度带宽ProjNo = (int)(X/1000000L) ; //查找带号longitude0 = (ProjNo-1) * ZoneWide + ZoneWide / 2;longitude0 = longitude0 * iPI ; //中央经线X0 = ProjNo*1000000L+500000L;Y0 = 0;xval = X-X0; yval = Y-Y0; //带内大地坐标e2 = 2*f-f*f;e1 = (1.0-sqrt(1-e2))/(1.0+sqrt(1-e2));ee = e2/(1-e2);M = yval;u = M/(a*(1-e2/4-3*e2*e2/64-5*e2*e2*e2/256));fai = u+(3*e1/2-27*e1*e1*e1/32)*sin(2*u)+(21*e1*e1/16-55*e1*e1*e1*e1/32)*sin( 4*u)+(151*e1*e1*e1/96)*sin(6*u)+(1097*e1*e1*e1*e1/512)*sin(8*u);C = ee*cos(fai)*cos(fai);T = tan(fai)*tan(fai);NN = a/sqrt(1.0-e2*sin(fai)*sin(fai));R = a*(1-e2)/sqrt((1-e2*sin(fai)*sin(fai))*(1-e2*sin(fai)*sin(fai))*(1-e2*sin (fai)*sin(fai)));D = xval/NN;//计算经度(Longitude) 纬度(Latitude)longitude1 = longitude0+(D-(1+2*T+C)*D*D*D/6+(5-2*C+28*T-3*C*C+8*ee+24*T*T)*D *D*D*D*D/120)/cos(fai);latitude1 = fai -(NN*tan(fai)/R)*(D*D/2-(5+3*T+10*C-4*C*C-9*ee)*D*D*D*D/24+(61+90*T+298*C+45*T*T-256*ee-3*C*C)*D*D*D*D*D*D/720);//转换为度 DD*longitude = longitude1 / iPI;*latitude = latitude1 / iPI;}//高斯投影由经纬度(Unit:DD)正算平面坐标(含带号，Unit:Metres)void GaussProjCal(double longitude, double latitude, double *X, double *Y) {int ProjNo=0; int ZoneWide; ////带宽double longitude1,latitude1, longitude0,latitude0, X0,Y0, xval,yval;double a,f, e2,ee, NN, T,C,A, M, iPI;iPI = 0.0174532925199433; ////3.1415926535898/180.0;ZoneWide = 6; ////6度带宽a=6378245.0; f=1.0/298.3; //54年北京坐标系参数////a=6378140.0; f=1/298.257; //80年西安坐标系参数ProjNo = (int)(longitude / ZoneWide) ;longitude0 = ProjNo * ZoneWide + ZoneWide / 2;longitude0 = longitude0 * iPI ;latitude0=0;longitude1 = longitude * iPI ; //经度转换为弧度latitude1 = latitude * iPI ; //纬度转换为弧度e2=2*f-f*f;ee=e2*(1.0-e2);NN=a/sqrt(1.0-e2*sin(latitude1)*sin(latitude1));T=tan(latitude1)*tan(latitude1);C=ee*cos(latitude1)*cos(latitude1);A=(longitude1-longitude0)*cos(latitude1);M=a*((1-e2/4-3*e2*e2/64-5*e2*e2*e2/256)*latitude1-(3*e2/8+3*e2*e2/32+45*e2*e2 *e2/1024)*sin(2*latitude1)+(15*e2*e2/256+45*e2*e2*e2/1024)*sin(4*latitude1)-(35 *e2*e2*e2/3072)*sin(6*latitude1));xval = NN*(A+(1-T+C)*A*A*A/6+(5-18*T+T*T+72*C-58*ee)*A*A*A*A*A/120);yval = M+NN*tan(latitude1)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24+(61-58*T+T*T+600*C-330*ee)*A*A*A*A*A*A/720);X0 = 1000000L*(ProjNo+1)+500000L;Y0 = 0;xval = xval+X0; yval = yval+Y0;*X = xval;*Y = yval;}NN卯酉圈曲率半径，测量学里面用N表示M为子午线弧长，测量学里用大X表示fai为底点纬度，由子午弧长反算公式得到，测量学里用Bf表示R为底点所对的曲率半径，测量学里用Nf表示。

高斯投影正反算学院：资源与环境工程工程学院专业：测绘工程 学号：X51414012 姓名：孙超一、高斯投影概述想象有一个椭圆柱面横套在地球椭球体外面，并与某一条子午线相切，椭圆柱的中心轴通过椭球体的中心，然后用一定投影方法，将中央子午线两侧各一定经差范围内的地区投影到椭圆柱面上，再将此柱面展开即成为投影面。

高斯投影由于是正形投影，故保证了投影的角度不变性，图形的相似性以及在某点各方向上长度比的同一性。

由于采用了同样法则的分带投影，这即限制了长度变形，又保证了在不同投影带中采用相同的简便公式和数表进行变形引起的各项改正的计算，并且带与带间的互相换算也能用相同的公式和方法进行。

高斯投影的这些优点必将使它得到广泛的推广和具有国际意义。

二、高斯投影坐标正算公式1.高斯投影必须满足以下三个条件 1）中央子午线投影后为直线 2）中央子午线投影后长度不变 3）投影具有正形性质，即正形投影条件2.高斯正算公式推导1）由第一个条件可知，由于地球椭球体是一个旋转椭球体，所以高斯投影必然有这样一个性质，即中央子午线东西两侧的投影必然对称于中央子午线。

2）由于高斯投影是换带投影，在每带内经差ｌ是不大的，ｌρ是一个微小量，所以可以将X=X （ｌ，q ），Y=Y （l ，q ）展开为经差为l 的幂级数，它可写成如下的形式X=m 0+m 2l 2+m 4l 4+…Y=m 1l+m 3l 2+m 5l 5+…式中m 0,m1,m2,…是待定系数，他们都是纬度B 的函数。

3）由第三个条件：∂y ∂l =∂x ∂q 和∂x ∂l =-∂y∂q ，将上式分别对l 和q 求偏导2340123423401234...........x m m l m l m l m l y n n l n l n l n l =+++++=+++++可得到下式0312123403121234111,,,, 234111,,,,234dm dm dm dm n n n n dq dq dq dq dn dn dn dn m m m m dq dq dq dq ⎧====⎪⎪⎨⎪=-=-=-=-⎪⎩经过计算可以得出232244524632235242225sin cos sin cos (594)224 sin cos (6158)720cos cos (1)6cos (5181458)120N N x X B B l B B t l NB B t t l Ny N B l B t l NB t t t l ηηηηη=+⋅+-+++-+=⋅+-++-++-三、高斯投影坐标反算公式推导1.思路：级数展开，应用高斯投影三个条件，待定系数法求解。

实验报告高斯正反算*名：**学号：********班级：测绘10-1班指导老师：***目录一、实验目的-----------------2二、实验内容及步骤--------2三、程序代码-----------------4四、流程图--------------------24五、运算结果-----------------26六、实验感想-----------------29一、实验目的1、了解高斯正反算的基本思想。

2.学会编写高斯正反算程序，加深理解。

二、实验内容及步骤高斯投影正算公式是由大地坐标（L，B）求解高斯平面坐标（x，y），而高斯投影反算公式则是由高斯平面坐标（x，y）求解大地坐标（L，B）。

现行的高斯投影用表都是采用克拉索夫斯基椭球参数，这次编程计算不仅采用这种椭球参数，还可以选择IAG椭球进行计算。

编程环境是在VC下，采用C++语言编写。

程序主要分为两部分，第一部分是高斯正反算函数，第二部分是主函数。

三、程序代码1.高斯投影正算// mydlg1.cpp : implementation file#include "stdafx.h"#include "高斯正反算.h"#include "mydlg1.h"#include"math.h"#define P 206264.8062471#ifdef _DEBUG#define new DEBUG_NEW#undef THIS_FILEstatic char THIS_FILE[] = __FILE__;#endifCmydlg1::Cmydlg1(CWnd* pParent /*=NULL*/) : CDialog(Cmydlg1::IDD, pParent){//{{AFX_DATA_INIT(Cmydlg1)m_num1 = 0;m_num2 = 0;m_num3 = 0.0;m_num4 = 0;m_num5 = 0;m_num6 = 0.0;m_num7 = 0;m_num8 = 0.0;m_num9 = 0.0;//}}AFX_DATA_INIT}void Cmydlg1::DoDataExchange(CDataExchange* pDX) {CDialog::DoDataExchange(pDX);//{{AFX_DATA_MAP(Cmydlg1)DDX_Text(pDX, IDC_EDIT1, m_num1);DDX_Text(pDX, IDC_EDIT2, m_num2);DDX_Text(pDX, IDC_EDIT3, m_num3);DDX_Text(pDX, IDC_EDIT4, m_num4);DDX_Text(pDX, IDC_EDIT5, m_num5);DDX_Text(pDX, IDC_EDIT6, m_num6);DDX_Text(pDX, IDC_EDIT7, m_num7);DDX_Text(pDX, IDC_EDIT8, m_num8);DDX_Text(pDX, IDC_EDIT9, m_num9);//}}AFX_DATA_MAP}BEGIN_MESSAGE_MAP(Cmydlg1, CDialog) //{{AFX_MSG_MAP(Cmydlg1)ON_BN_CLICKED(IDC_RADIO1, OnRadio1)ON_BN_CLICKED(IDC_RADIO2, OnRadio2)ON_BN_CLICKED(IDC_BUTTON2, OnButton2)ON_BN_CLICKED(IDC_BUTTON1, OnButton1)ON_BN_CLICKED(IDC_BUTTON3, OnButton3)//}}AFX_MSG_MAPEND_MESSAGE_MAP()/////////////////////////////////////////////////////////// //////////////////// Cmydlg1 message handlersvoid Cmydlg1::OnRadio1(){// TODO: Add your control notification handler code here}void Cmydlg1::OnRadio2(){// TODO: Add your control notification handler code here}void Cmydlg1::OnButton2(){// TODO: Add your control notification handler code here CEdit* pedt1 = (CEdit*)GetDlgItem(IDC_EDIT1);pedt1->SetWindowText("");CEdit* pedt2 = (CEdit*)GetDlgItem(IDC_EDIT2);pedt2->SetWindowText("");CEdit* pedt3 = (CEdit*)GetDlgItem(IDC_EDIT3);pedt3->SetWindowText("");CEdit* pedt4 = (CEdit*)GetDlgItem(IDC_EDIT4); pedt4->SetWindowText("");CEdit* pedt5= (CEdit*)GetDlgItem(IDC_EDIT5); pedt5->SetWindowText("");CEdit* pedt6 = (CEdit*)GetDlgItem(IDC_EDIT6); pedt6->SetWindowText("");CEdit* pedt7 = (CEdit*)GetDlgItem(IDC_EDIT7); pedt7->SetWindowText("");CEdit* pedt8 = (CEdit*)GetDlgItem(IDC_EDIT8); pedt8->SetWindowText("");CEdit* pedt9 = (CEdit*)GetDlgItem(IDC_EDIT9); pedt9->SetWindowText("");}void Cmydlg1::OnCancel(){// TODO: Add extra cleanup hereCDialog::OnCancel();}void Cmydlg1::OnButton1(){// TODO: Add your control notification handler code here double N,a0,a3,a4,a5,a6,B,L,l,x,y;UpdateData();B=m_num1*3600+m_num2*60+m_num3;L=m_num4*3600+m_num5*60+m_num6;B=B/P;l=(L-m_num7*3600)/P;N=6399698.902-(21562.267-(108.973-0.612*cos(B)*cos(B))*cos( B)*cos(B))*cos(B)*cos(B);a0=32140.404-(135.3302-(0.7092-0.0040*cos(B)*cos(B))*cos(B) *cos(B))*cos(B)*cos(B);a4=(0.25+0.00252*cos(B)*cos(B))*cos(B)*cos(B)-0.04166;a6=(0.166*cos(B)*cos(B)-0.084)*cos(B)*cos(B);a3=(0.3333333+0.001123*cos(B)*cos(B))*cos(B)*cos(B)-0.16666 67;a5=0.0083-(0.1667-(0.1968+0.0040*cos(B)*cos(B))*cos(B)*cos( B))*cos(B)*cos(B);x=6367558.4969*B-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*sin(B)*co s(B);y=(1+(a3+a5*l*l)*l*l)*l*N*cos(B);m_num8=x;m_num9=y;UpdateData(FALSE);}void Cmydlg1::OnButton3(){// TODO: Add your control notification handler code here double N,a0,a3,a4,a5,a6,B,L,l,x,y;UpdateData();B=m_num1*3600+m_num2*60+m_num3;L=m_num4*3600+m_num5*60+m_num6;B=B/P;l=(L-m_num7*3600)/P;N=6399596.652-(21565.045-(108.996-0.603*cos(B)*cos(B))*cos( B)*cos(B))*cos(B)*cos(B);a0=32144.5189-(135.3646-(0.7034-0.0041*cos(B)*cos(B))*cos(B )*cos(B))*cos(B)*cos(B);a4=(0.25+0.00253*cos(B)*cos(B))*cos(B)*cos(B)-0.04167;a6=(0.167*cos(B)*cos(B)-0.083)*cos(B)*cos(B);a3=(0.3333333+0.001123*cos(B)*cos(B))*cos(B)*cos(B)-0.16666 67;a5=0.00878-(0.1702-(0.20382*cos(B)*cos(B)))*cos(B)*cos(B);x=6367452.1328*B-(a0-(0.5+(a4+a6*l*l)*l*l)*l*l*N)*sin(B)*co s(B);y=(1+(a3+a5*l*l)*l*l)*l*N*cos(B);m_num8=x;m_num9=y;UpdateData(FALSE);}2.高斯投影反算// mydlg2.cpp : implementation file//#include "stdafx.h"#include "高斯正反算.h"#include "mydlg2.h"#include "math.h"#define P 206264.8062471#define PI 3.1415926535898#ifdef _DEBUG#define new DEBUG_NEW#undef THIS_FILEstatic char THIS_FILE[] = __FILE__;#endif/////////////////////////////////////////////////////////// //////////////////// Cmydlg2 dialogCmydlg2::Cmydlg2(CWnd* pParent /*=NULL*/): CDialog(Cmydlg2::IDD, pParent){//{{AFX_DATA_INIT(Cmydlg2)m_num1 = 0.0;m_num2 = 0.0;m_num4 = 0;m_num5 = 0;m_num6 = 0;m_num7 = 0.0;m_num8 = 0;m_num9 = 0;m_num10 = 0.0;//}}AFX_DATA_INIT}void Cmydlg2::DoDataExchange(CDataExchange* pDX) {CDialog::DoDataExchange(pDX);//{{AFX_DATA_MAP(Cmydlg2)DDX_Text(pDX, IDC_EDIT1, m_num1);DDX_Text(pDX, IDC_EDIT2, m_num2);DDX_Text(pDX, IDC_EDIT4, m_num4);DDX_Text(pDX, IDC_EDIT5, m_num5);DDX_Text(pDX, IDC_EDIT6, m_num6);DDX_Text(pDX, IDC_EDIT7, m_num7);DDX_Text(pDX, IDC_EDIT8, m_num8);DDX_Text(pDX, IDC_EDIT9, m_num9);DDX_Text(pDX, IDC_EDIT10, m_num10);//}}AFX_DATA_MAP}BEGIN_MESSAGE_MAP(Cmydlg2, CDialog) //{{AFX_MSG_MAP(Cmydlg2)ON_BN_CLICKED(IDC_RADIO2, OnRadio2)ON_BN_CLICKED(IDC_RADIO1, OnRadio1)ON_BN_CLICKED(IDC_BUTTON1, OnButton1)ON_BN_CLICKED(IDC_BUTTON2, OnButton2)ON_BN_CLICKED(IDC_BUTTON4, OnButton4)//}}AFX_MSG_MAPEND_MESSAGE_MAP()/////////////////////////////////////////////////////////// //////////////////// Cmydlg2 message handlersvoid Cmydlg2::OnRadio2(){// TODO: Add your control notification handler code here}void Cmydlg2::OnRadio1(){// TODO: Add your control notification handler code here }void Cmydlg2::OnButton1(){// TODO: Add your control notification handler code here CEdit* pedt1 = (CEdit*)GetDlgItem(IDC_EDIT1);pedt1->SetWindowText("");CEdit* pedt2 = (CEdit*)GetDlgItem(IDC_EDIT2);pedt2->SetWindowText("");CEdit* pedt4= (CEdit*)GetDlgItem(IDC_EDIT4);pedt4->SetWindowText("");CEdit* pedt5= (CEdit*)GetDlgItem(IDC_EDIT5);pedt5->SetWindowText("");CEdit* pedt6 = (CEdit*)GetDlgItem(IDC_EDIT6); pedt6->SetWindowText("");CEdit* pedt7 = (CEdit*)GetDlgItem(IDC_EDIT7); pedt7->SetWindowText("");CEdit* pedt8 = (CEdit*)GetDlgItem(IDC_EDIT8); pedt8->SetWindowText("");CEdit* pedt9 = (CEdit*)GetDlgItem(IDC_EDIT9); pedt9->SetWindowText("");CEdit* pedt10 = (CEdit*)GetDlgItem(IDC_EDIT10); pedt10->SetWindowText("");}void Cmydlg2::OnCancel(){// TODO: Add extra cleanup hereCDialog::OnCancel();}void Cmydlg2::OnButton2(){// TODO: Add your control notification handler code here double B,L,l,Bf,Nf,b2,b3,b4,b5,b,Z;UpdateData();b=(m_num1/6367558.4969)*P;b=(b/3600)*(PI/180);Bf=b+(50221746+(293622+(2350+22*cos(b)*cos(b))*cos(b)*cos(b ))*cos(b)*cos(b))*pow(10,-10)*sin(b)*cos(b)*P/3600*(PI/180) ;Nf=6399698.902-(21562.267-(108.973-0.612*cos(Bf)*cos(Bf))*c os(Bf)*cos(Bf))*cos(Bf)*cos(Bf);Z=m_num2/(Nf*cos(Bf));b2=(0.5+0.003369*cos(Bf)*cos(Bf))*sin(Bf)*cos(Bf);b3=0.333333-(0.1666667-0.001123*cos(Bf)*cos(Bf))*cos(Bf)*co s(Bf);b4=0.25+(0.16161+0.00562*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf); b5=0.2-(0.16667-0.0088*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);Bf=Bf*(180/PI);B=Bf-(1-(b4-0.12*Z*Z)*Z*Z)*Z*Z*b2*P/3600;l=(1-(b3-b5*Z*Z)*Z*Z)*Z*P;L=m_num4+l/3600;m_num5=int(B);m_num6=int((B-m_num5)*60);m_num7=((B-m_num5)*60-m_num6)*60;m_num8=int(L);m_num9=int((L-m_num8)*60);m_num10=((L-m_num8)*60-m_num9)*60;UpdateData(FALSE);}void Cmydlg2::OnButton4(){// TODO: Add your control notification handler code here double B,L,l,Bf,Nf,b2,b3,b4,b5,b,Z;UpdateData();b=m_num1/6367558.4969*P;b=(b/3600)*(PI/180);Bf=b+(50221746+(293622+(2350+22*cos(b)*cos(b))*cos(b)*cos(b ))*cos(b)*cos(b))*pow(10,-10)*sin(b)*cos(b)*P/3600*(PI/180) ;Nf=6399698.902-(21562.267-(108.973-0.612*cos(Bf)*cos(Bf))*c os(Bf)*cos(Bf))*cos(Bf)*cos(Bf);Z=m_num2/(Nf*cos(Bf));b2=(0.5+0.00336975*cos(Bf)*cos(Bf))*sin(Bf)*cos(Bf);b3=0.333333-(0.1666667-0.001123*cos(Bf)*cos(Bf))*cos(Bf)*co s(Bf);b4=0.25+(0.161612+0.005617*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf) ;b5=0.2-(0.16667-0.00878*cos(Bf)*cos(Bf))*cos(Bf)*cos(Bf);Bf=Bf*(180/PI);////////B=Bf-(1-(b4-0.147*Z*Z)*Z*Z)*Z*Z*b2*P/3600;l=(1-(b3-b5*Z*Z)*Z*Z)*Z*P;L=m_num4+l/3600;m_num5=int(B);m_num6=int((B-m_num5)*60);m_num7=((B-m_num5)*60-m_num6)*60; m_num8=int(L);m_num9=int((L-m_num8)*60);m_num10=((L-m_num8)*60-m_num9)*60; UpdateData(FALSE);}四、流程图1、高斯投影正算投影2.高斯投影反算五、运算结果1.主界面2.高斯投影正算（克氏椭球）3.高斯投影正算（IAG椭球）4.高斯投影反算（克氏椭球）5. 高斯投影反算（IAG椭球）六、实验感想我在这次实验编程中主要完成了在两种椭球下的高斯正反算。

高斯投影正反算学院：资源与环境工程工程学院专业：测绘工程学号：X51414012姓名：孙超一、高斯投影概述想象有一个椭圆柱面横套在地球椭球体外面，并与某一条子午线相切，椭圆柱的中心轴通过椭球体的中心，然后用一定投影方法，将中央子午线两侧各一定经差范围内的地区投影到椭圆柱面上，再将此柱面展开即成为投影面。

高斯投影由于是正形投影，故保证了投影的角度不变性，图形的相似性以及在某点各方向上长度比的同一性。

由于采用了同样法则的分带投影，这即限制了长度变形，又保证了在不同投影带中采用相同的简便公式和数表进行变形引起的各项改正的计算，并且带与带间的互相换算也能用相同的公式和方法进行。

高斯投影的这些优点必将使它得到广泛的推广和具有国际意义。

二、高斯投影坐标正算公式1.高斯投影必须满足以下三个条件1）中央子午线投影后为直线2）中央子午线投影后长度不变3）投影具有正形性质，即正形投影条件2.高斯正算公式推导1）由第一个条件可知，由于地球椭球体是一个旋转椭球体，所以高斯投影必然有这样一个性质，即中央子午线东西两侧的投影必然对称于中央子午线。

2）由于高斯投影是换带投影，在每带内经差ｌ是不大的，ｌρ是一个微小量，所以可以将 X=X （ｌ，q ），Y=Y （l ，q ）展开为经差为l 的幂级数，它可写成如下的形式X=m 0+m 2l 2+m 4l 4+…Y=m 1l+m 3l 2+m 5l 5+…式中m 0,m1,m2,…是待定系数，他们都是纬度B 的函数。

3）由第三个条件：∂y ∂l =∂x ∂q 和∂x ∂l =-∂y ∂q ，将上式分别对l 和q 求偏导2340123423401234...........x m m l m l m l m l y n n l n l n l n l =+++++=+++++可得到下式0312123403121234111,,,, 234111,,,,234dm dm dm dm n n n n dq dq dq dq dn dn dn dn m m m m dq dq dq dq ⎧====⎪⎪⎨⎪=-=-=-=-⎪⎩L L 经过计算可以得出232244524632235242225sin cos sin cos (594)224sin cos (6158)720cos cos (1) 6cos (5181458)120N N x X B B l B B t l N B B t t l N y N B l B t l N B t t t l ηηηηη=+⋅+-+++-+=⋅+-++-++-三、高斯投影坐标反算公式推导1.思路：级数展开，应用高斯投影三个条件，待定系数法求解。

「高斯投影坐标正反算公式及适合电算的高斯投影公式」高斯投影坐标正反算公式是用于计算高斯投影坐标的数学公式。

高斯投影坐标是一种地理坐标系统，常用于测量和测绘工作中。

高斯投影坐标正算是指已知一个点的经纬度坐标，通过公式计算出该点的高斯投影坐标。

而高斯投影坐标反算是指已知一个点的高斯投影坐标，通过公式计算出该点的经纬度坐标。

一、高斯投影坐标正算公式：已知一个点的经纬度坐标(φ,λ)，其中φ为纬度，λ为经度，以及椭球体参数a、f和中央经线经度L0，可以通过以下步骤计算出该点的高斯投影坐标(X,Y)：1.计算扁率f'：f'=(a-b)/a其中，b=a*(1-f)是椭球体的短半轴。

2.计算黄赤交角ε：ε = atan(b / a)3.计算辅助量t：t = tan(π/4 - φ/2) / [(1 - f' * sin²φ)⁰.⁵ * (1 + e' *sinφ)⁰.⁵]其中，e'=f'*(2-f')是椭球体的第一偏心率。

4.计算辅助量η：η = e'^2 * cos²φ5.计算系数A、B、C和D：A = (L - L0) * cosφC = (L - L0) * cos⁵φ * (5 - tan²φ + 9e'^² + 4e'^⁴ - 24e'^² * tan²φ - 45e'^⁴ * tan²φ)D = (L - L0) * cos⁷φ * (61 - 58tan²φ + tan⁴φ + 270e'^² - 330e'^² * tan²φ)6.计算高斯坐标X和Y：X=k0*a*(A+B/2+C/4+D/6)Y=k0*a*(C/2+D/8)其中，k0是比例系数，一般情况下取1二、高斯投影坐标反算公式：已知一个点的高斯投影坐标(X,Y)，以及椭球体参数a、f、中央经线经度L0、比例系数k0和起始经度L1，可以通过以下步骤计算出该点的经纬度坐标(φ,λ)：1.计算扁率f'：f'=(a-b)/a其中，b=a*(1-f)是椭球体的短半轴。

§8.3高斯投影坐标正反算公式任何一种投影①坐标对应关系是最主要的；②如果是正形投影，除了满足正形投影的条件外（C-R 偏微分方程），还有它本身的特殊条件。

8.3.1高斯投影坐标正算公式： B, x,yl ⇒高斯投影必须满足以下三个条件：①中央子午线投影后为直线；②中央子午线投影后长度不变；③投影具有正形性质，即正形投影条件。

由第一条件知中央子午线东西两侧的投影必然对称于中央子午线，即(8-10)式中，x 为的偶函数，y 为的奇函数；，即，l l 0330'≤l 20/1/≈''''ρl 如展开为的级数，收敛。

l （8-33）+++=++++=553316644220l m l m l m y l m l m l m m x 式中是待定系数，它们都是纬度B 的函数。

,,10m m 由第三个条件知：qyl x l y q x ∂∂-=∂∂∂∂=∂∂,(8-33)式分别对和q 求偏导数并代入上式l (8-34)----=++++++=+++5533156342442204523164253l dqdm l dq dm l dq dm l m l m l m l dqdm l dq dm dq dm l m l m m 上两式两边相等，其必要充分条件是同次幂前的系数应相等，即l(8-35)dq dm m dqdm m dqdm m 2312013121⋅=⋅-==(8-35)是一种递推公式，只要确定了就可依次确定其余各系数。

0m 由第二条件知:位于中央子午线上的点，投影后的纵坐标x 应等于投影前从赤道量至该点的子午线弧长X ，即(8-33)式第一式中，当时有：0=l(8-36)0m X x==顾及(对于中央子午线)B V Mr M B N dq dB M dBdXcos cos 2====得：(8-37,38) B Vc B N r dq dB dB dX dq dX dq dm m cos cos 01===⋅===(8-39)B B Ndq dB dB dm dq dm m cos sin 22121112=⋅-=⋅-=依次求得并代入(8-33)式，得到高斯投影正算公式6543,,,m m m m6425644223422)5861(cos sin 720)495(cos 24cos sin 2lt t B B N lt B simB N l B B N X x ''+-''+''++-''+''⋅''+=ρηηρρ (8-42)5222425532233)5814185(cos 120)1(cos 6cos l t t t B N lt B N l B N y ''-++-''+''+-''+''⋅''=ηηρηρρ8.3.2高斯投影坐标反算公式x,y B,⇒l投影方程：(8-43)),(),(21y x l y x B ϕϕ==满足以下三个条件：①x 坐标轴投影后为中央子午线是投影的对称轴；② x 坐标轴投影后长度不变；③投影具有正形性质，即正形投影条件。