留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

利用Molodensky理论求解第二大地边值问题

马健 魏子卿

马健, 魏子卿. 利用Molodensky理论求解第二大地边值问题[J]. 武汉大学学报 ● 信息科学版, 2019, 44(10): 1478-1483. doi: 10.13203/j.whugis20170420
引用本文: 马健, 魏子卿. 利用Molodensky理论求解第二大地边值问题[J]. 武汉大学学报 ● 信息科学版, 2019, 44(10): 1478-1483. doi: 10.13203/j.whugis20170420
MA Jian, WEI Ziqing. The Second Geodetic Boundary Value Problem Based on Molodensky Theory[J]. Geomatics and Information Science of Wuhan University, 2019, 44(10): 1478-1483. doi: 10.13203/j.whugis20170420
Citation: MA Jian, WEI Ziqing. The Second Geodetic Boundary Value Problem Based on Molodensky Theory[J]. Geomatics and Information Science of Wuhan University, 2019, 44(10): 1478-1483. doi: 10.13203/j.whugis20170420

利用Molodensky理论求解第二大地边值问题

doi: 10.13203/j.whugis20170420
基金项目: 

国家自然科学基金 41674025

国家自然科学基金 41674082

地理信息工程国家重点实验室开放研究基金 SKLGIE2016-M-1-5

详细信息
    作者简介:

    马健, 博士生, 主要从事物理大地测量研究。majian_19881006@163.com

  • 中图分类号: P223

The Second Geodetic Boundary Value Problem Based on Molodensky Theory

Funds: 

The National Natural Science Foundation of China 41674025

The National Natural Science Foundation of China 41674082

the Open Research Foundation of State Key Laboratory of Geo-information Engineering SKLGIE2016-M-1-5

More Information
    Author Bio:

    MA Jian, PhD candidate, specializes in the theories and methods of physical geodesy. E-mail: majian_19881006@163.com

图(5) / 表(3)
计量
  • 文章访问数:  764
  • HTML全文浏览量:  84
  • PDF下载量:  116
  • 被引次数: 0
出版历程
  • 收稿日期:  2018-08-19
  • 刊出日期:  2019-10-05

利用Molodensky理论求解第二大地边值问题

doi: 10.13203/j.whugis20170420
    基金项目:

    国家自然科学基金 41674025

    国家自然科学基金 41674082

    地理信息工程国家重点实验室开放研究基金 SKLGIE2016-M-1-5

    作者简介:

    马健, 博士生, 主要从事物理大地测量研究。majian_19881006@163.com

  • 中图分类号: P223

摘要: 过去由于无法获得大地高数据,传统的第三大地边值问题采用重力异常作为边值条件。GNSS技术的发展为第二边值问题的研究带来了契机。研究比较成熟的第三边值理论无疑为第二边值问题提供了很好的参考和借鉴,对此开展将第三边值问题中计算似大地水准面的Molodensky理论方法应用于第二边值问题的研究。首先推导了Hotine算子与梯度算子的关系,然后给出了基于Molodensky理论求解第二边值问题的算法。实验结果表明,该算法与传统第三边值问题中Molodensky理论的边值解精度相当,说明基于Molodensky理论求解第二大地边值问题是完全可行的。

English Abstract

