一种建立工程区域水准椭球的方法
【专利摘要】本发明公开了一种建立工程区域水准椭球的方法：在标准椭球方程的基础上构建由待定参数构成的工程区域水准椭球方程，利用测区各个控制点高程异常数据，在测区平均高程面基础上，将各个控制点高程换算到测区平均高程面后所对应的测区平均大地水准面，得到源椭球下测区椭球表面点空间直角坐标；在该坐标下求出平移参数和长短轴参数；通过平移和调整椭球长、短轴，构建出新的工程水准椭球E1，并求出目标椭球下的空间坐标，将转换后的坐标在工程水准椭球E1下进行椭球长、短轴可变的大地坐标转换，进而转换到工程平面坐标，即可用于进一步的精准归算。该方法构建的水准椭球面与工程测区平均高程面对应的大地水准面较好吻合，利于数据的精准归算。
【专利说明】
一种建立工程区域水准椭球的方法
技术领域
[0001] 本发明属于地理测绘科学与技术领域，具体涉及一种建立工程区域水准椭球的方 法。
【背景技术】
[0002] 工程椭球的理想状态是其参考椭球面与测区大地水准面重合，即，椭球面是正常 重力场中的一个等位面。工程参考椭球面是工程测量计算的基准面，又是研究地球形状和 地图投影的参考面，所以工程椭球一定要满足大地参考系统的基本常数(3，13,611，《)，而地 球的重力等位面是一个不规则的曲面，其不满足大地参考系统的基本常数，因此，不存在既 满足大地参考系统基本常数、又满足表面重力等位条件的参考椭球。
[0003] 目前，工程中习惯继承通用椭球模型参数，如1980西安坐标系及CGCS2000坐标系 所采用的椭球分别就是参心与地心椭球，对椭球模型长半轴进行简单调整，建立工程椭球 模型，这种参数的简单调整，远无法达到与工程区域大地水准面最为密合的目的，同时受方 法的限制，有着较大的误差。
【发明内容】
[0004] 本发明提供了一种建立工程区域水准椭球的方法，并构建出与测区在已有数据下 最优密合的椭球参数，解决了现有工程椭球无法与整个工程区域大地水准面较好密合的问 题。
[0005] 本发明采用的技术方案是，一种建立工程区域水准椭球的方法，包括以下步骤：
[0006] 步骤1，确定基于标准椭球的工程水准椭球
[0007] 在标准椭球方程的基础上引入由aE1、bE1、AX E1、AYE1、AZE15个待定参数构成的工 程区域水准椭球方程，
[0009] 式中，XE1、YE1、ZE1为在源椭球下测区椭球表面点空间直角坐标，AX E1、A YE1、AZE1 为平移参数，aEi、bEi分别为长轴和短轴；
[0010] 步骤2，求解空间直角坐标和工程椭球参数
[0011] 利用测区各个控制点高程异常数据，在测区平均高程面基础上，将各个控制点高 程换算到测区平均高程面后所对应的测区平均大地水准面，得到源椭球下测区椭球表面点 空间直角坐标X E1、YE1、ZE1;
[0012] 步骤3，求解工程椭球参数
[0013] 在空间直角坐标乂^士1、2^下，获取测区5组代表性坐标，然后采用最小二乘原理 求解出5个工程椭球参数A XE1、A YE1、A ZE1、aE1、bE1;
[0014] 步骤4，工程水准椭球坐标转换及椭球参数调整
[0015] 将5个工程椭球参数代入公式⑶，通过平移（AXE1、AYE1、AZ E1)、调整椭球长轴aE1 和短轴bE1，构建出新的工程水准椭球El;
[0016]步骤5,椭球空间直角坐标到工程平面坐标的转换
[0017]先将E1椭球的空间直角坐标转换为大地坐标，再将其转换到工程平面坐标，得到 椭球面与测区大地水准面密合的工程区域水准椭球。
[0018] 本发明的特点还在于：
[0019] 上述步骤2具体按照以下步骤实施：
[0020] 步骤2-1，求解大地高
[0021] 如果测区中采用GNSS联测IGS站获取了较高精度的ITRF框架下的Xitrf、Yit rf、Zitrf 坐标，则可以得到对应系统下控制点的相应精度的大地高He。
[0022] 如果测区中没有获取该坐标，可以在GNSS网三维约束平差时采用测站单点定位近 似坐标约束，通过解算得到空间直角坐标后，将各点空间直角坐标以GNSS测量所采用椭球 为目标椭球，利用大地正算公式得到各点位在WGS84或ITRF框架下的大地坐标(B，L，H)，进 而得到大地高He;
[0023] 步骤2-2,解算参考椭球参数的测区各点对应的大地高
[0024] 定义控制点相对源椭球的大地高为He，正常高为Hn，则该点的高程异常|E1 = He-HN;
[0025] 设选定的测区平均高程或者根据工程需求而指定的基准面为H，工程椭球面拟合 所对应的大地高为：馬 < 所1與1…W…尽匕f高程异常为 ... 则工程椭球面拟合所对应的大地高出=H+ CE1；
[0026] 求出He后，根据大地坐标与空间直角坐标转换公式将大地坐标Zf 转换成空间直角坐标[If导到的If 即为空间直角坐标XE1、YE1、 Zei0
[0027] 上述步骤3具体按照以下步骤实施：
[0028] 获取测区5组代表性坐标，由间接平差模型得知：
[0029] = X-X° , I = L- f(Xl') = L-ll\ V :抿-1
[0030] 非线性误差方程为：
[0031] y= f(X)-L
[0032] 式中：V为观测值的改正数向量（残差向量）；X为包含5个未知向量的估值： 又=饫El h ^*^￡1 ^El ；
[0033]将非线性模型(8)式在#处用泰勒级数展开，取至一次项，再由#为测区已知的5个 控制点求得的第一初始值为：
[0034] X°=[a° b° AX° AY0 AZ°]T (11)
[0035] 1为带入X13后式子
的值，有n个;P为单位 权矩阵，得到
[0036] x=0rPBylJfPl (12)
[0037] 通过解算式（12)就可以算出I，由此通过式（^彡算出尤。
[0038] X=：X°+J= a b -hK AI IsZ
[0039] 至此得到了待求的未知参数AXE1、A YE1、AZE1、aE1、bE1。
[0040]上述步骤5中将大地坐标转换到工程平面坐标具体按照以下步骤实施：
[0041] (1)建立大地子午线弧长)(模型：
[0042] X(B)\；\lX = ^(X\iX (16)
[0043] (2)根据大地子午线弧长)(模型建立垂足炜度Bf模型：
[0044] Bf称为垂足炜度或底点炜度，过待算点P做中央子午线的垂线PP' .0P'弧所对应的 弧段即为&，则有
[_]心蕭（丨9)
[0046]对公式(19)变形后得到
[0048]即为求解垂足炜度Bf的模型；
[0049] (3)按照大地正算公式计算工程平面坐标：
[0052] 式中：1(B)为子午线弧长，
为卯酉圈半径，t = tan B，l = i- 为经差，Lo为中央子午线经度，1(B)为从赤道到投影点的椭球面弧长。
[0053]本发明的有益效果是:本发明的建立工程区域水准椭球的方法构建的工程区域水 准椭球面与工程测区平均高程面对应的大地水准面较好吻合，能使测区大地测量中的几何 参考面和物理参考面统一起来，易于实际的计算。与现有其他方法构建的可保证在选定高 程面上，椭球面与测区大地水准面的最好密合，利于数据的精准归算。
【具体实施方式】
[0054] 下面结合【具体实施方式】对本发明作进一步的详细说明。
[0055] 步骤1，确定基于标准椭球的工程水准椭球
[0056] 基于全球性通用椭球，通过改变椭球中心位置及椭球长短轴，使得椭球面与测区 平均高程面对应的大地水准面最优吻合。工程水准椭球的构建必须首先确定测区平均高程 面所对应的大地水准面，该大地水准面可以采用分布于测区的若干具有正常高成果的控制 点近似表达出来。
[0057]标准椭球方程见式(1)。
[0059]式中，a、b为赤道半径，c为极半径，X、Y、Z为标准椭球的空间直角坐标。
[0060]按照空间直角坐标系的定义，旋转椭球体在X和Y方向半径相等，故a = b，得到：
[0062] 理想的工程水准椭球体是椭球面与测区平均大地水准面重合，但受椭球数学表达 的限制，椭球面与测区平均大地水准面重合无法实现，但可以在一定规则下达到椭球面与 大地水准面的最优吻合。因此，在公式(2)的基础上引入平移参数，形成了由长轴、短轴、平 移参数构成的椭球方程，于是式(2)可表达为：
[0064] 式中，Xei、Yei、Zei为在源椭球下测区椭球表面点空间直角坐标，A Xei、A Yei、A Zei 为平移参数，aEi、bEi分别为长轴和短轴。其中，A Xei、A Yei、A ZEi、aEi、bEi为五个待定工程椭 球参数。
[0065] 步骤2，求解空间直角坐标XE1、YE1、ZE1
[0066] XE1、YE1、ZE1可以理解为控制点在所求工程椭球表面投影点的位置在源椭球下的空 间直角坐标。其求解方法为，利用测区各个控制点高程异常数据，在测区平均高程面基础 上，将各个控制点高程换算到测区平均高程面后所对应的测区平均大地水准面，得到源椭 球下测区椭球表面点空间直角坐标Xei、Yei、Zei。具体方法如下：
[0067] 步骤2-1，求解大地坐标
[0068] 如果测区中采用GNSS联测IGS站获取了较高精度的ITRF框架下的Xitrf、Yit rf、Zitrf 坐标，则可以得到对应系统下控制点的相应精度的大地高He。
[0069] 如果测区中没有获取该坐标，可以在GNSS网三维约束平差时采用测站单点定位近 似坐标约束，通过解算得到空间直角坐标后，将各点空间直角坐标以GNSS测量所采用椭球 为目标椭球，利用公式(4)得到点各位在WGS84或ITRF框架下的大地坐标(B，L，H)。 N表不卯酉圈的半径。
[0072]由式(4)得到控制点的经度B、炜度L和大地高He。
[0073]步骤2-2，解算参考椭球参数的测区各点对应大地高的确定方法：
[0074]定义控制点相对源椭球的大地高为He，正常高为Hn，则该点的高程异常|E1求解见 式(5)。
[0075] CE1 = Hb-Hn (5)
[0076] 原则上用于求定参考椭球参数的测区各点的大地高宜采用测区的平均高程或者 选定的某一高程面，但由于各点的高程异常不一致，宜在选定的高程基础上加上对应点位 的高程异常值。
[0077] 设选定的测区平均高程或者根据工程需求而指定的基准面为H，工程椭球 面拟合所对应的大地高为：A …Hy…H;二高程异常为 4El=[4im ^ …C则工程椭球面拟合所对应的大地高由式 (6)得到
[0078] Hg = H+CE1 (6)
[0079] 求出HC后，将大地坐标if //ff按照式（7)转换成空间直角坐标 X° =(N + ) cos Bfl cos Lf'
[0080] Y° ={N+IIa)cosBf\mLm > Z：=[N{\-n + !I^]smB；1 \ (?)
[0081] 式中，N为卯酉圈的半径
[0082] 式(7)得到的Z；^即为空间直角坐标XE1、YE1、ZE1。
[0083] 步骤3、求解工程椭球参数
[0084] 式(3)的五个待定工程椭球参数分别为A Xei、A Yei、A 2￡1、8￡1、匕￡1，要求解该方程， 至少需要均匀分布于测区的5组代表性坐标，然后采用最小二乘原理求解该5个参数。具体 方法如下：
[0085]步骤3-1，由间接平差模型得知：
[0086] J = X-r, t = L - f(X°) = L-ff，V=^-t [0087]非线性误差方程为：
[0088] V = f{X)-L (：g)
[0089 ]式中：V为观测值的改正数向量(残差向量）；龙为包含5个未知向量的估值：
[0090] X = p C
[0091] 由非线性误差方程式(8)知，式中至少有5个方程，且有5个未知数(3个坐标平移参 数以及椭球长、短轴参数）。因此非线性误差方程式(8)是非线性不定方程组，有无穷组解。 在这无穷组解中，必然有一组解能使
[0092] V' PV = (f(X)- L)P(f(X)- L) = min
[0093] 推出：
[0094] x={B'PB) 'B1 PI (9)
[0095] 将非线性模型(8)式在X*3处用泰勒级数展开，取至一次项，得：
[0098] 由X13为测区已知的5个控制点求得的第一初始值为：
[0099] X°=[a° b° AX° AY0 AZ°]T (11)
[0100] 1为带入X13后式子 的值，有n个。P为单位 权矩阵。得到
[0101] X=(BIPB) 'B'P! ( 12)
[0102] 通过解算式（12)就可以算出$。由此通过式（13)算出f。
[0103] X = X° + x - a b AX AF AZ 匕 」 C i ?3 )
[0104] 至此得到了待求的未知参数A XE1、A YE1、A ZE1、aE1、bE1。
[0105]通过这样转换的参数对空间三维坐标进行A X、A Y、A Z平移转换及对椭球参数进 行变换，可以保证椭球面与测区平均大地水准面的较好吻合，可满足常见的工程精密坐标 归算及高程拟合需求。
[0106]步骤4，工程水准椭球坐标转换及椭球参数调整
[0107] 根据前述步骤求解的平移参数么乂￡1、么￥￡1、么2 ￡1、&￡1、13￡1，通过平移（八乂￡1、么丫 ￡1、 A ZE1)和调整椭球长、短轴aE1、bE1，就可以构建出新的工程水准椭球El。具体方法如下：
[0108] 步骤4-1，根据式（14)求出平移后的工程水准椭球E1的空间坐标[XE1 YE1 ZE1]T， 「職'「if
[0109] Ym = AYei + Y A'i」声拟」Lz」（i4)
[0110] 式中，[X Y Z]T表示原坐标系下（如WGS84或ITRF框架）的坐标；[XE1 YE1 ZE1]T表示 工程水准椭球对应坐标系下的坐标;[A X A Y E Z ]T为平移三参数。
[0111] 步骤4-2,采用前述步骤求解得到的aE1、bE1调整椭球参数，在使用中一般通过参数 aEi与扁率eEi表达椭球参数，其中：
(15)
[0113] 经过以上步骤，建立出工程区域水准椭球E1，但由于工程网和城市网的特殊性，导 致了参考椭球坐标到平面坐标及从平面坐标到参考椭球坐标变换的必要性，因此需要将E1 椭球的空间直角坐标转换为大地坐标。
[0114] 步骤5,椭球空间直角坐标到工程平面坐标的转换
[0115] 步骤5-1，E1椭球的空间直角坐标转换为大地坐标
[0116] 将上一步公式（15)所解算的[XE1 YE1 ZE1]T转换为大地坐标，转换公式参照公式 (4)得到各点大地坐标(B，L，H)。
[0120] 步骤5-1，基于El椭球大地坐标向工程平面坐标转换
[0121] 由于上述步骤所确定的椭球参数中椭球的aE1、bE1可变的，而参考椭球的 变带来的直接问题便是大地正反算计算模型中X和Bf参数难于求解。因此需要先建立aE1、b E1 可变的X和Bf参数的通用计算模型，然后进行坐标转换。具体方法如下：
[0122] (1)建立大地子午线弧长)(模型：
[0123] X{B)^c!X = ^iXyiX (16)
[0124] 其中B为大地炜度，X(B)为大地子午线弧长函数，M(B)为大地子午线曲率半径，M (B)为一随着B值而改变的函数。 (18) (17)
[0126]由（16)和（17)可知：
[0128] (2)根据大地子午线弧长)(模型建立垂足炜度Bf模型：
[0129] Bf称为垂足炜度或底点炜度，过待算点P做中央子午线的垂线PP' .0P'弧所对应的 弧段即为Bf。
[0130] M9)
[0131] 公式（19)为XojPBf之间的数学关系，为了求Bf对公式做如下变形：
[0132] 对(19)式两端求微分，则：
[0133] dX〇P=MdBf (20)
[0134] M移到等式的左边，得：
(21 )
[0136] 对(21)式两边同求积分，得
C22)
[0138] 上式即为求解垂足炜度Bf的模型。
[0139] (3)按照大地正算公式计算工程平面坐标：
[0142] 式中：1(B)为子午线弧长，
为卯酉圈半径，t = tan B，€ Aj 为经差，Lo为中央子午线经度，1(B)为从赤道到投影点的椭球面弧长，可用下式计算：1(B) = a[B+0sin2B+ 丫 sin4B+5sin6B+esin8B+.??]。
[0143] 经过以上步骤，建立出工程区域水准椭球E1并进行了坐标转换，该椭球面与工程 测区平均高程面对应的大地水准面较为吻合，使测区大地测量中的几何参考面和物理参考 面统一起来，易于实际的计算。
[0144]本发明以上描述只是一种实施例，但是本发明并不局限于上述的【具体实施方式】。 上述的【具体实施方式】是示意性的，并不是限制性的。凡是采用本发明的材料和方法，在不脱 离本发明宗旨和权利要求所保护的范围情况下，所有具体拓展均属本发明的保护范围之 内。
【主权项】
1. 一种建立工程区域水准椭球的方法，其特征在于，包括以下步骤： 步骤1，确定基于标准椭球的工程水准椭球 在标准椭球方程的基础上引入由aE1、bE1、ΔΧΕ1、ΔΥΕ1、ΔΖ Ε15个待定参数构成的工程区 域水准椭球方程，式中，Χει、Υει、Ζει为在源椭球下测区椭球表面点空间直角坐标，ΔΧει、Δ Υει、Δ Ζει为平 移参数，aEi、bEi分别为长轴和短轴； 步骤2,求解空间直角坐标 利用测区各个控制点高程异常数据，在测区平均高程面基础上，将各个控制点高程换 算到测区平均高程面后所对应的测区平均大地水准面，得到源椭球下测区椭球表面空间直 角坐标 Χει、Yei、Zei ; 步骤3,求解工程椭球参数 在空间直角坐标乂^、￥^、2^下，获取测区5组代表性坐标，然后采用最小二乘原理求解 出5个工程椭球参数Δ ΧΕ1、Δ ΥΕ1、Δ ZE1、aE1、bE1; 步骤4，工程水准椭球坐标转换及椭球参数调整 将5个工程椭球参数代入公式(3)，通过平移（Δ ΧΕ1、Δ ΥΕ1、Δ ZE1)、调整椭球长轴aE1和短 轴bE1，构建出新的工程水准椭球El; 步骤5,椭球空间直角坐标到工程平面坐标的转换 先将E1椭球的空间直角坐标转换为大地坐标，再将其转换到工程平面坐标，得到椭球 面与测区大地水准面密合的工程区域水准椭球。2. 根据权利要求1所述的建立工程区域水准椭球的方法，其特征在于，所述步骤2具体 按照以下步骤实施： 步骤2-1，求解大地高 如果测区中采用GNSS联测IGS站获取了较高精度的ITRF框架下的Xitrf、Yitrf、Z itrf坐标， 则可以得到对应系统下控制点的相应精度的大地高He; 如果测区中没有获取该坐标，可以在GNSS网三维约束平差时采用测站单点定位近似坐 标约束，通过解算得到空间直角坐标后，将各点空间直角坐标以GNSS测量所采用椭球为目 标椭球，利用大地正算公式得到各点位在WGS84或ITRF框架下的大地坐标(B，L，H)，进而得 到大地高He; 步骤2-2，解算参考椭球参数的测区各点对应的大地高 定义控制点相对源椭球的大地高为He，正常高为Hn，则该点的高程异常ξΕ1 = Ηβ-ΗΝ; 设选定的测区平均高程或者根据工程需求而指定的基准面为Η，工程椭球面 拟合所对应的大地高为...Η;'…1? 高程异常为 U -1 ， 二Clflf '·· …d'则工程椭球面拟合所对应的大地高恥=肝 ξΕ1; 求出出后，根据大地坐标与空间直角坐标转换公式将大地坐标Lf /iff转换成 空间直角坐标得到的f 即为空间直角坐标χΕ1、γΕ1、ζΕ1。3. 根据权利要求1所述的建立工程区域水准椭球的方法，其特征在于，所述步骤3具体 按照以下步骤实施： 获取测区5组代表性坐标，由间接平差模型得知： x=X-X° ,l = L-f(X0)=L-L°,F = 5x-/ 非线性误差方程为： V = f{X)-L (：8) 式中：V为观测值的改正数向量（残差向量）；x为包含5个未知向量的估值：将非线性模型(8)式在X*3处用泰勒级数展开，取至一次项，再由X*3为测区已知的5个控制 点求得的第一初始值为： X°=[a° b° ΔΧ° ΔΥ° ΔΖ°]Τ (11)，有η个;Ρ为单位权矩 阵，得到 x={B' PB) 'B1 PI (12) 通过解算式（12)就可以算出i，由此通过式（13)算出X; = 无二卩 S S &? ΔΖ]7 ,、 L 」 （13) 至此得到了待求的未知参数Δ Χει、Δ ΥΕ1、Δ ZEi、aEi、bEi。4. 根据权利要求1所述的建立工程区域水准椭球的方法，其特征在于，步骤5中所述将 大地坐标转换到工程平面坐标具体按照以下步骤实施： (1) 建立大地子午线弧长对莫型：(2) 根据大地子午线弧长)(模型建立垂足炜度Bf模型： Bf称为垂足炜度或底点炜度，过待算点P做中央子午线的垂线PP'.OP'弧所对应的弧段 即为Bf，则有对公式(19)变形后得到即为求解垂足炜度Bf的模型； (3) 按照大地正算公式计算工程平面坐标：式中：1(B)为子午线弧长，为卯酉圈半径，t = tanB夕=1-人为经差， Lo为中央子午线经度，1(B)为从赤道到投影点的椭球面弧长。
【文档编号】G06F17/11GK105930306SQ201610231775
【公开日】2016年9月7日
【申请日】2016年4月14日
【发明人】李祖锋, 赵彦华
【申请人】中国电建集团西北勘测设计研究院有限公司