留言板

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

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

航空重力向下延拓的多参数正则化方法

徐新强 赵俊

徐新强, 赵俊. 航空重力向下延拓的多参数正则化方法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(7): 956-963, 973. doi: 10.13203/j.whugis20180335
引用本文: 徐新强, 赵俊. 航空重力向下延拓的多参数正则化方法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(7): 956-963, 973. doi: 10.13203/j.whugis20180335
XU Xinqiang, ZHAO Jun. A Multi-Parameter Regularization Method in Downward Continuation for Airborne Gravity Data[J]. Geomatics and Information Science of Wuhan University, 2020, 45(7): 956-963, 973. doi: 10.13203/j.whugis20180335
Citation: XU Xinqiang, ZHAO Jun. A Multi-Parameter Regularization Method in Downward Continuation for Airborne Gravity Data[J]. Geomatics and Information Science of Wuhan University, 2020, 45(7): 956-963, 973. doi: 10.13203/j.whugis20180335

航空重力向下延拓的多参数正则化方法

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

大地测量与地球动力学国家重点实验室开放基 SKLGED 2017-3-2-E

详细信息
    作者简介:

    徐新强,高级工程师,主要从事地球重力场方面的数据处理工作及研究。xaxcq@126.com

    通讯作者: 赵俊,博士,工程师。Zhaojun4368@163.com
  • 中图分类号: P223

A Multi-Parameter Regularization Method in Downward Continuation for Airborne Gravity Data

Funds: 

State Key Laboratory of Geodesy and Earth's Dynamics SKLGED 2017-3-2-E

More Information
    Author Bio:

    XU Xinqiang, senior engineer, specializes in the data processing of earth gravity survey. E-mail: xaxcq@126.com

图(5) / 表(2)
计量
  • 文章访问数:  697
  • HTML全文浏览量:  145
  • PDF下载量:  42
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-04-07
  • 刊出日期:  2020-07-30

航空重力向下延拓的多参数正则化方法

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

    大地测量与地球动力学国家重点实验室开放基 SKLGED 2017-3-2-E

    作者简介:

    徐新强,高级工程师,主要从事地球重力场方面的数据处理工作及研究。xaxcq@126.com

    通讯作者: 赵俊,博士,工程师。Zhaojun4368@163.com
  • 中图分类号: P223

摘要: 为了克服航空重力向下延拓解算的病态性影响,介绍了一种多参数正则化方法,以均方误差最小为目标函数,设计了选取正则化参数的迭代算法,并比较了基于L曲线法、广义交叉核实(generalized cross-validation,GCV)方法选取正则化参数的Tikhonov正则化方法,同时给出了均方误差意义下多参数正则化解优于最小二乘估计的条件。基于EGM2008地球重力场模型进行了仿真试验,计算结果表明,多参数正则化方法能够保证向下延拓结果的可靠性和稳定性,并优于现有的Tikhonov正则化方法,验证了多参数方法在航空重力向下延拓中的可行性。

English Abstract