马健, 魏子卿. 利用Molodensky理论求解第二大地边值问题[J]. 武汉大学学报 ● 信息科学版, 2019, 44(10): 1478-1483. doi: 10.13203/j.whugis20170420
引用本文: 马健, 魏子卿. 利用Molodensky理论求解第二大地边值问题[J]. 武汉大学学报 ● 信息科学版, 2019, 44(10): 1478-1483. doi: 10.13203/j.whugis20170420
MA Jian, WEI Ziqing. The Second Geodetic Boundary Value Problem Based on Molodensky Theory[J]. Geomatics and Information Science of Wuhan University, 2019, 44(10): 1478-1483. doi: 10.13203/j.whugis20170420
Citation: MA Jian, WEI Ziqing. The Second Geodetic Boundary Value Problem Based on Molodensky Theory[J]. Geomatics and Information Science of Wuhan University, 2019, 44(10): 1478-1483. doi: 10.13203/j.whugis20170420
  • Stokes理论和Molodensky理论是以重力异常为基本物理输入量、使用正高/正常高高程系统的第三大地边值问题的两大理论体系。在第三边值问题中,重力异常为大地水准面上的重力与参考椭球面上的正常重力之差(Stokes理论)或地表重力与近似地形面上的正常重力之差(Molodensky理论)。1969年,英国学者Hotine提出了Hotine积分公式,提出了以扰动重力为边值条件的第二大地边值问题。由于扰动重力是同一点的重力与正常重力之差,因此扰动重力较重力异常具有更明确的物理意义。

    过去由于无法获得大地高数据,因此第二大地边值理论的发展和应用受到限制。现阶段,SRTM(shuttle radar topography mission)数据在确定高精度的(似)大地水准面中获得越来越广泛的应用[1-4]。由于SRTM数据以EGM96(earth gravitational model 96)大地水准面起伏为高程基准,因此通过EGM96重力场模型可方便地将SRTM数据转化为大地高。与此同时,卫星重力技术的发展使重力场模型的精度不断提高,因而可借助高精度的重力场模型将离散重力点的正常高数据也转化为大地高。SRTM数据和高精度重力场模型的出现,使得大地高数据能够间接地获得,因此也初步具备了研究与应用第二大地边值问题的条件。

    在第二边值问题的理论层面,Moritz最先提出GPS边值问题及其一般解式。文献[5]分析了重力异常和扰动重力的特点,指出以扰动重力作为边值条件更加简单、严密。文献[6]详细论述了以地心参考椭球面为边界面的第二大地边值问题的理论意义和实际意义,进而引出了该问题的原理和内涵。文献[7-12]给出了GPS边值问题的定义、数学表达,并通过类似于Molodensky理论中的分析方法,得到了第二边值问题以地表为参考面的实用公式。在上述研究成果的基础上,本文研究将第三边值问题中发展和应用均比较成熟的Molodensky理论方法应用于第二边值问题。为此,本文推导了扰动重力梯度与重力异常梯度的关系、Hotine算子[13]与梯度算子的关系,给出了针对第二边值问题的Molodensky算法,并通过实验验证在第二边值问题中应用Molodensky算法的可行性。

    • 传统第三边值问题中,用于解算似大地水准面的Molodensky算法公式为:

      $$ \begin{array}{l} {\zeta _A}{\rm{ = }}\frac{R}{{4{\rm{ \mathsf{ π} }}{\gamma _A}}}\iint {\left( {\Delta g + C} \right)} S\left( \psi \right)d\mathit{\Omega }'-\\ {\rm{ }} \frac{{{h_A}\Delta {g_A}}}{{{\gamma _A}}} + \frac{{{\rm{ \mathsf{ π} }}G\rho }}{{{\gamma _A}}}\left( {h_A^2 - \delta h_A^2} \right) \end{array} $$ (1)

      式中,dΩ′为球面微分单元;ζAγAhA、ΔgA分别为计算点A处的高程异常、正常重力、高程、重力异常;R为地球平均半径;ψ为积分点和计算点间的角距;Δg为Stokes积分中流动积分点的重力异常;C为局部地形改正;S(ψ)为Stokes积分的核函数,定义为:

      $$ \begin{array}{l} S(\psi ){\rm{ = }}\sum\limits_{n = 2}^\infty {\frac{{2n + 1}}{{n - 1}}{P_n}\left( {\cos \psi } \right)} {\rm{ = }}\frac{1}{t} - 6t - 4+\\ ~~~~~~~~~~~~~ 10{t^2} - 3\left( {1 - 2{t^2}} \right)\ln \left( {t + {t^2}} \right){\rm{ }} \end{array} $$ (2)

      其中,t=sin(ψ/2);Pn(cosψ)为勒让德函数;δhA2hA2(计算点高程的平方)的零阶项与一阶项之和,通常采用下式计算:

      $$ \begin{array}{l} \delta h_A^2 = 0.453 - 0.018\sin {\varphi _A} + \\ 0.087\cos {\varphi _A}\cos {\lambda _A} + 0.204\cos {\varphi _A}\sin {\lambda _A} \end{array} $$ (3)

      其中,φAλA分别为计算点的球面纬度和经度。

    • 类似于Stokes积分的Moritz解析延拓解的形式,李斐等[7-9]给出了以地面为参考面的Hotine积分的实用公式:

      $$ \begin{array}{c} {\zeta _A}{\rm{ = }}\frac{R}{{4{\rm{ \mathsf{ π} }}{\gamma _A}}} \times \\ \iint {\left[ {\delta g - \left( {h - {h_A}} \right)L\left( {\delta g} \right)} \right]} H\left( \psi \right)d\mathit{\Omega }' \end{array} $$ (4)

      式中,δg为Hotine积分中流动积分点的扰动重力;h为Hotine积分中流动积分点的高程;H(ψ)为Hotine积分的核函数;L表示梯度算子;L(δg)为地面扰动重力沿地心向径的梯度,其计算公式为:

      $$ \begin{gathered} L\left( {\delta g} \right) = \frac{{\partial \delta g}}{{\partial r}} =\hfill \\ \frac{{{R^2}}}{{2{{\rm{ \mathsf{ π} }}}}}\iint\limits_{{\mathit{\Omega }_0}} {\frac{1}{{\ell _0^3}}\left( {\delta g' - \delta g} \right){\text{d}}\mathit{\Omega }'} - \frac{{2\delta g}}{R} \hfill \\ \end{gathered} $$ (5)

      式中,δg′为梯度算子中流动积分点的扰动重力;r为计算点的地心向径;R为地球平均半径;ℓ0为计算点与积分点的平面距离,计算公式为ℓ0=2Rsin(ψ/2)。

      Hotine积分核函数H(ψ)在实际应用中有多种形式,其中包含零阶、一阶项的标准形式为:

      $$ \begin{gathered} H\left( \psi \right) = \sum\limits_{n = 0}^\infty {\frac{{2n + 1}}{{n + 1}}{P_n}\left( {\cos \psi } \right)}= \hfill \\ \frac{1}{t} - \ln \left( {1 + \frac{1}{t}} \right) \hfill \\ \end{gathered} $$ (6)

      式(6)所示的标准Hotine核函数的Legendre级数包含了零阶、一阶项。由于通常认为扰动位不含零阶、一阶项,因此标准Hotine核函数也可以是如下不含零阶、一阶项的形式:

      $$ \begin{gathered} H\left( \psi \right) = \sum\limits_{n = 2}^\infty {\frac{{2n + 1}}{{n + 1}}{P_n}\left( {\cos \psi } \right)}= \hfill \\ \frac{1}{t} - \ln \left( {1 + \frac{1}{t}} \right) - 1 - \frac{3}{2}\cos \psi \hfill \\ \end{gathered} $$ (7)

      此外还可将Hotine核函数除去与参考的重力场模型阶数一致的Legendre级数项,可称为截断Hotine核函数:

      $$ \begin{gathered} H\left( \psi \right) = \sum\limits_{n = L + 1}^\infty {\frac{{2n + 1}}{{n + 1}}{P_n}\left( {\cos \psi } \right)}= \hfill \\ \frac{1}{t} - \ln \left( {1 + \frac{1}{t}} \right) - \sum\limits_{n = 0}^L {\frac{{2n + 1}}{{n + 1}}{P_n}\left( {\cos \psi } \right)} \hfill \\ \end{gathered} $$ (8)

      式中,L为移去Legendre级数的最高阶数。

    • 下面给出第二边值问题中的Molodensky算法的推导过程。

      球近似下,重力异常和扰动重力与扰动位T的关系分别为:

      $$ \left\{\begin{array}{l}{\Delta g=-\frac{\partial T}{\partial r}-\frac{2}{R} T} \\ {\delta g=-\frac{\partial T}{\partial r}}\end{array}\right. $$ (9)

      由式(9)可知Δg=δg-2T/R,进而可推得扰动重力梯度与重力异常梯度存在下列关系式:

      $$ L\left( {\delta g} \right) = L\left( {\Delta g} \right) + 2\delta g/R $$ (10)

      式中,2δg/R量级很小,可忽略(梯度算子式(5)也包含2δg/R项,该项在梯度算子中通常也被忽略),即有L(δg)≈Lg)。在传统的Molodensky理论中,通过假设重力异常和高程呈线性关系,得到如下公式[14]

      $$ hL\left( {\Delta g} \right) \approx - C + {{\rm{ \mathsf{ π} }}}G\rho L\left( {{h^2}} \right) $$ (11)

      因此有hL(δg)≈-CL(h2),将其代入式(4)得:

      $$ \begin{gathered} \zeta_{A}= \frac{R}{4 {\rm{ \mathsf{ π} }} \gamma_{A}} \iint\left[\delta g+C-{\rm{ \mathsf{ π} }} G_{\rho} L\left(h^{2}\right)+\right.\\ h_{A} L(\delta g) ] H(\psi) \mathrm{d} \mathit{\Omega }^{\prime} \end{gathered} $$ (12)

      下面推导Hotine算子与梯度算子的关系式。对于球面上的连续函数f,设g为其梯度,则有:

      $$ g = L\left( f \right) = \sum\limits_{n = 0}^\infty {{g_n}} $$ (13)

      gngn阶球谐项,根据球谐展开原理,有[14]:

      $$ {g_n} = - \frac{1}{R}\left( {n + 1} \right){f_n} $$ (14)

      式中,fnfn阶球谐项。根据球函数的展开理论,L(f)的Hotine积分可进行如下推导:

      $$ \begin{gathered} {\text{ }}HL\left( f \right) = H\left( g \right){\text{ = }} \hfill \\ \frac{R}{{4{{\rm{ \mathsf{ π} }}}\gamma }}\iint g\sum\limits_{n = 0}^\infty {\frac{{2n + 1}}{{n + 1}}{P_n}\left( {\cos \psi } \right)} {\text{d}}\mathit{\Omega }' = \hfill \\ \frac{R}{{4{{\rm{ \mathsf{ π} }}}\gamma }}\sum\limits_{n = 0}^\infty {\frac{{2n + 1}}{{n + 1}}} \iint g{P_n}\left( {\cos \psi } \right){\text{d}}\mathit{\Omega }' = \hfill \\ \frac{R}{\gamma }\sum\limits_{n = 0}^\infty {\frac{{{g_n}}}{{n + 1}}} = - \frac{1}{\gamma }\sum\limits_{n = 0}^\infty {{f_n}} = - \frac{1}{\gamma }f \hfill \\ \end{gathered} $$ (15)

      式中,H表示包含零阶、一阶项的标准Hotine算子(见式(6))。式(15)的结果说明,当Hotine积分核为包含零阶、一阶项的标准形式时,Hotine算子与梯度算子严格满足关系式HL(f)=-γ-1I(f)(I表示单位算子,I(f)=f)。当Hotine核函数为式(7)中不含零阶、一阶项的形式时,根据式(15)可很容易得出:

      $$ HL\left( f \right) = - \frac{1}{\gamma }\sum\limits_{n = 2}^\infty {{f_n}} = - \frac{1}{\gamma }\left( {f - {f_0} - {f_1}} \right) $$ (16)

      式中,f0f1分别为f的零阶和一阶球谐项。

      在传统的Molodensky理论中推导了梯度算子与Stokes算子间的关系式为:

      $$ SL\left( f \right) = - \frac{1}{\gamma }\left( {f - {f_0} - {f_1}} \right) - \frac{{2S\left( f \right)}}{R} $$ (17)

      式中,S表示Stokes算子,该式等号右端最后一项量级很小,式(1)实际上省略了该项的计算。比较式(15)、(16)与式(17)可知,Hotine算子与梯度算子的关系式无需忽略微小项,因而比传统Molodensky理论中采用的Stokes算子与梯度算子的近似关系式SL(f)=-γ-1I(f-f0-f1)更为严格准确。

      将式(15)代入式(12),得到Hotine积分核包括零阶、一阶项时的Molodensky算法为:

      $$ \begin{gathered} {\zeta _A}{\text{ = }}\frac{R}{{4{{\rm{ \mathsf{ π} }}}{\gamma _A}}}\iint {\left( {\delta g + C} \right)} H\left( \psi \right){\text{d}}\mathit{\Omega }' - \hfill \\ {\text{ }}\frac{{{h_A}\delta {g_A}}}{{{\gamma _A}}} + \frac{{{{\rm{ \mathsf{ π} }}}G\rho }}{{{\gamma _A}}}h_A^2 \hfill \\ \end{gathered} $$ (18)

      将式(16)代入式(12),即为Hotine积分核不包括零阶、一阶项时的Molodensky算法:

      $$ \begin{gathered} {\zeta _A}{\text{ = }}\frac{R}{{4{{\rm{ \mathsf{ π} }}}{\gamma _A}}}\iint {\left( {\delta g + C} \right)} H\left( \psi \right){\text{d}}\Omega' - \hfill \\ {\text{ }}\frac{{{h_A}\delta {g_A}}}{{{\gamma _A}}} + \frac{{{{\rm{ \mathsf{ π} }}}G\rho }}{{{\gamma _A}}}\left( {h_A^2 - \delta h_A^2} \right) \hfill \\ \end{gathered} $$ (19)

      式中,δhA2的计算参见式(3)。类似于传统的Molodensky理论,式(19)的结果也是基于重力异常不包含零阶和一阶项这一前提。

    • SRTM数据以EGM96地球重力场模型计算的大地水准面为高程基准面。为了分别进行第二、第三边值实验,本实验将SRTM3数据通过以下方式分别转化为大地高和正常高:①将SRTM3正高数据与EGM96大地水准面高相加以得到大地高数据;②通过EGM96模型的高程异常——大地水准面高改正模型(改正模型的球谐系数[15]与EGM96模型同时发布),将SRTM3正高数据转化为正常高数据。

      实验区实测重力离散点共有70 379个,范围为105.5°~116.5°E,25.5°~34.5°N。将离散重力点格网化为106°~116°E、26°~34°N区域范围的1.5′×1.5′分辨率的扰动重力和重力异常数据。该区域地形高范围为0~3 370 m,平均515 m,图 1为区域地形分布图。

      图  1  区域地形图

      Figure 1.  Topographic Map of the Experiment Zone

      表 1对格网扰动重力和重力异常的最小值、最大值、平均值及均方根(root mean square, RMS)进行了统计,其分布如图 2图 3所示。

      表 1  格网扰动重力和重力异常统计/mGal

      Table 1.  Statistics of the Grid Gravity Disturbances and Gravity Anomalies/mGal

      统计量 min max mean RMS
      扰动重力 -141.70 194.74 -28.21 41.97
      重力异常 -134.69 194.58 -21.44 36.89

      图  2  扰动重力分布

      Figure 2.  Distribution of the Gravity Disturbances

      图  3  重力异常分布

      Figure 3.  Distribution of the Gravity Anomalies

      表 1中扰动重力和重力异常的统计结果相近,反映了该两类物理量具有相近的频谱特性。图 2扰动重力的分布与图 3重力异常的分布比较相似,说明本文推导的扰动重力梯度和重力异常梯度近似相等的关系式与实际情况是比较符合的。

      下面以Eigen-6C4重力场模型的前360阶为参考重力场进行边值解算实验。生成的重力似大地水准面范围为108°~114°E、28°~32°N。该区域GPS水准检核点共计68个,其分布如图 4所示。

      图  4  GPS水准点分布

      Figure 4.  Distribution of the GPS Levelling Points

      首先利用Moritz提出的梯度算子计算Hotine积分,Hotine积分核分别采用以下3种形式:①包含零阶和一阶项的标准形式;②不含零阶和一阶项的标准形式;③截断形式,Hotine积分半径为2°。计算的似大地水准面精度统计见表 2,其中包含RMS和标准差(standard deriation, SD)两项精度指标。

      表 2  Hotine积分使用梯度算子计算的重力似大地水准面精度/m

      Table 2.  Accuracy of the Gravimetric Quasi-geoid by the Hotine Integral Using Gradient Operator/m

      积分核 min max mean RMS SD
      -1.537 -0.169 -0.770 0.870 0.407
      -1.505 -0.164 -0.749 0.848 0.400
      -0.120 0.296 0.100 0.130 0.083

      表 2中可看出,两种形式的标准Hotine核函数计算的似大地水准面的精度较差。因此当使用梯度算子直接计算似大地水准面时,应使用截断形式的核函数,此时解算的似大地水准面精度为13.0 cm(均方根)、8.3 cm(标准差)。

      下面分别利用传统第三边值问题中的Molodensky理论方法和本文推导的针对第二边值问题的Molodensky算法求解似大地水准面。其中局部地形改正的内区积分半径为0.5°,采用3″地形计算,外区积分半径为1°,采用9″地形计算。表 3统计了传统第三边值的Molodensky算法与第二边值问题的Molodensky算法计算的重力似大地水准面的精度,其中Hotine积分核采用两种形式:①包含零阶、一阶项;②不含零阶、一阶项,Hotine积分半径为2°。

      表 3  Molodensky算法计算的重力似大地水准面精度/m

      Table 3.  Accuracy of the Gravimetric Quasi-geoid by Molodensky Algorithms/m

      积分核 min max mean RMS SD
      Stokes -0.217 0.108 -0.062 0.087 0.062
      Hotine① -0.251 0.131 -0.042 0.086 0.076
      Hotine② -0.210 0.080 -0.079 0.100 0.062

      表 3可以看出,第二边值问题的Molodensky算法与第三边值问题的Molodensky算法计算的重力似大地水准面达到同等精度水平,且高于使用梯度算子计算的重力似大地水准面精度,说明了Molodensky理论方法可以适用于第二边值问题。图 5所示为第二边值问题的Molodensky算法计算的重力似大地水准面。

      图  5  重力似大地水准面

      Figure 5.  The Gravimetric Quasi-geoid

    • 本文在李斐等人给出的Hotine积分实用公式的基础上,借鉴第三边值问题中计算似大地水准面的Molodensky理论方法,给出了第二边值问题的Molodensky解算方法。结论如下:

      1) 重力异常的梯度与扰动重力的梯度相差2δg/R,量级很小,因此可近似认为重力异常的梯度与扰动重力的梯度相等。

      2) Stokes算子与梯度算子间的关系式SL(f)=-γ-1I(f-f0-f1)是种近似关系;包含零阶、一阶项的Hotine算子与梯度算子严格满足HL(f)=-γ-1I(f);不包含零阶、一阶项的Hotine算子与梯度算子的关系式为HL(f)=-γ-1I(f-f0-f1)。

      3) 通过实验得出第二边值问题的Molodensky算法计算的重力似大地水准面精度与传统第三边值问题中Molodensky算法计算的重力似大地水准面精度相当,且高于直接用梯度算子计算的重力似大地水准面的精度,这说明基于Molodensky理论对第二边值问题也是适用的。

      本文借助于SRTM数据和高精度的重力场模型,证明了Molodensky算法在第二大地边值问题解算中的适用性。现阶段第二边值和第三边值已经能够以同样的精度水平计算似大地水准面,随着未来GPS对大地高测量精度的提高,第二边值必将具有更广阔的发展前景。

参考文献 (15)

目录

    /

    返回文章
    返回