CN105930306B - 一种建立工程区域水准椭球的方法 - Google Patents

Abstract

本发明公开了一种建立工程区域水准椭球的方法：在标准椭球方程的基础上构建由待定参数构成的工程区域水准椭球方程，利用测区各个控制点高程异常数据，在测区平均高程面基础上，将各个控制点高程换算到测区平均高程面后所对应的测区平均大地水准面，得到源椭球下测区椭球表面点空间直角坐标；在该坐标下求出平移参数和长短轴参数；通过平移和调整椭球长、短轴，构建出新的工程水准椭球E1，并求出目标椭球下的空间坐标，将转换后的坐标在工程水准椭球E1下进行椭球长、短轴可变的大地坐标转换，进而转换到工程平面坐标，即可用于进一步的精准归算。该方法构建的水准椭球面与工程测区平均高程面对应的大地水准面较好吻合，利于数据的精准归算。

Description

一种建立工程区域水准椭球的方法

技术领域

本发明属于地理测绘科学与技术领域，具体涉及一种建立工程区域水准椭球的方法。

背景技术

工程椭球的理想状态是其参考椭球面与测区大地水准面重合，即，椭球面是正常重力场中的一个等位面。工程参考椭球面是工程测量计算的基准面，又是研究地球形状和地图投影的参考面，所以工程椭球一定要满足大地参考系统的基本常数(a,b,GM,w)，而地球的重力等位面是一个不规则的曲面，其不满足大地参考系统的基本常数，因此，不存在既满足大地参考系统基本常数、又满足表面重力等位条件的参考椭球。

目前，工程中习惯继承通用椭球模型参数，如1980西安坐标系及CGCS2000坐标系所采用的椭球分别就是参心与地心椭球，对椭球模型长半轴进行简单调整，建立工程椭球模型，这种参数的简单调整，远无法达到与工程区域大地水准面最为密合的目的，同时受方法的限制，有着较大的误差。

发明内容

本发明提供了一种建立工程区域水准椭球的方法，并构建出与测区在已有数据下最优密合的椭球参数，解决了现有工程椭球无法与整个工程区域大地水准面较好密合的问题。

本发明采用的技术方案是，一种建立工程区域水准椭球的方法，包括以下步骤：

步骤1，确定基于标准椭球的工程水准椭球

在标准椭球方程的基础上引入由a_E1、b_E1、ΔX_E1、ΔY_E1、ΔZ_E15个待定参数构成的工程区域水准椭球方程，

式中，X_E1、Y_E1、Z_E1为在源椭球下测区椭球表面点空间直角坐标，ΔX_E1、ΔY_E1、ΔZ_E1为平移参数，a_E1、b_E1分别为长轴和短轴；

步骤2，求解空间直角坐标和工程椭球参数

利用测区各个控制点高程异常数据，在测区平均高程面基础上，将各个控制点高程换算到测区平均高程面后所对应的测区平均大地水准面，得到源椭球下测区椭球表面点空间直角坐标X_E1、Y_E1、Z_E1；

步骤3，求解工程椭球参数

在空间直角坐标X_E1、Y_E1、Z_E1下，获取测区5组代表性坐标，然后采用最小二乘原理求解出5个工程椭球参数ΔX_E1、ΔY_E1、ΔZ_E1、a_E1、b_E1；

步骤4，工程水准椭球坐标转换及椭球参数调整

将5个工程椭球参数代入公式(3)，通过平移(ΔX_E1、ΔY_E1、ΔZ_E1)、调整椭球长轴a_E1和短轴b_E1，构建出新的工程水准椭球E1；

步骤5，椭球空间直角坐标到工程平面坐标的转换

先将E1椭球的空间直角坐标转换为大地坐标，再将其转换到工程平面坐标，得到椭球面与测区大地水准面密合的工程区域水准椭球。

本发明的特点还在于：

上述步骤2具体按照以下步骤实施：

步骤2-1，求解大地高

