留言板

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

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

基于逆Vening-Meinesz公式的测高重力中央区效应精密计算

李厚朴 边少锋 纪兵 陈永冰

李厚朴, 边少锋, 纪兵, 陈永冰. 基于逆Vening-Meinesz公式的测高重力中央区效应精密计算[J]. 武汉大学学报 ● 信息科学版, 2019, 44(2): 200-205. doi: 10.13203/j.whugis20150744
引用本文: 李厚朴, 边少锋, 纪兵, 陈永冰. 基于逆Vening-Meinesz公式的测高重力中央区效应精密计算[J]. 武汉大学学报 ● 信息科学版, 2019, 44(2): 200-205. doi: 10.13203/j.whugis20150744
LI Houpu, BIAN Shaofeng, JI Bing, CHEN Yongbing. Precise Calculation of Innermost Area Effects in Altimetry Gravity Based on the Inverse Vening-Meinesz Formula[J]. Geomatics and Information Science of Wuhan University, 2019, 44(2): 200-205. doi: 10.13203/j.whugis20150744
Citation: LI Houpu, BIAN Shaofeng, JI Bing, CHEN Yongbing. Precise Calculation of Innermost Area Effects in Altimetry Gravity Based on the Inverse Vening-Meinesz Formula[J]. Geomatics and Information Science of Wuhan University, 2019, 44(2): 200-205. doi: 10.13203/j.whugis20150744

基于逆Vening-Meinesz公式的测高重力中央区效应精密计算

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

国家自然科学基金 41631072

国家自然科学基金 41771487

国家自然科学基金 41474061

详细信息
    作者简介:

    李厚朴, 博士, 副教授, 主要从事大地测量数学分析研究。lihoupu1985@126.com

    通讯作者: 边少锋, 博士, 教授。sfbian@sina.com
  • 中图分类号: P223

Precise Calculation of Innermost Area Effects in Altimetry Gravity Based on the Inverse Vening-Meinesz Formula

Funds: 

The National Natural Science Foundation of China 41631072

The National Natural Science Foundation of China 41771487

The National Natural Science Foundation of China 41474061

More Information
    Author Bio:

    LI Houpu, PhD, associate professor, specializes in the mathematical analysis of geodesy. E-mail: lihoupu1985@126.com

    Corresponding author: BIAN Shaofeng, PhD, professor. E-mail:sfbian@sina.com
图(4) / 表(1)
计量
  • 文章访问数:  1289
  • HTML全文浏览量:  85
  • PDF下载量:  229
  • 被引次数: 0
出版历程
  • 收稿日期:  2018-06-16
  • 刊出日期:  2019-02-05

基于逆Vening-Meinesz公式的测高重力中央区效应精密计算

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

    国家自然科学基金 41631072

    国家自然科学基金 41771487

    国家自然科学基金 41474061

    作者简介:

    李厚朴, 博士, 副教授, 主要从事大地测量数学分析研究。lihoupu1985@126.com

    通讯作者: 边少锋, 博士, 教授。sfbian@sina.com
  • 中图分类号: P223

摘要: 为提高利用逆Vening-Meinesz公式反演测高重力中央区效应的精度,视中央区为矩形域,将垂线偏差分量表示成双二次多项式插值形式,引入非奇异变换,推导出了重力异常的计算公式。以低纬度区域2'×2'的垂线偏差实际数据为背景场进行了计算,结果表明,当中央区包含4个网格时,传统公式与推导出的重力异常计算公式误差的最大值大于1 mGal。推导出的公式可为高精度测高重力中央区效应的计算提供理论依据。

English Abstract