徐新强, 赵俊. 航空重力向下延拓的多参数正则化方法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(7): 956-963, 973. doi: 10.13203/j.whugis20180335
引用本文: 徐新强, 赵俊. 航空重力向下延拓的多参数正则化方法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(7): 956-963, 973. doi: 10.13203/j.whugis20180335
XU Xinqiang, ZHAO Jun. A Multi-Parameter Regularization Method in Downward Continuation for Airborne Gravity Data[J]. Geomatics and Information Science of Wuhan University, 2020, 45(7): 956-963, 973. doi: 10.13203/j.whugis20180335
Citation: XU Xinqiang, ZHAO Jun. A Multi-Parameter Regularization Method in Downward Continuation for Airborne Gravity Data[J]. Geomatics and Information Science of Wuhan University, 2020, 45(7): 956-963, 973. doi: 10.13203/j.whugis20180335
  • 航空重力测量是获取地球重力场观测数据的重要方式,而要将其应用到大地测量等相关领域,则需将航空重力数据进行格网化并归算到大地水准面上,故而要采用向下延拓技术。目前,关于航空重力数据向下延拓方法的讨论已成为国内外学者关注的焦点和热点[1-4]

    航空重力向下延拓方法主要有直接代表法、虚拟点质量法、正则化方法以及利用球内Dirich‐ let问题解向下延拓等。陈欣等[5]对这些方法进行了比较分析,结合大同地区实测航空重力数据,应用密集地面实测重力数据进行检核,得出基于Poisson积分的正则化方法精度最高的结论。但Poisson积分逆问题的解算方程存在严重的病态性,如果不采取有效措施对其进行处理,其解算结果会严重偏离参考值,进而得到不可靠的结论。为此,许多学者将正则化方法引入航空重力测量数据向下延拓的数据处理中,并得到了一些有意义的结论,验证了该方法的可行性[6-7]。为了进一步提高航空重力数据向下延拓的精度,在经典的正则化方法的基础上,许多学者进一步讨论了相对应的正则化解算策略,例如顾勇为等[8-9]研究了多种Tikhonov正则化策略,通过信噪比方法诊断系数矩阵的复共线性,有针对性地构造正则化矩阵,有效地改进了先前正则化方法的不足;邓凯亮等[10]、蒋涛等[11]分别应用双参数正则化方法和广义岭估计方法在航空重力向下延拓中的应用进行了探讨,均取得了良好应用效果。然而,以上研究大多数均以单个正则化参数为基础。近年来,相关学者讨论了多个正则化参数,并给出了相关参数选取方法和解算策略。但现有方法存在正则化参数选取较为复杂、计算困难的不足[12]。为此,本文介绍了一种简洁和方便计算的多参数正则化方法,结合MSE(mean square of estimation)标准选取正则化参数,并将其应用于航空重力数据向下延拓,验证了该方法的可行性。

    • 根据位理论和Dirichlet边值理论,重力异常可通过Poisson积分公式计算,表达式为:

      $$ \Delta {{g}^{h}}(r, \phi , \lambda )=\frac{{{R}^{2}}({{r}^{2}}-{{R}^{2}})}{4\pi r}\iint\limits_{\sigma }{\frac{\Delta g(R, {{\phi }^{\prime }}, {{\lambda }^{\prime }})}{{{\rho }^{3}}}}\text{d}\sigma $$ (1)

      式中,R为地球半径;r = R + H,为计算点的地心距离,H为延拓高度;ϕλ为计算点的纬度和经度;ϕ'、λ'为单位流动面元dσ的纬度、经度;ρ为计算点与流动面元间的距离。其计算公式为:

      $$ {\rho = \sqrt {{r^2} + {R^2} - 2rR{\kern 1pt} {\kern 1pt} {\rm{cos}}\psi } } $$ (2)
      $$ {{\rm{cos}}\psi = {\rm{sin}}{\phi ^\prime }{\kern 1pt} {\rm{sin}}\phi + {\rm{cos}}{\phi ^\prime }{\kern 1pt} {\rm{cos}}\phi {\kern 1pt} {\rm{cos}}({\lambda ^\prime } - \lambda )} $$ (3)

      鉴于航空重力异常和地面重力异常均为离散值,于是将式(1)离散化,有:

      $$ \begin{array}{*{20}{c}} {\Delta {g^h}(r, {\varphi _i}, {\lambda _i}) - \Delta {g_{\Pi - \varOmega }}(r, {\varphi _j}, {\lambda _j}) = }\\ {\sum\limits_{i = 1}^t {{A_{ij}}} \Delta g(R, \varphi _j^\prime , \lambda _j^\prime )} \end{array} $$ (4)

      式中,Ω为有效积分范围;Δg∏ - Ω(rφ jλ j)为远区影响,可通过重力场模型近似计算得到。为了抑制截断误差的影响,保证向下延拓精度的可靠性,这里采用移去‐恢复技术进行处理,具体计算步骤见文献[10]。进而可将式(4)写成矩阵形式:

      $$ {\mathit{\boldsymbol{L}}_{\Delta {g^h}}} = \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{\Delta g}} $$ (5)

      其中,LΔg h = Δg h(r) - Δg∏ - Ω (r);A=[A ij]n × tXΔg = Δg (R)。

      利用最小二乘(least-square,LS)方法,可得到地面重力估计值为:

      $$ {\mathit{\boldsymbol{\hat X}}_{\Delta g}} = {({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}})^{ - 1}}\mathit{\boldsymbol{A}}{\mathit{\boldsymbol{L}}_{\Delta {g^{\rm{h}}}}} $$ (6)

      上述公式中AT A的条件数较大,呈现出严重的病态性,微小扰动就会对解造成巨大影响,因而需采用正则化方法进行求解。许多学者就相关方法进行了深入研究,主要集中于单个正则化参数的Tikhonov正则化方法。为了进一步保证向下延拓的可靠性和精度,本文介绍了一种多参数正则化方法。

    • 考虑如下Gauss-Markov模型[13]

      $$ \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{AX}} + \mathit{\boldsymbol{ \boldsymbol{\varDelta} \boldsymbol{\varDelta} }}\backsim \mathit{\boldsymbol{N}}(0, \sigma _0^2{\mathit{\boldsymbol{P}}^{ - 1}}) $$ (7)

      式中,Ln × 1观测向量;An × t系数矩阵;Xt × 1未知参数向量;Δn × 1随机误差向量。通常仅考虑权阵为P = I,则有:

      $$ \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{AX}} + \varepsilon $$ (8)

      式中,εn × 1随机误差向量,ε ~ N (0,σ02 In)。

      模型(7)的LS估计为:

      $$ {\mathit{\boldsymbol{\hat X}}_{{\rm{LS}}}} = {({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}})^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{L}} $$ (9)

      如果对系数矩阵作奇异值分解(singular value decomposition,SVD),有:

      $$ \mathit{\boldsymbol{A}} = \mathit{\boldsymbol{US}}{\mathit{\boldsymbol{V}}^{\rm{T}}} $$ (10)

      式中,Un× t的正交矩阵;S为由系数矩阵A的特征值组成的对角矩阵,即S=diag(λ1λ2λ t);Vt × t的正交矩阵。

      对式(9)变形,得:

      $$ \mathit{\boldsymbol{\hat X}} = {\mathit{\boldsymbol{A}}^ + }\mathit{\boldsymbol{L}} $$ (11)

      其中,A+ = VS-1U T为影响矩阵。显然式(11)跟式(9)等价。

      当系数矩阵A的条件数非常大时,式(11)呈现出严重的病态性,得到的未知参数估值$\mathit{\boldsymbol{\hat X}}$极不稳定。采用Tikhonov正则化方法,即最小化如下目标函数[14-15]

      $$ \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} = {\left\| {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{AX}}} \right\|^2} + \alpha {\left\| \mathit{\boldsymbol{X}} \right\|^2} $$ (12)

      可得未知参数估计为:

      $$ \mathit{\boldsymbol{\hat X}} = {({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}} + \alpha {\mathit{\boldsymbol{I}}_t})^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{L}} $$ (13)

      上述正则化参数可通过L曲线法确定,其本质是由离散点$\left( {{\rm{ln}}\left( {{{\left\| {\mathit{\boldsymbol{A\hat X}} - \mathit{\boldsymbol{L}}} \right\|}^2}} \right), {\rm{ln}}\left( {{{\left\| {\mathit{\boldsymbol{\hat X}}} \right\|}^2}} \right)} \right)$形成L曲线,然后确定相应的拐点,对应的正则化参数即为所求[16]

      定义如下GCV(generalized cros‐ validation)函数[17]

      $$ {\rm{GCV}} (k) = \frac{{{{\left\| {\mathit{\boldsymbol{A\hat X}} - \mathit{\boldsymbol{L}}} \right\|}^2}}}{{[ {\rm{tr}} (\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{A}}^ + })]}} $$ (14)

      GCV方法选取正则化参数的本质即求使式(14)达到最小的正则化参数。

      通常称上述方法为单参数Tikhonov正则化方法,原因在于只包含一个正则化参数。而多参数正则化方法的目标函数为[12]

      $$ \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} = {\left\| {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{AX}}} \right\|^2} + \sum\limits_{i = }^k {{{\left\| {{\mathit{\boldsymbol{B}}_i}\mathit{\boldsymbol{X}}} \right\|}^2}} + \beta {\left\| \mathit{\boldsymbol{X}} \right\|^2} $$ (15)

      其中正则化参数采用偏差原理求得:

      $$ {\left\| {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}^\delta }({\alpha _1}, {\alpha _2} \cdots {\alpha _k}, \beta )} \right\|^2} = c\delta , c \ge 1 $$ (16)

      于是构造出如下形式的模型函数:

      $$ m({\alpha _1}, {\alpha _2} \cdots {\alpha _k}, \beta ) = {\left\| \mathit{\boldsymbol{L}} \right\|^2} + \sum\limits_{i = 1}^k {\frac{{{C_i}}}{{{\alpha _i}}}} + \frac{D}{{\beta + T}} $$ (17)

      然后利用偏差原理迭代求解正则化参数,并令

      $$ \begin{array}{*{20}{l}} {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {G^m}({\alpha _1}, {\alpha _2} \cdots {\alpha _k}, \beta ) = m({\alpha _1}, {\alpha _2} \cdots {\alpha _k}, \beta ) - }\\ {\beta \frac{{\partial m({\alpha _1}, {\alpha _2} \cdots {\alpha _k}, \beta )}}{{\partial \beta }} - \sum\limits_{i = 1}^k {{\alpha _i}} \frac{{\partial m({\alpha _1}, {\alpha _2} \cdots {\alpha _k}, \beta )}}{{\partial {\alpha _i}}}} \end{array} $$ (18)

      如果(α1,lα2,lαklβ l)为第l次正则化参数解,可通过如下方程:

      $$ {G^m}({\alpha _{1, l + 1}}, {\alpha _{2, l + 1}} \cdots {\alpha _{j - 1, l + 1, }}, \alpha , {\alpha _{j + 1, l, }}{\alpha _{k, l}}, {\beta _l}) = {c^2}{\delta ^2} $$ (19)

      求解得到第k个正则化参数α = αjl + 1,其中j = 1,2…k - 1。因而,在每次迭代计算完毕α j(j = 1…k),同理可以通过如下方程:

      $$ {G^m}({\alpha _{1, l + 1}}, {\alpha _{2, l + 1}} \cdots {\alpha _{k, l}}, \beta ) = {c^2}{\delta ^2} $$ (20)

      计算得到β = βl + 1

      上述方法的设计和实现存在两个方面的问题:①k值取多大合适;②如何同时确定多个正则化参数α i。而且每次迭代需要单个确定正则化参数,然后通过此参数确定下一个参数,显然计算量相对过大。为此,本文结合SVD分解,设计相对简便的多参数正则化方法。

      对系数矩阵的特征值进行修正,以得到修正其条件数的目的,即

      $$ $\begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{A}}^*} = \\ \mathit{\boldsymbol{U}}\left( {\begin{array}{*{20}{c}} {{\lambda _1} + {\alpha _1}}&{}&{}&{}&0\\ {}&{{\lambda _2} + {\alpha _2}}&{}&{}&{}\\ {}&{}& \ddots &{}&{}\\ {}&{}&{}&{{\lambda _{t - 1}} + {\alpha _{t - 1}}}&{}\\ 0&{}&{}&{}&{{\lambda _t} + {\alpha _t}} \end{array}} \right){\mathit{\boldsymbol{V}}^{\rm{T}}} \end{array}$ $$ (21)

      于是,通过式(11)的广义逆解法,设计如下多参数正则化解:

      $$ {\mathit{\boldsymbol{\hat X}}_{{\rm{mul}}}} = {\mathit{\boldsymbol{A}}^{* + }}\mathit{\boldsymbol{b}} = \mathit{\boldsymbol{V}}{\left[ {\begin{array}{*{20}{c}} {{\lambda _1} + {\alpha _1}}&{}&{}&{}&0\\ {}&{{\lambda _2} + {\alpha _2}}&{}&{}&{}\\ {}&{}& \ddots &{}&{}\\ {}&{}&{}&{{\lambda _{t - 1}} + {\alpha _{t - 1}}}&{}\\ 0&{}&{}&{}&{{\lambda _t} + {\alpha _t}} \end{array}} \right)^{ - 1}}{\mathit{\boldsymbol{U}}^{\rm{T}}}\mathit{\boldsymbol{L}} $$ (22)

      对应的协方差矩阵为:

      $$ \mathit{\boldsymbol{D}}({\mathit{\boldsymbol{\hat X}}_{{\rm{mul}}}}) = \sigma _0^2\mathit{\boldsymbol{V}}{\left[ {\begin{array}{*{20}{c}} {{{({\lambda _1} + {\alpha _1})}^2}}&{}&{}&{}&0\\ {}&{{{({\lambda _2} + {\alpha _2})}^2}}&{}&{}&{}\\ {}&{}& \ddots &{}&{}\\ {}&{}&{}&{{{({\lambda _{t - 1}} + {\alpha _{t - 1}})}^2}}&{}\\ 0&{}&{}&{}&{{{({\lambda _t} + {\alpha _t})}^2}} \end{array}} \right)^{ - 1}}{\mathit{\boldsymbol{V}}^{\rm{T}}} $$ (23)

      对上述协方差矩阵求迹,得到:

      $$ {\rm{ tr}} [\mathit{\boldsymbol{D}}({\mathit{\boldsymbol{\hat X}}_{{\rm{ mul }}}})] = \sigma _0^2\sum\limits_{i = 1}^t {\frac{1}{{{{({\lambda _t} + {\alpha _i})}^2}}}} $$ (24)

      显然,基于式(22)的未知参数估计为有偏估计,即$\mathit{\boldsymbol{E}}(\mathit{\boldsymbol{\hat X}}) \ne \mathit{\boldsymbol{X}}$。于是有:

      $$ $\begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{Bias}}\left( {{{\mathit{\boldsymbol{\hat X}}}_{{\rm{mul}}}}} \right) = \\ \mathit{\boldsymbol{E}}\left( {\mathit{\boldsymbol{\hat X}}} \right) - \mathit{\boldsymbol{X}} = \left[ {\mathit{\boldsymbol{V}}{{\left( {\begin{array}{*{20}{c}} {\frac{{{\lambda _1} + {\alpha _1}}}{{{\lambda _1}}}}&{}&{}&{}&0\\ {}&{\frac{{{\lambda _2} + {\alpha _2}}}{{{\lambda _2}}}}&{}&{}&{}\\ {}&{}& \ddots &{}&{}\\ {}&{}&{}&{\frac{{{\lambda _{t - 1}} + {\alpha _{t - 1}}}}{{{\lambda _{t - 1}}}}}&{}\\ 0&{}&{}&{}&{\frac{{{\lambda _t} + {\alpha _t}}}{{{\lambda _t}}}} \end{array}} \right)}^{ - 1}}{\mathit{\boldsymbol{V}}^{\rm{T}}} - {\mathit{\boldsymbol{I}}_n}} \right]\mathit{\boldsymbol{X}} \end{array}$ $$ (25)

      结合如下MSE的计算公式[18]

      $$ {\rm{MSE}}({\mathit{\boldsymbol{\hat X}}_{{\rm{mul}}}}) = \mathit{\boldsymbol{E}}{({\mathit{\boldsymbol{\hat X}}_{{\rm{mul}}}} - \mathit{\boldsymbol{X}})^{\rm{T}}}({\mathit{\boldsymbol{\hat X}}_{{\rm{mul}}}} - \mathit{\boldsymbol{X}}) = {\rm{tr}}[D({\mathit{\boldsymbol{\hat X}}_{{\rm{mul}}}})] + {\rm{Bias}}{({\mathit{\boldsymbol{\hat X}}_{{\rm{mul}}}})^{\rm{T}}}{\rm{Bias}}({\mathit{\boldsymbol{\hat X}}_{{\rm{mul}}}}) $$ (26)

      于是多参数正则化解${\mathit{\boldsymbol{\hat X}}_{{\rm{mul}}}}$的MSE为:

      $$ {\rm{ MSE}} ({\mathit{\boldsymbol{\hat X}}_{{\rm{mul}}}}) = \sum\limits_{i = 1}^t {\left( {\frac{{\sigma _0^2}}{{{{({\lambda _i} + {\alpha _i})}^2}}} + {{\left( {\frac{{{\alpha _i}}}{{{\lambda _i} + {\alpha _i}}}} \right)}^2}{{(\mathit{\boldsymbol{v}}_i^{\rm{T}}\mathit{\boldsymbol{X}})}^2}} \right)} $$ (27)

      式中,vi为正交矩阵V的第i列。事实上,针对有偏估计,通常采用MSE作为评价标准[18],故目标是求正则化参数α i(i = 1,2…t),使得${\rm{MSE}} ({\mathit{\boldsymbol{\hat X}}_{{\rm{mul}}}})$达到极小值。对正则化矩阵,求${\rm{MSE}} ({\mathit{\boldsymbol{\hat X}}_{{\rm{mul}}}})$的偏导数并令其为零,得到:

      $$ $\begin{array}{*{20}{c}} {\frac{{\partial {\rm{MSE}}\left( {{{\mathit{\boldsymbol{\hat X}}}_{{\rm{mul}}}}} \right)}}{{\partial {\alpha _i}}} = }\\ {\left[ { - \frac{{2\sigma _0^2}}{{{{\left( {{\lambda _i} + {\alpha _i}} \right)}^3}}} + \frac{{2{\alpha _i}{\lambda _i}}}{{{{\left( {{\lambda _i} + {\alpha _i}} \right)}^3}}}{{\left( {\mathit{\boldsymbol{v}}_i^{\rm{T}}\mathit{\boldsymbol{X}}} \right)}^2}} \right] = }\\ {{\rm{}}\frac{{2\left[ {{\alpha _i}{\lambda _i}{{\left( {\mathit{\boldsymbol{v}}_i^{\rm{T}}\mathit{\boldsymbol{X}}} \right)}^2} - \sigma _0^2} \right]}}{{{{\left( {{\lambda _i} + {\alpha _i}} \right)}^3}}} = 0{\rm{}}} \end{array}$ $$ (28)

      进而得到第i个正则化参数为:

      $$ {\alpha _i} = \frac{{\sigma _0^2}}{{{\lambda _i}{{(\mathit{\boldsymbol{v}}_i^{\rm{T}}\mathit{\boldsymbol{X}})}^2}}} $$ (29)

      显然式(29)是存在的,原因在于λ i ≠ 0和v iT X ≠ 0。假设v iT X = 0,T = A* X = USV T X与系数矩阵是列满秩相矛盾,则具体的算法实现步骤如下:

      1)给定初值X (0)和(σ02)(0);

      2)针对第i次迭代,由计算得到的X (i)和(σ02) (i)计算正则化参数αi(i)(i = 1,2…t);

      3)通过计算出的正则化参数,由式(22)计算多参数正则化解$\mathit{\boldsymbol{\hat X}}_{{\rm{mul}}}^{\left( {i + 1} \right)}$;

      4)由多参数正则化解$\mathit{\boldsymbol{\hat X}}_{{\rm{mul}}}^{\left( {i + 1} \right)}$计算验后单位权方差(σ02) (i + 1)

      5)判断$\left\| {\alpha _{{\rm{ mul }}}^{(i + 1)} - \alpha _{{\rm{ mul }}}^{(i)}} \right\| < \varepsilon $是否满足,其中$\alpha _{{\rm{ mul }}}^{(i)} = {[\alpha _1^{(i)}, \alpha _2^{(i)} \cdots \alpha _t^{(i)}]^{\rm{T}}}$。如果满足条件,则计算结束;否则,返回步骤2)。

      定理当${\lambda _i} < \mathit{\boldsymbol{v}}_i^{\rm{T}}\mathit{\boldsymbol{X}}/{\sigma _0}{\rm{ 或者 }}{\lambda _i} > \mathit{\boldsymbol{v}}_i^{\rm{T}}\mathit{\boldsymbol{X}}/{\sigma _0}$和$0 < {\alpha _i} < - 2{\lambda _i}/[1 - \lambda _i^2{(\mathit{\boldsymbol{v}}_i^{\rm{T}}\mathit{\boldsymbol{X}})^2}/\sigma _0^2]$时,均方误差意义下多参数正则化解优于LS估计。

      通过式(26),可以得到LS估计的MSE:

      $$ {\rm{MSE}} ({\mathit{\boldsymbol{\hat X}}_{{\rm{LS}}}}) = {\rm{tr}} [D({\mathit{\boldsymbol{\hat X}}_{{\rm{LS}}}})] = \sigma _0^2\sum\limits_{i = 1}^t {\frac{1}{{\lambda _i^2}}} $$ (30)

      结合多参数正则化解的MSE,它们之间的差值为:

      $$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{MSE }}({{\mathit{\boldsymbol{\hat X}}}_{{\rm{LS}}}}) - {\rm{MSE }} ({{\mathit{\boldsymbol{\hat X}}}_{{\rm{mul}}}}) = \\ \sigma _0^2\mathop \sum \limits_{i = 1}^t \frac{1}{{\lambda _i^2}} - \mathop \sum \limits_{i = 1}^t \left[ {\frac{{\sigma _0^2}}{{{{({\lambda _i} + {\alpha _i})}^2}}} + {{(\frac{{{\alpha _i}}}{{{\lambda _i} + {\alpha _i}}})}^2}{{(\mathit{\boldsymbol{v}}_i^{\rm{T}}\mathit{\boldsymbol{X}})}^2}} \right] = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sigma _0^2\sum\limits_{i = 1}^t {\frac{{{{({\lambda _t} + {\alpha _i})}^2} - \lambda _i^2 - \lambda _i^2\alpha _i^2{{(\mathit{\boldsymbol{v}}_i^{\rm{T}}\mathit{\boldsymbol{X}})}^2}/\sigma _0^2}}{{\lambda _i^2{{({\lambda _i} + {\alpha _i})}^2}}}} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sigma _0^2\sum\limits_{i = 1}^t {\frac{{\alpha _i^2 + 2{\alpha _i}{\lambda _i} - \lambda _i^2\alpha _i^2{{(\mathit{\boldsymbol{v}}_i^{\rm{T}}\mathit{\boldsymbol{X}})}^2}/\sigma _0^2}}{{\lambda _i^2{{({\lambda _i} + {\alpha _i})}^2}}}} \end{array} $$ (31)

      当满足${\lambda _i} < {\sigma _0}/\mathit{\boldsymbol{v}}_i^{\rm{T}}\mathit{\boldsymbol{X}}$时,${({\lambda _i} + {\alpha _i})^2} - \lambda _i^2 - \lambda _i^2\alpha _i^2{(\mathit{\boldsymbol{v}}_i^{\rm{T}}\mathit{\boldsymbol{X}})^2}/\sigma _0^2 > 0$;而当满足${\lambda _i} > {\sigma _0}/\mathit{\boldsymbol{v}}_i^{\rm{T}}\mathit{\boldsymbol{X}}$和$0 < {\alpha _i} < - 2{\lambda _i}/[1 - \lambda _i^2{(\mathit{\boldsymbol{v}}_i^{\rm{T}}\mathit{\boldsymbol{X}})^2}/\sigma _0^2]$时,${({\lambda _i} + )^2} - \lambda _i^2 - \lambda _i^2\alpha _i^2{(\mathit{\boldsymbol{v}}_i^{\rm{T}}\mathit{\boldsymbol{X}})^2}/\sigma _0^2 > 0$,从而有:

      $$ {\rm{MSE}} ({\mathit{\boldsymbol{\hat X}}_{{\rm{LS}}}}) - {\rm{MSE}} ({\mathit{\boldsymbol{\hat X}}_{{\rm{mul}}}}) > 0 $$ (32)

      得证。

    • 为了验证多参数正则化方法在航空重力向下延拓数据处理中的可靠性,本文基于EGM2008地球重力场模型设计航行高度为5 km的仿真试验。通常来说,向下延拓的病态性程度跟航行高度有关,故本文选择相对适宜的高度,许多文献也采用此飞行高度。将EGM2008模型计算得到的地面重力异常值当作参考值,然后将其向上延拓到空中。仿真数据区域大小为1.5°× 1.5°,格网间距为2′×2′。为了使模拟数据更加符合真实场景,在模拟的航空重力数据中加入中误差为3 mGal和5 mGal的随机误差,共设计以下4种方案进行航空重力向下延延拓计算:

      方案1:LS方法;

      方案2:基于GCV方法选取正则化参数的Tikhonov正则化方法;

      方案3:基于L曲线方法选取正则化参数的Tikhonov正则化方法;

      方案4:多参数正则化方法。

      这里选用了两个地势不同的区域进行数据仿真,具体位置如图 1所示。其中第一个区域地势相对较高,第二个区域相对平缓,地势相对较低。在计算中发现,采用LS方法时,其法方程中法矩阵的条件数大于10 000,存在严重的病态性,得到的解不能够如实地反映真实情况。在3 mGal和5 mGal误差条件下,对于区域1、区域2,各种方案向下延拓的统计结果分别列于表 1表 2。从表 1中结果可以看出,LS估计的向下延拓结果严重失真,而正则化解则有效地削弱了法方程的病态性影响,得到了相对可靠的结果。3 mGal误差条件下,基于L曲线法和GCV方法选取正则化参数的两种正则化方法的向下延拓结果相当,均方根(root mean square,RMS)均在8 mGal左右,而本文建议的多参数正则化方法则在5 mGal,精度得到显著提高。从图 2中给出的各种方案的向下延拓结果与参考值之间的距离绝对值来看,构造多个正则化参数的正则化方法在大多数情形下,误差为1~2 mGal,向下延拓结果明显优于基于L曲线法和GCV方法的正则化方法。为了进一步分析在不同误差条件下各种方法的稳健性,继续加入5 mGal的随机误差,此时方案2和方案3的RMS分别在10 mGal和8 mGal左右,而方案4处于6 mGal,说明噪声越大,越会影响各种方法的向下延拓结果。

      图  1  数据仿真区域地形

      Figure 1.  Topography of Data Simulation Area

      表 1  各种方案大地水准面重力异常与参考值之间距离绝对值的统计结果(区域1) /mGal

      Table 1.  Statistical Results of Computing Geoid Gravity Gravity Anomaly and Reference Value for Each Scheme(First Area) /mGal

      误差 方案 最大值 最小值 平均值 RMS
      3 1 520.3 0.2 109.2 134.7
      2 29.1 0.0 6.7 8.5
      3 27.8 0.0 6.4 8.1
      4 18.3 0.0 4.2 5.3
      5 1 804.3 0.8 194.5 245.9
      2 43.8 0.0 8.2 10.4
      3 36.0 0.0 6.9 8.7
      4 23.0 0.0 5.3 6.8

      表 2  各种方案大地水准面重力异常与参考值之间距离绝对值的统计结果(区域2)/mGal

      Table 2.  Statistical Results of Computing Geoid Gravity Gravity Anomaly and Reference Value for Each Scheme(Second Area) /mGal

      误差 方案 最大值 最小值 平均值 RMS
      3 1 866.8 0.5 164.6 209.9
      2 11.9 0.0 2.5 3.1
      3 9.9 0.0 2.2 2.7
      4 5.6 0.0 1.2 1.6
      5 1 1 027.1 0.7 259.0 327.0
      2 16.1 0.0 3.3 4.2
      3 11.7 0.0 2.7 3.3
      4 3.3 0.0 1.9 2.4

      图  2  3 mGal误差条件下航空重力向下延拓值与地面参考重力值的差值绝对值(区域1)

      Figure 2.  Absolute Difference of Down‐Continue Results and Reference Value with 3 mGal Random Error (First Area)

      图 3可以得出,多参数正则化方法的大多数向下延拓结果与参考值之间的距离绝对值在0~2 mGal,优于常规的正则化方法,其优越性得到了充分体现。而从区域2的向下延拓分析来看,LS的向下延拓结果仍然较差。而其他的正则化方法则较为有效地削弱了系数矩阵病态性的影响,得到了较好的延拓结果。综合表 2图 4图 5的结果来看,方案2的RMS在4 mGal左右,方案3的RMS在3 mGal左右,而多参数正则化方法的RMS在2 mGal左右。可以看出,多参数正则化方法向下延拓结果相对较优,能够抵御系数矩阵的病态性影响。

      图  3  5 mGal误差条件下航空重力向下延拓值与地面参考重力值的差值绝对值(区域1)

      Figure 3.  Absolute Difference of Down‐Continue Results and Reference Value with 5 mGal Random Error (First Area)

      图  4  3 mGal误差条件下航空重力向下延拓值与地面参考重力值的差值绝对值(区域2)

      Figure 4.  Absolute Difference of Down‐Continue Results and Reference Value with 3 mGal Random Error (Second Area)

      图  5  5 mGal误差条件下航空重力向下延拓值与地面参考重力值的差值绝对值(区域2)

      Figure 5.  Absolute Difference of Down‐Continue Results and Reference Value with 5 mGal Random Error (Second Area)

    • 1)基于Poisson积分问题进行航空重力向下延拓,其解算方程的法矩阵的条件数较大,呈现出严重的病态性,采用LS估计方法进行法方程解算,延拓结果较差。

      2)本文建议了一种简单和方便计算的多参数正则化方法,并采用MSE标准设计了相应的正则化参数选取算法。在均方误差意义下,给出了多参数正则化方法优于最小二乘估计的条件。

      3)基于EGM2008重力场模型进行了仿真试验,并比较了单参数正则化方法的L曲线法和GCV方法。数值结果表明,多参数正则化方法能进一步提高向下延拓的精度和可靠性,验证了该方法在航空重力向下延拓的可应用性和可行性。

参考文献 (18)

目录

    /

    返回文章
    返回