如果测区中采用GNSS联测IGS站获取了较高精度的ITRF框架下的X_itrf、Y_itrf、Z_itrf坐标，则可以得到对应系统下控制点的相应精度的大地高H_e。

如果测区中没有获取该坐标，可以在GNSS网三维约束平差时采用测站单点定位近似坐标约束，通过解算得到空间直角坐标后，将各点空间直角坐标以GNSS测量所采用椭球为目标椭球，利用大地正算公式得到各点位在WGS84或ITRF框架下的大地坐标(B，L，H)，进而得到大地高H_e；

步骤2-2，解算参考椭球参数的测区各点对应的大地高

定义控制点相对源椭球的大地高为H_e，正常高为H_N，则该点的高程异常ξ^E1＝H_e-H_N；

设选定的测区平均高程或者根据工程需求而指定的基准面为H，工程椭球面拟合所对应的大地高为:高程异常为则工程椭球面拟合所对应的大地高H_G＝H+ξ^E1；

求出H_G后，根据大地坐标与空间直角坐标转换公式将大地坐标转换成空间直角坐标得到的即为空间直角坐标X_E1、Y_E1、Z_E1。

上述步骤3具体按照以下步骤实施：

获取测区5组代表性坐标，由间接平差模型得知：

非线性误差方程为：

式中：V为观测值的改正数向量(残差向量)；为包含5个未知向量的估值：

将非线性模型(8)式在X⁰处用泰勒级数展开，取至一次项，再由X⁰为测区已知的5个控制点求得的第一初始值为：

X⁰＝[a⁰ b⁰ ΔX⁰ ΔY⁰ ΔZ⁰]^T (11)

l为带入X⁰后式子的值，有n个；P为单位权矩阵，得到

通过解算式(12)就可以算出由此通过式(13)算出

至此得到了待求的未知参数ΔX_E1、ΔY_E1、ΔZ_E1、a_E1、b_E1。

上述步骤5中将大地坐标转换到工程平面坐标具体按照以下步骤实施：

(1)建立大地子午线弧长X模型：

(2)根据大地子午线弧长X模型建立垂足纬度B_f模型：

B_f称为垂足纬度或底点纬度，过待算点P做中央子午线的垂线PP'.OP'弧所对应的弧段即为B_f，则有

对公式(19)变形后得到

即为求解垂足纬度B_f的模型；

(3)按照大地正算公式计算工程平面坐标：

式中：l(B)为子午线弧长，为卯酉圈半径，t＝tan B，为经差，L₀为中央子午线经度，l(B)为从赤道到投影点的椭球面弧长。

本发明的有益效果是：本发明的建立工程区域水准椭球的方法构建的工程区域水准椭球面与工程测区平均高程面对应的大地水准面较好吻合，能使测区大地测量中的几何参考面和物理参考面统一起来，易于实际的计算。与现有其他方法构建的可保证在选定高程面上，椭球面与测区大地水准面的最好密合，利于数据的精准归算。

具体实施方式

下面结合具体实施方式对本发明作进一步的详细说明。

步骤1，确定基于标准椭球的工程水准椭球

基于全球性通用椭球，通过改变椭球中心位置及椭球长短轴，使得椭球面与测区平均高程面对应的大地水准面最优吻合。工程水准椭球的构建必须首先确定测区平均高程面所对应的大地水准面，该大地水准面可以采用分布于测区的若干具有正常高成果的控制点近似表达出来。

标准椭球方程见式(1)。

式中，a、b为赤道半径，c为极半径，X、Y、Z为标准椭球的空间直角坐标。

按照空间直角坐标系的定义，旋转椭球体在X和Y方向半径相等，故a＝b，得到：

理想的工程水准椭球体是椭球面与测区平均大地水准面重合，但受椭球数学表达的限制，椭球面与测区平均大地水准面重合无法实现，但可以在一定规则下达到椭球面与大地水准面的最优吻合。因此，在公式(2)的基础上引入平移参数，形成了由长轴、短轴、平移参数构成的椭球方程，于是式(2)可表达为：