李厚朴, 边少锋, 纪兵, 陈永冰. 基于逆Vening-Meinesz公式的测高重力中央区效应精密计算[J]. 武汉大学学报 ● 信息科学版, 2019, 44(2): 200-205. doi: 10.13203/j.whugis20150744
引用本文: 李厚朴, 边少锋, 纪兵, 陈永冰. 基于逆Vening-Meinesz公式的测高重力中央区效应精密计算[J]. 武汉大学学报 ● 信息科学版, 2019, 44(2): 200-205. doi: 10.13203/j.whugis20150744
LI Houpu, BIAN Shaofeng, JI Bing, CHEN Yongbing. Precise Calculation of Innermost Area Effects in Altimetry Gravity Based on the Inverse Vening-Meinesz Formula[J]. Geomatics and Information Science of Wuhan University, 2019, 44(2): 200-205. doi: 10.13203/j.whugis20150744
Citation: LI Houpu, BIAN Shaofeng, JI Bing, CHEN Yongbing. Precise Calculation of Innermost Area Effects in Altimetry Gravity Based on the Inverse Vening-Meinesz Formula[J]. Geomatics and Information Science of Wuhan University, 2019, 44(2): 200-205. doi: 10.13203/j.whugis20150744
  • 20世纪70年代,卫星测高技术的出现以及多个卫星测高任务的顺利实施为人们获取全球海域高精度、高分辨率的海面高数据提供了极大的便利。利用卫星测高数据推求海域重力异常是卫星测高在大地测量研究中的主要应用,国内外许多学者对此进行了大量的研究,取得了丰富的研究成果[1-9]。其中,逆Vening-Meinesz公式法的输入量为垂线偏差,而垂线偏差是由测高观测值的一次差分求得,可以消除与地理位置相关的径向轨道误差和长波海面地形等类似系统误差,同时含有丰富的重力场高频成分,非常有利于高分辨率海洋重力场的恢复。因此,逆Vening-Meinesz公式在测高重力反演中的应用较为广泛[10]

    在计算点本身及其附近区域计算重力异常时,积分区域包含计算点,导致逆Vening-Meinesz公式中的积分奇异,本文称该积分区域为中央区。Hwang[5]将中央区视为圆域,推导出了中央区重力异常的计算公式;常晓涛等[9]将中央区视为方形域,推导出了中央区效应的计算公式,指出中央区效应与奇异区内垂线偏差分量的梯度及奇异区面积的大小有关。由于实际数据通常为网格化分布,受子午线收敛的影响,中央区更接近于矩形域,因此,将中央区视为圆域和方形域的处理方法与数据真实分布情况并不相符,由此产生的误差在高精度重力反演中能否忽略值得深入研究。鉴于此,本文推导出了视中央区为矩形域时的重力异常计算公式,并对导出公式和圆域与方形域下传统公式的误差进行了分析比较。

    • 略去详细的推导过程,根据文献[5],直接写出逆Vening-Meinesz公式如下:

      $$ \Delta g = \frac{{{\gamma _0}}}{{4{\pi }}}\iint\limits_\sigma {H'\left( {{\psi _{PQ}}} \right)\left( {{\xi _Q}\cos {\alpha _{QP}} + {\eta _Q}\sin {\alpha _{QP}}} \right){\text{d}}\sigma } $$ (1)

      式中,γ0=979.8Gal为地球平均重力;αQP代表流动点Q至计算点P的方位角;ψPQPQ之间的球面距离;ξQηQ分别代表流动点Q处的垂线偏差在子午圈方向和卯酉圈方向上的分量;H′(ψPQ)为积分核函数,其表达式为:

      $$ H'\left( {{\psi _{PQ}}} \right) = - \frac{{\cos \frac{{{\psi _{PQ}}}}{2}}}{{2{{\sin }^2}\frac{{{\psi _{PQ}}}}{2}}} + \frac{{\cos \frac{{{\psi _{PQ}}}}{2}\left( {3 + 2\sin \frac{{{\psi _{PQ}}}}{2}} \right)}}{{2\sin \frac{{{\psi _{PQ}}}}{2}\left( {1 + \sin \frac{{{\psi _{PQ}}}}{2}} \right)}} $$ (2)

      为对中央区进行计算,首先定义以计算点P为原点的局部切平面直角坐标系,如图 1所示,x轴指向北极,y轴指向东,且xy平面在P点与地球表面相切,Q(x, y)为中央区内流动的积分点,l为计算点与流动点的平面距离。由图 1可知:

      图  1  局部坐标系示意图

      Figure 1.  The Sketch Map of the Local Coordinate System

      $$ \left\{ \begin{array}{l} \cos {\alpha _{QP}} = - \frac{x}{l}\\ \sin {\alpha _{QP}} = - \frac{y}{l} \end{array} \right. $$ (3)

      式中,$l = \sqrt {{x^2} + {y^2}} $。

      当积分区域比较小时,球面角距ψPQ可以简化为:

      $$ {\psi _{PQ}} \approx \frac{l}{R} $$ (4)

      则式(1)中的积分核函数可以简化为:

      $$ \begin{array}{*{20}{c}} {H'\left( {{\psi _{PQ}}} \right) \approx - \frac{{\cos \frac{{{\psi _{PQ}}}}{2}}}{{2{{\sin }^2}\frac{{{\psi _{PQ}}}}{2}}} + \frac{{3\cos \frac{{{\psi _{PQ}}}}{2}}}{{2\sin \frac{{{\psi _{PQ}}}}{2}}} \approx - }\\ {\frac{1}{{2{{\left( {\frac{{{\psi _{PQ}}}}{2}} \right)}^2}}} = - \frac{2}{{\psi _{PQ}^2}} \approx - \frac{{2{R^2}}}{{{l^2}}}} \end{array} $$ (5)

      将式(3)和式(5)代入式(1),顾及平面近似下积分面元为R2dσ=dxdy,可得中央区重力异常的计算公式为:

      $$ \begin{array}{*{20}{c}} {\Delta g = \frac{{{\gamma _0}}}{{4{\pi}}}\iint\limits_\sigma {\left( { - \frac{{2{R^2}}}{{{l^2}}}} \right)\left( {\left( { - \frac{x}{l}} \right){\xi _Q} + \left( { - \frac{y}{l}} \right){\eta _Q}} \right)} \cdot } \\ {\frac{1}{{{R^2}}}{\text{d}}x{\text{d}}y + \frac{{{\gamma _0}}}{{2{\pi}}}\iint\limits_\sigma {\frac{{{\xi _Q}x + {\eta _Q}y}}{{{{\left( {{x^2} + {y^2}} \right)}^{3/2}}}}{\text{d}}x{\text{d}}y}} \end{array} $$ (6)

      可以看出,式(6)中的积分在计算点P处奇异。为解决这一问题,文献[5]将中央区视为圆域,将中央区内的垂线偏差分量ξQηQ展开为泰勒级数形式:

      $$ \begin{array}{*{20}{c}} {\xi \left( {x,y} \right) = {\xi _P} + {\xi _x}x + {\xi _y}y + \frac{1}{{2!}}\left( {{\xi _{xx}}{x^2} + } \right.}\\ {\left. {2{\xi _{xy}}xy + {\xi _{yy}}{y^2}} \right) + \cdots } \end{array} $$ (7)
      $$ \begin{array}{*{20}{c}} {\eta \left( {x,y} \right) = {\eta _P} + {\eta _x}x + {\eta _y}y + \frac{1}{{2!}}\left( {{\eta _{xx}}{x^2} + } \right.}\\ {\left. {2{\eta _{xy}}xy + {\eta _{yy}}{y^2}} \right) + \cdots } \end{array} $$ (8)

      式中,${\xi _x} = \frac{{\partial \xi }}{{\partial x}}, {\xi _y} = \frac{{\partial \xi }}{{\partial y}}$,${\xi _{xx}} = \frac{{{\partial ^2}\xi }}{{\partial {x^2}}}, {\xi _{yy}} = \frac{{{\partial ^2}\xi }}{{\partial {y^2}}}$,${\xi _{xy}} = \frac{{{\partial ^2}\xi }}{{\partial x\partial y}}, {\eta _x} = \frac{{\partial \eta }}{{\partial x}}$,${\eta _y} = \frac{{\partial \eta }}{{\partial y}}, {\eta _{xx}} = \frac{{{\partial ^2}\eta }}{{\partial {x^2}}}$,${\eta _{yy}} = \frac{{{\partial ^2}\eta }}{{\partial {y^2}}}, {\eta _{xy}} = \frac{{{\partial ^2}\eta }}{{\partial x\partial y}}$。

      导出的中央区重力异常计算公式为:

      $$ \Delta {g_1} = \frac{{{\gamma _0}}}{2}\sqrt {\frac{{\Delta x\Delta y}}{{\rm{ \mathsf{ π} }}}} \left( {{\xi _x} + {\eta _y}} \right) $$ (9)

      式中,Δx、Δy分别表示x方向和y方向的网格间距。

      文献[9]将中央区视为方形域,将垂线偏差分量同样表示为式(7)和式(8),导出的重力异常计算公式为:

      $$ \Delta {g_2} = \frac{{2\ln \left( {1 + \sqrt 2 } \right)}}{{\rm{ \mathsf{ π} }}}s{\gamma _0}\left( {{\xi _x} + {\eta _y}} \right) $$ (10)

      式中,s为方形域长度的一半。

    • 首先将中央区垂线偏差分量表示为双二次多项式插值形式,之后利用文献[11-13]提出的非奇异变换对逆Vening-Meinesz公式中的奇异积分进行了处理,推导出了中央区重力异常的精密计算公式。

    • 图 2所示,设中央区为σ[-a < x < a, -b < y < b],其中a=1,b=cosφ。中央区共包含4个网格单元,9个网格节点。将垂线偏差子午分量ξQ、卯酉分量ηQ表示成双二次多项式插值形式:

      图  2  垂线偏差双二次多项式插值表示时的中央区

      Figure 2.  The Innermost Area when the Components of Deflections of the Vertical are Expressed as Bi-quad- ratic Polynomials

      $$ \mathit{\boldsymbol{\xi }}\left( {x,y} \right) = \sum\limits_{i = 0}^2 {\sum\limits_{j = 0}^2 {{\alpha _{ij}}{x^i}{y^j}} } = \sum\limits_{i = 0}^2 {{x^i}\sum\limits_{j = 0}^2 {{\alpha _{ij}}{y^j}} } $$ (11)
      $$ \eta \left( {x,y} \right) = \sum\limits_{i = 0}^2 {\sum\limits_{j = 0}^2 {{\beta _{ij}}{x^i}{y^j}} } = \sum\limits_{i = 0}^2 {{x^i}\sum\limits_{j = 0}^2 {{\beta _{ij}}{y^j}} } $$ (12)

      为确定出待定系数αijβij,将式(11)和式(12)改写为:

      $$ \begin{array}{*{20}{c}} {\xi \left( {x,y} \right) = \sum\limits_{i = 0}^2 {{x^i}\sum\limits_{j = 0}^2 {{c_{ij}}{{\left( {y/\cos \varphi } \right)}^j}} } = }\\ {\left( {1,x,{x^2}} \right)\left( {\begin{array}{*{20}{c}} {{c_{00}}}&{{c_{01}}}&{{c_{02}}}\\ {{c_{10}}}&{{c_{11}}}&{{c_{12}}}\\ {{c_{20}}}&{{c_{21}}}&{{c_{22}}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} 1\\ {y/\cos \varphi }\\ {{y^2}/{{\cos }^2}\varphi } \end{array}} \right)} \end{array} $$ (13)
      $$ \begin{array}{*{20}{c}} {\eta \left( {x,y} \right) = \sum\limits_{i = 0}^2 {{x^i}\sum\limits_{j = 0}^2 {{d_{ij}}{{\left( {y/\cos \varphi } \right)}^j}} } = }\\ {\left( {1,x,{x^2}} \right)\left( {\begin{array}{*{20}{c}} {{d_{00}}}&{{d_{01}}}&{{d_{02}}}\\ {{d_{10}}}&{{d_{11}}}&{{d_{12}}}\\ {{d_{20}}}&{{d_{21}}}&{{d_{22}}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} 1\\ {y/\cos \varphi }\\ {{y^2}/{{\cos }^2}\varphi } \end{array}} \right)} \end{array} $$ (14)

      则有:

      $$ {\alpha _{ij}} = {c_{ij}}{\left( {1/\cos \varphi } \right)^j} $$ (15)
      $$ {\beta _{ij}} = {d_{ij}}{\left( {1/cos\varphi } \right)^j} $$ (16)

      图 2中网格节点处的垂线偏差子午分量ξij=ξ(i, jcosφ)和卯酉分量ηij=η(i, jcosφ)(i, j=-1, 0, 1)作为插值条件代入式(13)和式(14),可得:

      $$ \begin{array}{*{20}{c}} {\left( {\begin{array}{*{20}{c}} 1&{ - 1}&1\\ 1&0&0\\ 1&1&1 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{c_{00}}}&{{c_{01}}}&{{c_{02}}}\\ {{c_{10}}}&{{c_{11}}}&{{c_{12}}}\\ {{c_{20}}}&{{c_{21}}}&{{c_{22}}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} 1&1&1\\ { - 1}&0&1\\ 1&0&1 \end{array}} \right) = }\\ {\left( {\begin{array}{*{20}{c}} {{\xi _{ - 1 - 1}}}&{{\xi _{ - 10}}}&{{\xi _{ - 11}}}\\ {{\xi _{0 - 1}}}&{{\xi _{00}}}&{{\xi _{01}}}\\ {{\xi _{1 - 1}}}&{{\xi _{10}}}&{{\xi _{11}}} \end{array}} \right)} \end{array} $$ (17)
      $$ \begin{array}{*{20}{c}} {\left( {\begin{array}{*{20}{c}} 1&{ - 1}&1\\ 1&0&0\\ 1&1&1 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{d_{00}}}&{{d_{01}}}&{{d_{02}}}\\ {{d_{10}}}&{{d_{11}}}&{{d_{12}}}\\ {{d_{20}}}&{{d_{21}}}&{{d_{22}}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} 1&1&1\\ { - 1}&0&1\\ 1&0&1 \end{array}} \right) = }\\ {\left( {\begin{array}{*{20}{c}} {{\eta _{ - 1 - 1}}}&{{\eta _{ - 10}}}&{{\eta _{ - 11}}}\\ {{\eta _{0 - 1}}}&{{\eta _{00}}}&{{\eta _{01}}}\\ {{\eta _{1 - 1}}}&{{\eta _{10}}}&{{\eta _{11}}} \end{array}} \right)} \end{array} $$ (18)

      式(17)、(18)可以简化为:

      $$ \mathit{\boldsymbol{A}}\left( {{c_{ij}}} \right){\mathit{\boldsymbol{A}}^{\rm{T}}} = \left( {{\xi _{ij}}} \right) $$ (19)
      $$ \mathit{\boldsymbol{A}}\left( {{d_{ij}}} \right){\mathit{\boldsymbol{A}}^{\rm{T}}} = \left( {{\eta _{ij}}} \right) $$ (20)

      其中,

      $$ \mathit{\boldsymbol{A}} = \left( {\begin{array}{*{20}{c}} 1&{ - 1}&1\\ 1&0&0\\ 1&1&1 \end{array}} \right) $$ (21)

      系数矩阵A非奇异、可逆,并且:

      $$ {\mathit{\boldsymbol{A}}^{ - 1}} = \left( {\begin{array}{*{20}{c}} 0&1&0\\ { - \frac{1}{2}}&0&{\frac{1}{2}}\\ {\frac{1}{2}}&{ - 1}&{\frac{1}{2}} \end{array}} \right) $$ (22)

      因此有:

      $$ \left( {{c_{ij}}} \right) = {\mathit{\boldsymbol{A}}^{ - 1}}\left( {{\xi _{ij}}} \right){\left( {{\mathit{\boldsymbol{A}}^{ - 1}}} \right)^{\rm{T}}} $$ (23)
      $$ \left( {{d_{ij}}} \right) = {\mathit{\boldsymbol{A}}^{ - 1}}\left( {{\eta _{ij}}} \right){\left( {{\mathit{\boldsymbol{A}}^{ - 1}}} \right)^{\rm{T}}} $$ (24)

      将式(23)和式(24)代入式(15)和式(16)即可确定出待定系数αijβij

    • 为推导方便,记式(6)中的奇异积分为:

      $$ IV = \iint\limits_\sigma {\frac{{{\xi _Q}x + {\eta _Q}y}}{{{{\left( {{x^2} + {y^2}} \right)}^{3/2}}}}{\text{d}}x{\text{d}}y} = I{V_\xi } + I{V_\eta } $$ (25)

      其中,

      $$ I{V_\xi } = \iint\limits_\sigma {\frac{{{\xi _Q}x}}{{{{\left( {{x^2} + {y^2}} \right)}^{3/2}}}}{\text{d}}x{\text{d}}y} $$ (26)
      $$ I{V_\eta } = \iint\limits_\sigma {\frac{{{\eta _Q}y}}{{{{\left( {{x^2} + {y^2}} \right)}^{3/2}}}}{\text{d}}x{\text{d}}y} $$ (27)

      将式(11)和式(12)分别代入式(26)和式(27),并考虑到奇偶函数的积分性质,则有:

      $$ I{V_\xi } = \iint\limits_\sigma {\frac{{\left( {{\alpha _{10}} + {\alpha _{12}}{y^2}} \right){x^2}}}{{{{\left( {{x^2} + {y^2}} \right)}^{3/2}}}}{\text{d}}x{\text{d}}y} $$ (28)
      $$ I{V_\eta } = \iint\limits_\sigma {\frac{{\left( {{\beta _{01}} + {\beta _{21}}{x^2}} \right){y^2}}}{{{{\left( {{x^2} + {y^2}} \right)}^{3/2}}}}{\text{d}}x{\text{d}}y} $$ (29)

      由计算点P向右上顶点作连线分右上象限为σ1[0 < x < 1, 0 < y < bx]和σ2[0 < x < b-1y, 0 < y < b](见图 2)。

      σ1引入非奇异变换[11-13]

      $$ \left\{ \begin{array}{l} x = x\\ y = kx \end{array} \right. $$ (30)

      σ2引入非奇异变换[11-13]

      $$ \left\{ \begin{array}{l} x = \lambda y\\ y = y \end{array} \right. $$ (31)

      则三角形σ1映射为矩形σ1[0 < x < 1, 0 < k < b],三角形σ2映射为矩形σ2[0 < λ < b-1, 0 < y < b]。将式(30)和式(31)分别代入式(28)和式(29),并注意到两式中的被积函数均为偶函数,因此有:

      $$ \begin{gathered} I{V_\xi } = 4\iint\limits_{{\sigma _1}} {\frac{{\left( {{\alpha _{10}} + {\alpha _{12}}{y^2}} \right){x^2}}}{{{{\left( {{x^2} + {y^2}} \right)}^{3/2}}}}{\text{d}}x{\text{d}}y} + 4\iint\limits_{{\sigma _2}} {\frac{{\left( {{\alpha _{10}} + {\alpha _{12}}{y^2}} \right){x^2}}}{{{{\left( {{x^2} + {y^2}} \right)}^{3/2}}}}{\text{d}}x{\text{d}}y} =\\ 4\int_0^b {\int_0^1 {\frac{{{\alpha _{10}} + {\alpha _{12}}{k^2}{x^2}}}{{{{\left( {1 + {k^2}} \right)}^{3/2}}}}{\text{d}}x{\text{d}}k} } + \hfill \\ 4\int_0^{{b^{ - 1}}} {\int_0^b {\frac{{\left( {{\alpha _{10}} + {\alpha _{12}}{y^2}} \right){\lambda ^2}}}{{{{\left( {1 + {\lambda ^2}} \right)}^{3/2}}}}{\text{d}}y{\text{d}}\lambda } } = 4\int_0^b {\frac{{{\alpha _{10}} + \frac{1}{3}{\alpha _{12}}{k^2}}}{{{{\left( {1 + {k^2}} \right)}^{3/2}}}}{\text{d}}k} +\\ 4b\left( {{\alpha _{10}} + \frac{1}{3}{\alpha _{12}}{b^2}} \right)\int_0^{{b^{ - 1}}} {\frac{{{\lambda ^2}}}{{{{\left( {1 + {\lambda ^2}} \right)}^{3/2}}}}{\text{d}}\lambda } \hfill \\ \end{gathered} $$ (32)
      $$ \begin{gathered} I{V_\eta } = 4\iint\limits_{{\sigma _1}} {\frac{{\left( {{\beta _{01}} + {\beta _{21}}{x^2}} \right){y^2}}}{{{{\left( {{x^2} + {y^2}} \right)}^{3/2}}}}{\text{d}}x{\text{d}}y} + 4\iint\limits_{{\sigma _2}} {\frac{{\left( {{\beta _{01}} + {\beta _{21}}{x^2}} \right){y^2}}}{{{{\left( {{x^2} + {y^2}} \right)}^{3/2}}}}{\text{d}}x{\text{d}}y} =\\ 4\int_0^b {\int_0^1 {\frac{{\left( {{\beta _{01}} + {\beta _{21}}{x^2}} \right){k^2}}}{{{{\left( {1 + {k^2}} \right)}^{3/2}}}}{\text{d}}x{\text{d}}k} } + \hfill \\ 4\int_0^{{b^{ - 1}}} {\int_0^b {\frac{{{\beta _{01}} + {\beta _{21}}{\lambda ^2}{y^2}}}{{{{\left( {1 + {\lambda ^2}} \right)}^{3/2}}}}{\text{d}}y{\text{d}}\lambda } } = 4\left( {{\beta _{01}} + \frac{1}{3}{\beta _{21}}} \right)\int_0^b {\frac{{{k^2}}}{{{{\left( {1 + {k^2}} \right)}^{3/2}}}}{\text{d}}k} +\\ 4b\int_0^{{b^{ - 1}}} {\frac{{{\beta _{01}} + \frac{1}{3}{\beta _{21}}{\lambda ^2}{b^2}}}{{{{\left( {1 + {\lambda ^2}} \right)}^{3/2}}}}{\text{d}}\lambda } \hfill \\ \end{gathered} $$ (33)

      可以看出,在非奇异变换下,原来含有xy两变量的二维积分转换为只含k(λ)变量的一维积分,利用Mathematica计算机代数系统[14-15]可求得:

      $$ \begin{array}{l} I{V_\xi } = 4\left( {\frac{b}{{\sqrt {1 + {b^2}} }}{\alpha _{10}} + \frac{1}{3}{\alpha _{12}}\left( {{\rm{arsh}}\left( b \right) - \frac{b}{{\sqrt {1 + {b^2}} }}} \right) + } \right.\\ \;\;\;\left. {b\left( {{\alpha _{10}} + \frac{1}{3}{\alpha _{12}}{b^2}} \right)\left( {{\rm{arsh}}\left( {{b^{ - 1}}} \right) - \frac{1}{{\sqrt {1 + {b^2}} }}} \right)} \right) \end{array} $$ (34)
      $$ \begin{array}{l} I{V_\eta } = 4\left( {\frac{b}{{\sqrt {1 + {b^2}} }}{\beta _{01}} + \left( {{\beta _{01}} + \frac{1}{3}{\beta _{21}}} \right)\left( {{\rm{arsh}}\left( b \right) - } \right.} \right.\\ \;\;\;\left. {\left. {\frac{b}{{\sqrt {1 + {b^2}} }}} \right) + \frac{1}{3}{\beta _{21}}{b^3}\left( {{\rm{arsh}}\left( {{b^{ - 1}}} \right) - \frac{1}{{\sqrt {1 + {b^2}} }}} \right)} \right) \end{array} $$ (35)
      $$ {\rm{arsh}}\left( b \right) = \ln \left( {b + \sqrt {1 + {b^2}} } \right) $$ (36)
      $$ {\rm{arsh}}\left( {{b^{ - 1}}} \right) = \ln \left( {{b^{ - 1}} + \sqrt {1 + {b^{ - 2}}} } \right) $$ (37)

      因此,式(6)中的奇异积分可以表示为:

      $$ \begin{array}{*{20}{c}} {IV = 4\left( {\frac{b}{{\sqrt {1 + {b^2}} }}\left( {{\alpha _{10}} + {\beta _{01}}} \right) + \left( {\frac{1}{3}{\alpha _{12}} + {\beta _{01}} + \frac{1}{3}{\beta _{21}}} \right) \cdot } \right.}\\ {\left( {{\rm{arsh}}\left( b \right) - \frac{b}{{\sqrt {1 + {b^2}} }}} \right) + b\left( {{\alpha _{10}} + \frac{1}{3}{\alpha _{12}}{b^2} + \frac{1}{3}{\beta _{21}}{b^2}} \right) \cdot }\\ {\left. {\left( {{\rm{arsh}}\left( {{b^{ - 1}}} \right) - \frac{1}{{\sqrt {1 + {b^2}} }}} \right)} \right)} \end{array} $$ (38)

      将式(38)代入式(6),可得中央区包含4个网格时,垂线偏差双二次多项式插值表示下重力异常的计算公式为:

      $$ \begin{array}{l} \Delta {g_3} = \frac{{2{\gamma _0}}}{{\rm{ \mathsf{ π} }}}\left( {\frac{b}{{\sqrt {1 + {b^2}} }}\left( {{\alpha _{10}} + {\beta _{01}}} \right) + \left( {\frac{1}{3}{\alpha _{12}} + {\beta _{01}} + } \right.} \right.\\ \left. {\frac{1}{3}{\beta _{21}}} \right)\left( {{\rm{arsh}}\left( b \right) - \frac{b}{{\sqrt {1 + {b^2}} }}} \right) + b\left( {{\alpha _{10}} + \frac{1}{3}{\alpha _{12}}{b^2} + } \right.\\ \left. {\left. {\frac{1}{3}{\beta _{21}}{b^2}} \right)\left( {{\rm{arsh}}\left( {{b^{ - 1}}} \right) - \frac{b}{{\sqrt {1 + {b^2}} }}} \right)} \right) \end{array} $$ (39)

      若将积分区域取为$\sigma '\left[{-\frac{1}{2} < x < \frac{1}{2}, -\frac{1}{2}b < y < \frac{1}{2}b} \right]$,则计算点P所在的1个网格对重力异常的贡献可采用类似方法得到:

      $$ \begin{array}{l} \Delta {g_4} = \frac{{{\gamma _0}}}{{2\rm{ \mathsf{ π} }}}\left( {\frac{{2b}}{{\sqrt {1 + {b^2}} }}\left( {{\alpha _{10}} + {\beta _{01}}} \right) + \left( {\frac{1}{6}{\alpha _{12}} + 2{\beta _{01}} + } \right.} \right.\\ \left. {\frac{1}{6}{\beta _{21}}} \right)\left( {{\rm{arsh}}\left( b \right) - \frac{b}{{\sqrt {1 + {b^2}} }}} \right) + b\left( {2{\alpha _{10}} + \frac{1}{6}{\alpha _{12}}{b^2} + } \right.\\ \;\;\;\;\;\;\;\;\;\;\left. {\left. {\frac{1}{6}{\beta _{21}}{b^2}} \right)\left( {{\rm{arsh}}\left( {{b^{ - 1}}} \right) - \frac{1}{{\sqrt {1 + {b^2}} }}} \right)} \right) \end{array} $$ (40)
    • 选定中国南海12°N~18°N,112°E~118°E海域作为试算区,以EGM2008地球重力场模型计算得到的2′×2′分辨率的垂线偏差数据作为背景场,计算了该区域的中央区效应。

      试算区共179×179个网格,假定中央区为4′×4′,即包含4个网格,利用式(9)、式(10)、式(39)计算了该区域的中央区效应,结果分别为Δg1、Δg2、Δg3, 计算结果之间的比较情况如表 1所示。Δg1-Δg3、Δg2-Δg3计算结果分别如图 3图 4所示。

      表 1  3种方法计算结果之间的比较/mGal

      Table 1.  Comparisons of the Innermost Effects Computed by Three Methods/mGal

      差值 最小值 最大值 标准差
      Δg1-Δg3 -1.229 1.276 0.142
      Δg2-Δg3 -1.192 1.156 0.129
      Δg2-Δg1 -0.199 0.128 0.024

      图  3  Δg1-Δg3示意图

      Figure 3.  The Sketch Map of Δg1-Δg3

      图  4  Δg2-Δg3示意图

      Figure 4.  The Sketch Map of Δg2-Δg3

      表 1图 3图 4可以看出,当中央区包含4个网格时,文献[5]给出的式(9)与本文导出的式(39)计算得到的中央区重力异常差值的标准差均为0.142 mGal,最大值达1.276 mGal;文献[9]给出的式(10)与本文导出的式(39)计算得到的中央区重力异常差值的标准差均为0.129 mGal,最大值达1.156 mGal。

    • 为提高利用逆Vening-Meinesz公式反演中央区重力异常的精度,本文将中央区垂线偏差表示成双二次多项式插值形式,推导出了与实际数据分布更为相符的矩形域中央区效应精密计算公式,并利用模型数据计算分析了导出公式和传统公式的误差。研究表明:

      1) 非奇异变换可化逆Vening-Meinesz公式中的奇异积分为非奇异积分,可在其他地球重力场奇异积分问题中推广使用。

      2) 低纬度试算区2′×2′分辨率垂线偏差数据下的中央区效应分析表明,当中央区包含4个网格时,传统公式与本文导出公式计算结果差值的均方差和标准差大于0.1 mGal,最大值大于1 mGal,需要在高精度测高重力计算中加以考虑。

参考文献 (15)

目录

    /

    返回文章
    返回