式中，X_E1、Y_E1、Z_E1为在源椭球下测区椭球表面点空间直角坐标，ΔX_E1、ΔY_E1、ΔZ_E1为平移参数，a_E1、b_E1分别为长轴和短轴。其中，ΔX_E1、ΔY_E1、ΔZ_E1、a_E1、b_E1为五个待定工程椭球参数。

步骤2，求解空间直角坐标X_E1、Y_E1、Z_E1

X_E1、Y_E1、Z_E1可以理解为控制点在所求工程椭球表面投影点的位置在源椭球下的空间直角坐标。其求解方法为，利用测区各个控制点高程异常数据，在测区平均高程面基础上，将各个控制点高程换算到测区平均高程面后所对应的测区平均大地水准面，得到源椭球下测区椭球表面点空间直角坐标X_E1、Y_E1、Z_E1。具体方法如下：

步骤2-1，求解大地坐标

如果测区中采用GNSS联测IGS站获取了较高精度的ITRF框架下的X_itrf、Y_itrf、Z_itrf坐标，则可以得到对应系统下控制点的相应精度的大地高H_e。

如果测区中没有获取该坐标，可以在GNSS网三维约束平差时采用测站单点定位近似坐标约束，通过解算得到空间直角坐标后，将各点空间直角坐标以GNSS测量所采用椭球为目标椭球，利用公式(4)得到点各位在WGS84或ITRF框架下的大地坐标(B，L，H)。

式中，N表示卯酉圈的半径。

由式(4)得到控制点的经度B、纬度L和大地高H_e。

步骤2-2，解算参考椭球参数的测区各点对应大地高的确定方法：

定义控制点相对源椭球的大地高为H_e，正常高为H_N，则该点的高程异常ξ^E1求解见式(5)。

ξ^E1＝H_e-H_N (5)

原则上用于求定参考椭球参数的测区各点的大地高宜采用测区的平均高程或者选定的某一高程面，但由于各点的高程异常不一致，宜在选定的高程基础上加上对应点位的高程异常值。

设选定的测区平均高程或者根据工程需求而指定的基准面为H，工程椭球面拟合所对应的大地高为:高程异常为则工程椭球面拟合所对应的大地高由式(6)得到

H_G＝H+ξ^E1 (6)

求出H_G后，将大地坐标按照式(7)转换成空间直角坐标

式中，N为卯酉圈的半径，

式(7)得到的即为空间直角坐标X_E1、Y_E1、Z_E1。

步骤3、求解工程椭球参数

式(3)的五个待定工程椭球参数分别为ΔX_E1、ΔY_E1、ΔZ_E1、a_E1、b_E1，要求解该方程，至少需要均匀分布于测区的5组代表性坐标，然后采用最小二乘原理求解该5个参数。具体方法如下：

步骤3-1，由间接平差模型得知：

非线性误差方程为：

式中：V为观测值的改正数向量(残差向量)；为包含5个未知向量的估值：

由非线性误差方程式(8)知，式中至少有5个方程，且有5个未知数(3个坐标平移参数以及椭球长、短轴参数)。因此非线性误差方程式(8)是非线性不定方程组，有无穷组解。在这无穷组解中，必然有一组解能使

推出：

将非线性模型(8)式在X⁰处用泰勒级数展开，取至一次项，得：

由X⁰为测区已知的5个控制点求得的第一初始值为：

X⁰＝[a⁰ b⁰ ΔX⁰ ΔY⁰ ΔZ⁰]^T (11)

l为带入X⁰后式子的值，有n个。P为单位权矩阵。得到

通过解算式(12)就可以算出由此通过式(13)算出

至此得到了待求的未知参数ΔX_E1、ΔY_E1、ΔZ_E1、a_E1、b_E1。

通过这样转换的参数对空间三维坐标进行ΔX、ΔY、ΔZ平移转换及对椭球参数进行变换，可以保证椭球面与测区平均大地水准面的较好吻合，可满足常见的工程精密坐标归算及高程拟合需求。

步骤4，工程水准椭球坐标转换及椭球参数调整

根据前述步骤求解的平移参数ΔX_E1、ΔY_E1、ΔZ_E1、a_E1、b_E1，通过平移(ΔX_E1、ΔY_E1、ΔZ_E1)和调整椭球长、短轴a_E1、b_E1，就可以构建出新的工程水准椭球E1。具体方法如下：

步骤4-1，根据式(14)求出平移后的工程水准椭球E1的空间坐标[X_E1 Y_E1 Z_E1]^T，

式中，[X Y Z]^T表示原坐标系下(如WGS84或ITRF框架)的坐标；[X_E1 Y_E1 Z_E1]^T表示工程水准椭球对应坐标系下的坐标；[ΔX ΔY ∑Z]^T为平移三参数。

步骤4-2，采用前述步骤求解得到的a_E1、b_E1调整椭球参数，在使用中一般通过参数a_E1与扁率e_E1表达椭球参数，其中：

经过以上步骤，建立出工程区域水准椭球E1，但由于工程网和城市网的特殊性，导致了参考椭球坐标到平面坐标及从平面坐标到参考椭球坐标变换的必要性，因此需要将E1椭球的空间直角坐标转换为大地坐标。

步骤5，椭球空间直角坐标到工程平面坐标的转换

步骤5-1，E1椭球的空间直角坐标转换为大地坐标

将上一步公式(15)所解算的[X_E1 Y_E1 Z_E1]^T转换为大地坐标，转换公式参照公式(4)得到各点大地坐标(B，L，H)。

其中：

步骤5-1，基于E1椭球大地坐标向工程平面坐标转换

由于上述步骤所确定的椭球参数中椭球的a_E1、b_E1可变的，而参考椭球的a_E1、b_E1可变带来的直接问题便是大地正反算计算模型中X和B_f参数难于求解。因此需要先建立a_E1、b_E1可变的X和B_f参数的通用计算模型，然后进行坐标转换。具体方法如下：

(1)建立大地子午线弧长X模型：

其中B为大地纬度，X(B)为大地子午线弧长函数，M(B)为大地子午线曲率半径，M(B)为一随着B值而改变的函数。

由(16)和(17)可知:

(2)根据大地子午线弧长X模型建立垂足纬度B_f模型：

B_f称为垂足纬度或底点纬度，过待算点P做中央子午线的垂线PP'.OP'弧所对应的弧段即为B_f。

公式(19)为X_op和B_f之间的数学关系，为了求B_f对公式做如下变形：

对(19)式两端求微分,则：

dX_op＝MdB_f (20)

M移到等式的左边,得：

对(21)式两边同求积分,得

上式即为求解垂足纬度B_f的模型。

(3)按照大地正算公式计算工程平面坐标：

式中：l(B)为子午线弧长，为卯酉圈半径，t＝tan B，为经差，L₀为中央子午线经度，l(B)为从赤道到投影点的椭球面弧长，可用下式计算：l(B)＝α[B+βsin2B+γsin4B+δsin6B+εsin8B+...]。

经过以上步骤，建立出工程区域水准椭球E1并进行了坐标转换，该椭球面与工程测区平均高程面对应的大地水准面较为吻合，使测区大地测量中的几何参考面和物理参考面统一起来，易于实际的计算。

本发明以上描述只是一种实施例，但是本发明并不局限于上述的具体实施方式。上述的具体实施方式是示意性的，并不是限制性的。凡是采用本发明的材料和方法，在不脱离本发明宗旨和权利要求所保护的范围情况下，所有具体拓展均属本发明的保护范围之内。

Claims (4)

1.一种建立工程区域水准椭球的方法，其特征在于，包括以下步骤：

步骤1，确定基于标准椭球的工程水准椭球

在标准椭球方程的基础上引入由a_E1、b_E1、ΔX_E1、ΔY_E1、ΔZ_E15个待定参数构成的工程区域水准椭球方程，

式中，X_E1、Y_E1、Z_E1为在源椭球下测区椭球表面点空间直角坐标，ΔX_E1、ΔY_E1、ΔZ_E1为平移参数，a_E1、b_E1分别为长轴和短轴；

步骤2，求解空间直角坐标

利用测区各个控制点高程异常数据，在测区平均高程面基础上，将各个控制点高程换算成测区平均正常高程所对应的各个控制点大地高，得到源椭球下测区椭球表面空间直角坐标X_E1、Y_E1、Z_E1；

步骤3，求解工程椭球参数

在空间直角坐标X_E1、Y_E1、Z_E1下，获取测区5组代表性坐标，然后采用最小二乘原理求解出5个工程椭球参数ΔX_E1、ΔY_E1、ΔZ_E1、a_E1、b_E1；

步骤4，工程水准椭球坐标转换及椭球参数调整

将5个工程椭球参数代入公式(3)，通过平移(ΔX_E1、ΔY_E1、ΔZ_E1)、调整椭球长轴a_E1和短轴b_E1，构建出新的工程水准椭球E1；

步骤5，椭球空间直角坐标到工程平面坐标的转换

先将E1椭球的空间直角坐标转换为大地坐标，再将其转换到工程平面坐标，得到椭球面与测区大地水准面密合的工程区域水准椭球。

2.根据权利要求1所述的建立工程区域水准椭球的方法，其特征在于，所述步骤2具体按照以下步骤实施：

步骤2-1，求解大地高

如果测区中采用GNSS联测IGS站获取了较高精度的ITRF框架下的X_itrf、Y_itrf、Z_itrf坐标，则可以得到对应系统下控制点的相应精度的大地高H_e；

如果测区中没有获取该坐标，可以在GNSS网三维约束平差时采用测站单点定位近似坐标约束，通过解算得到空间直角坐标后，将各点空间直角坐标以GNSS测量所采用椭球为目标椭球，利用大地正算公式得到各点位在WGS84或ITRF框架下的大地坐标(B，L，H)，进而得到大地高H_e；

步骤2-2，解算参考椭球参数的测区各点对应的大地高

定义控制点相对源椭球的大地高为H_e，正常高为H_N，则该点的高程异常ξ^E1＝H_e-H_N；

设选定的测区平均高程或者根据工程需求而指定的基准面为H，工程椭球面拟合所对应的大地高为:高程异常为则工程椭球面拟合所对应的大地高H_G＝H+ξ^E1；

求出H_G后，根据大地坐标与空间直角坐标转换公式将大地坐标转换成空间直角坐标得到的即为空间直角坐标X_E1、Y_E1、Z_E1。

3.根据权利要求1所述的建立工程区域水准椭球的方法，其特征在于，所述步骤3具体按照以下步骤实施：

获取测区5组代表性坐标，由间接平差模型得知：

l＝L-f(X⁰)＝L-L⁰，

非线性误差方程为：

式中：V为观测值的改正数向量(残差向量)；为包含5个未知向量的估值：

将非线性模型(8)式在X⁰处用泰勒级数展开，取至一次项，再由X⁰为测区已知的5个控制点求得的第一初始值为：

X⁰＝[a⁰ b⁰ ΔX⁰ ΔY⁰ ΔZ⁰]^T (11)

l为带入X⁰后式子的值，有n个；P为单位权矩阵，得到

通过解算式(12)就可以算出由此通过式(13)算出

至此得到了待求的未知参数ΔX_E1、ΔY_E1、ΔZ_E1、a_E1、b_E1。

4.根据权利要求1所述的建立工程区域水准椭球的方法，其特征在于，步骤5中所述将大地坐标转换到工程平面坐标具体按照以下步骤实施：

(1)建立大地子午线弧长X模型：

(2)根据大地子午线弧长X模型建立垂足纬度B_f模型：

B_f称为垂足纬度或底点纬度，过待算点P做中央子午线的垂线PP'.OP'弧所对应的弧段即为B_f，则有

对公式(19)变形后得到

即为求解垂足纬度B_f的模型；

(3)按照大地正算公式计算工程平面坐标：

式中：为卯酉圈半径，t＝tanB，为经差，L₀为中央子午线经度。