留言板

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

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

基于χ2准则的磁梯度张量3D聚焦反演方法

李金朋 张英堂 范红波 李志宁 尹刚 刘敏

李金朋, 张英堂, 范红波, 李志宁, 尹刚, 刘敏. 基于χ2准则的磁梯度张量3D聚焦反演方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(2): 255-261. doi: 10.13203/j.whugis20160063
引用本文: 李金朋, 张英堂, 范红波, 李志宁, 尹刚, 刘敏. 基于χ2准则的磁梯度张量3D聚焦反演方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(2): 255-261. doi: 10.13203/j.whugis20160063
LI Jinpeng, ZHANG Yingtang, FAN Hongbo, LI Zhining, YIN Gang, LIU Min. Three-dimensional Focusing Inversion of Magnetic Gradient Tensor Data Based on the χ2 Principle[J]. Geomatics and Information Science of Wuhan University, 2018, 43(2): 255-261. doi: 10.13203/j.whugis20160063
Citation: LI Jinpeng, ZHANG Yingtang, FAN Hongbo, LI Zhining, YIN Gang, LIU Min. Three-dimensional Focusing Inversion of Magnetic Gradient Tensor Data Based on the χ2 Principle[J]. Geomatics and Information Science of Wuhan University, 2018, 43(2): 255-261. doi: 10.13203/j.whugis20160063

基于χ2准则的磁梯度张量3D聚焦反演方法

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

国家自然科学基金 51305454

详细信息
    作者简介:

    李金朋, 博士生, 主要从事磁测数据处理及解释方法研究。18626648671@163.com

  • 中图分类号: P318

Three-dimensional Focusing Inversion of Magnetic Gradient Tensor Data Based on the χ2 Principle

Funds: 

The National Natural Science Foundation of China 51305454

More Information
    Author Bio:

    LI Jinpeng, PhD candidate, specializes in magnetic data processing and interpretation method. E-mail: 18626648671@163.com

图(7) / 表(2)
计量
  • 文章访问数:  1078
  • HTML全文浏览量:  67
  • PDF下载量:  242
  • 被引次数: 0
出版历程
  • 收稿日期:  2016-05-23
  • 刊出日期:  2018-02-05

基于χ2准则的磁梯度张量3D聚焦反演方法

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

    国家自然科学基金 51305454

    作者简介:

    李金朋, 博士生, 主要从事磁测数据处理及解释方法研究。18626648671@163.com

  • 中图分类号: P318

摘要: 针对铁磁性物质反演中正则化参数自适应选择的问题,提出了基于χ2准则的磁梯度张量3D聚焦反演方法。利用深度加权矩阵和最小支撑矩阵对经典Tikhonov正则化理论框架下的反演模型进行约束得到目标函数,避免了由于反演参数多于采集点数而导致反演解的多解性,并有效解决了核函数随深度增大而快速衰减的问题。通过对目标函数进行迭代奇异值分解获得最佳物性参数,并根据χ2准则自适应地确定目标函数在迭代过程中的正则化参数,提高了迭代速度和求解精度。仿真和实验结果表明:该方法能准确还原磁性异常体的轮廓形态,具有较好的模型分辨率。

English Abstract

李金朋, 张英堂, 范红波, 李志宁, 尹刚, 刘敏. 基于χ2准则的磁梯度张量3D聚焦反演方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(2): 255-261. doi: 10.13203/j.whugis20160063
引用本文: 李金朋, 张英堂, 范红波, 李志宁, 尹刚, 刘敏. 基于χ2准则的磁梯度张量3D聚焦反演方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(2): 255-261. doi: 10.13203/j.whugis20160063
LI Jinpeng, ZHANG Yingtang, FAN Hongbo, LI Zhining, YIN Gang, LIU Min. Three-dimensional Focusing Inversion of Magnetic Gradient Tensor Data Based on the χ2 Principle[J]. Geomatics and Information Science of Wuhan University, 2018, 43(2): 255-261. doi: 10.13203/j.whugis20160063
Citation: LI Jinpeng, ZHANG Yingtang, FAN Hongbo, LI Zhining, YIN Gang, LIU Min. Three-dimensional Focusing Inversion of Magnetic Gradient Tensor Data Based on the χ2 Principle[J]. Geomatics and Information Science of Wuhan University, 2018, 43(2): 255-261. doi: 10.13203/j.whugis20160063
  • 近年来,磁力勘探主要应用在军事侦察、水文及工程地质勘探等领域。运用磁探测技术可以对地下或者水下小尺度磁性目标进行定位与识别[1],为下一步的工作提供数据支持。反演方法是解释实测数据的重要一步,其主要目的是利用地表磁测数据来评估地下未知磁异常体的轮廓形态[2]。在磁测数据线性反演方法中,假设观测面以下的实测区域可以离散化为具有固定物性参数值的规则立方体模型,由于地下待求未知参数远远多于观测数据点数[3],因此会导致线性方程求解的模糊性,难以得到方程的最优解。测量数据的噪声和固有的场源不唯一性,导致磁测数据的反演是病态多解的。为了降低测量误差等因素的影响,常常使用L2或者Tikhonov正则化思想[4, 5]。正则化思想的两个重要的方面是模型约束器和正则化参数的选择,前者可以使数据获得稳定解,后者是数据项和正则项之间的折中系数。在反演问题中正则化参数的选择方法是多年来研究的热点问题之一。确定正则化参数的主要方法有LC曲线(L-curve, LC)、广义交叉验证(generalized cross validation, GCV)及莫罗佐夫偏差原则(Morozov discrepancy principle, MDP)[6-8]。其中LC曲线法计算量较大,需要计算多个正则化参数来确定最优解; GCV方法难以保证在在每次迭代过程中都能获得GCV功能矩阵的极小值;相比前两种方法,MDP目前应用较为广泛,但是在求解过程中会产生过高的正则化参数值,导致结果过于光滑[9-11]

    针对上述反演过程存在的问题,本文提出了基于χ2准则的磁梯度张量3D聚焦反演方法。通过向目标函数中加入深度加权矩阵和最小支撑矩阵,避免了反演解的不稳定性;采用奇异值分解法将反演解转化到合理的值域范围内;利用χ2准则对正则化参数进行自适应选择,提高了迭代速度和求解精度。

    • 磁梯度张量正演公式如下:

      $$ \boldsymbol{d} = \boldsymbol{Am},p \ll n $$ (1)

      式中,核函数矩阵A的元素Aij为第i(i=1, 2, 3…p)个观测点观测的第j(j=1, 2, 3…n)个单元的响应;向量dRp为观测数据矩阵;向量mRn为待测模型参数矩阵。在磁梯度张量数据实际测量中,由于操作误差等影响导致实测数据中包含一定的噪声。设eRp代表实测数据中的误差,则式(1)表示为:

      $$ {\boldsymbol{d}_{\rm obs}} = \boldsymbol{Am} + \boldsymbol{e} $$ (2)

      根据式(2)可知,磁梯度张量反演的过程就是利用包含噪声的实际测量数据dobs求解模型参数m。式(2)表示的病态反问题可以表示为Tikhonov最小二乘公式:

      $$\begin{align} \boldsymbol{m}\left( \alpha \right) = &\arg \mathop {\min }\limits_m \{ {\left\| {{\boldsymbol{W}_d}\left( {\boldsymbol{Am} - {\boldsymbol{d}_{\rm obs}}} \right)} \right\|}^2_2 + \\ & {\alpha ^2}\left\| {\boldsymbol{D}(\boldsymbol{m} - {\boldsymbol{m}_{\rm apr}})} \right\|^2_2 \}\end{align} $$ (3)

      式中,‖Wd(Amdobs)‖22为加权数据保真度,‖D(mmapr)‖22为正则项,Am为预测数据向量,Wd=diag(1/η1, …, 1/ηp)为数据加权矩阵,ηi为第i个基准噪声的标准偏差,D为正则化矩阵,mapr为模型参数先验信息向量,α为正则化参数。设Ã=WdA,$\boldsymbol{\tilde d} $obs=WddobsZ=mmapr。转换后式(3)表示为:

      $$ \begin{array}{c} \boldsymbol{Z}\left( \alpha \right) = \arg \mathop {\min }\limits_y \left\{ {{\left\| {\boldsymbol{\tilde A}\boldsymbol{m} - \boldsymbol{\tilde r}} \right\|}^2_2 + {\alpha ^2}{\left\| {\boldsymbol{DZ}} \right\|}^2_2} \right\},\\ \boldsymbol{\tilde r }= ({{\boldsymbol{\tilde d}}_{\rm obs}} - \boldsymbol{\tilde A}{\boldsymbol{m}_{\rm apr}}) \end{array} $$ (4)

      此时m(α)仅由正则化参数α决定,用矩阵Ã(α)表示为:

      $$ \begin{align} \boldsymbol{Z}\left( \alpha \right) &= {({\boldsymbol{\tilde A}^{\rm T}}\boldsymbol{\tilde A} + {\alpha ^2}{\boldsymbol{D}^{\rm T}}\boldsymbol{D})^{ - 1}}{\boldsymbol{\tilde A}^{\rm T}}\boldsymbol{\tilde r} = \boldsymbol{\tilde A}\left( \alpha \right)\boldsymbol{\tilde r},\\ &\boldsymbol{\tilde A}\left( \alpha \right) = {({\boldsymbol{\tilde A}^{\rm T}}\boldsymbol{\tilde A} + {\alpha ^2}{\boldsymbol{D}^{\rm T}}\boldsymbol{D})^{ - 1}}{\boldsymbol{\tilde A}^{\rm T}} \end{align} $$ (5)
      $$ \boldsymbol{m}\left( \alpha \right) = {\boldsymbol{m}_{\rm apr}} + \boldsymbol{Z}\left( \alpha \right) = {\boldsymbol{m}_{\rm apr}} + \boldsymbol{\tilde A}\left( \alpha \right)\boldsymbol{\tilde r} $$ (6)

      Hansen指出矩阵D为可逆矩阵[10],因此:

      $$ \begin{array}{c} ({\boldsymbol{\tilde A}^{\rm T}}\boldsymbol{\tilde A} + {\alpha ^2}{\boldsymbol{D}^T}\boldsymbol{D}) = \\ {\boldsymbol{D}^T}({({\boldsymbol{D}^T})^{ - 1}}{\boldsymbol{\tilde A}^{\rm T}}\boldsymbol{\tilde A}{\boldsymbol{D}^{ - 1}} + {\alpha ^2}{I_n})\boldsymbol{D} \end{array} $$ (7)

      定义右预处理矩阵$ \boldsymbol{\mathop A}\limits^ \approx$=ÃD-1,正则化反演矩阵$ \boldsymbol{\mathop A}\limits^ \approx(\alpha)=({\boldsymbol{\mathop A}\limits^ \approx }^{\rm T}\boldsymbol{\mathop A}\limits^ \approx+\alpha^2I_n)^{-1}{\boldsymbol{\mathop A}\limits^ \approx}^{\rm T}$,此时J(α)=DZ(α),则式(4)表示为:

      $$ \boldsymbol{J}\left( \alpha \right) = \arg \mathop {\min }\limits_J \left\{ {{\left\| {\boldsymbol{\mathop A}\limits^ \approx \boldsymbol{J} - \boldsymbol{\tilde r}} \right\|}^2_2 + {\alpha ^2}{\left\| \boldsymbol{J} \right\|}^2_2} \right\} $$ (8)
      $$ \boldsymbol{m}\left( \alpha \right) = {\boldsymbol{m}_{\rm apr}} + {\boldsymbol{D}^{ - 1}}\boldsymbol{J}\left( \alpha \right) $$ (9)

      正则化矩阵D(k)Rn×n可以表示为D(k)=We(k)Wdepth。文献[6]给出深度加权矩阵Wdepth=diag(1/(zj+ξ)β)来抵消核函数随深度增加所带来的衰减,其中zj为地下第j个剖分单元的中心埋深,ξ>0是为了避免出现奇点而引入的数值,β是可以控制加权函数大小的正数。Vatankhah指出,利用最小支撑矩阵We(k)=diag((m(k)m(k-1))+ε2)-1/2获得的反演结果具有较好的非光滑特征,能够提高块状异常体的分辨率,k>0,W(0)=I, m(0)=maprε>0为聚焦参数。通过正则化矩阵对正则项进行约束,减小模型解的多解性,提高了反演解的稳定性。在反演过程中,为了得到与实际情况更加吻合的物性参数分布,利用惩罚函数在反演迭代过程中去除病态解,在迭代过程中设置物性参数上下限约束函数[mmin, mmax],对物性参数进行约束使反演结果控制在合理的范围内:

      $$ \left\{ \begin{array}{c} m\left( j \right) = {m_{\max}},若m\left( j \right) > {m_{\max}}\\ m\left( j \right) = {m_{\min}},若m\left( j \right) < {m_{\min}}\\ j = 1,2,3, \ldots ,n; \end{array} \right.{\rm{ }} $$ (10)

      设置终止迭代条件:当χComputed2=‖(dobs)i-(dpre)i/ηi22p+$\sqrt {2p} $或达到最大迭代次数k时迭代终止。迭代过程为:

      $$ \boldsymbol{J}({\alpha ^{(k + 1)}}) = ({\boldsymbol{\mathop A}\limits^ \approx}^{\rm T}\boldsymbol{\mathop A}\limits^ \approx + {({\alpha ^{(k)}})^2}{\boldsymbol{I}_n})^{ - 1}\boldsymbol{\mathop A}\limits^ \approx{\boldsymbol{\tilde r}^{(k)}} $$ (11)
      $$ {\boldsymbol{m}^{(k + 1)}} = {\boldsymbol{m}^{(k)}} + {({\boldsymbol{D}^{\left( {k + 1} \right)}})^{ - 1}}\boldsymbol{J}({\alpha ^{(k + 1)}}) $$ (12)

      在求解J(α)时,对右处理矩阵$ \boldsymbol{\mathop A}\limits^ \approx$∈Rp×n进行奇异值分解(singular value decomposition,SVD)。矩阵$ \boldsymbol{\mathop A}\limits^ \approx$可以表示为$ \boldsymbol{\mathop A}\limits^ \approx$=UVT,其中,UV分别为p×pn×n的正交矩阵。∑=diag([δ1, δ2, …δp]),其中奇异值δ1δ2≥…≥δp>0。为了防止高频观测误差被放大,引入正则化滤波因子fi(α)=$\frac{{\delta ^2_i}}{{\delta ^2_i + {\alpha ^2}}} $, 1≤ip。因此,对式(11)进行奇异值分解得:

      $$ \boldsymbol{J}\left( \alpha \right) = \sum\limits_{i = 1}^p {\frac{{\delta ^2_i}}{{\delta ^2_i + {\alpha ^2}}}\frac{{\boldsymbol{u}^{\rm T}_i\boldsymbol{\tilde r}}}{{{\delta _i}}}} {\boldsymbol{v}_i} = \sum\limits_{i = 1}^p {{f_i}\left( \alpha \right)\frac{{\boldsymbol{u}^{\rm T}_i\boldsymbol{\tilde r}}}{{{\delta _i}}}{\boldsymbol{v}_i}} $$ (13)

      式中,uivi分别为矩阵UV的第i列。

    • 徐新禹等讨论了在Tikhonov正则化方法中如何有效的选择正则化参数[12, 13]。Zhdanov提出一种自适应正则化参数的选择方法,定义α(k)=α(1)q(k)q为衰减因子[14]。但是这种方法不能保证每次迭代得到的正则化参数为最佳解,影响求解精度。还有一种常用的方法是MDP法,但Kilmer指出该方法会导致反演结果过度光滑[15]。因此,本文提出基于χ2准则的正则化参数选择方法。

      χ2准则实际上是MDP方法的一种扩展算法。Mead指出MDP是基于后验策略的偏差原理[11],通过利用先验误差水平信息,使得实测数据与预测数据的残差为最小值的正则化参数即为最优正则化参数,MDP实际上就是遵循自由度为p-nχ2分布,可以表示为:

      $$ \left\| {\boldsymbol{\mathop A}\limits^ \approx \boldsymbol{J}\left( \alpha \right) - \boldsymbol{\tilde r}} \right\|^2_2 = p - n,p \ge n $$ (14)

      由式(14)可知,当pn时,上述方法是不可用的。而基于χ2准则的正则化参数确定方法是通过确定Tikhonov正则化目标函数的极小值来确定最优正则化参数,且无需误差先验信息。该方程遵循自由度为pχ2分布,可以表示为:

      $$ \left\| {{\boldsymbol{W}_d}\left( {\boldsymbol{Am} - {\boldsymbol{d}_{\rm obs}}} \right)} \right\|^2_2 + {\alpha ^2}\left\| {\boldsymbol{D}\left( {\boldsymbol{m} - {\boldsymbol{m}_{\rm apr}}} \right)} \right\|^2_2 = p $$ (15)

      定义一个单直立长方体理论模型,纵横向切片图如图 1(a)(b)所示。将地下待测空间划分为22×22×10=4 840个单元格,每个单元格均为边长为0.1 m的正方体。在仅考虑感应磁化条件下,异常体的磁化强度为20 A/m,磁倾角为55°,磁偏角为-18°。根据d=Am计算右端项的理论值,并加入0.01×std(d, 1)×rand(484, 1)的误差作为实测数据。定义初始正则化参数α(1)=(n/p)γ(max(δi)/mean(δi)),其中0≤γ≤2,δi为对应核函数矩阵$\boldsymbol{\mathop A}\limits^ \approx $的奇异值,定义模型重构相对误差rms=‖mexactm(k)2/‖mexact2,当χComputed2≤506或最高迭代次数K=100时迭代终止。本文仅选择垂直向和水平向分辨率相对较高的Bzz分量进行方法的有效性验证[16-17], 反演结果如图 1(c)~1(f)所示。

      图  1  不同迭代原则下长方体反演结果

      Figure 1.  Cuboid Inversion Results in Different Iterative Principles

      从反演结果可以看出,两种方法都使得反演结果具有较好的横向和纵向分辨率,但是从表 1可知,基于χ2准则的正则化参数确定方法比MDP迭代次数少,求得的反演结果精度更高。由图 2可以看出,正则化参数初始值很大,表明在迭代初期正则项起到较大的约束作用,随后迅速衰减,最后随迭代次数增加趋于稳定。

      表 1  不同迭代原则下各部分参数值

      Table 1.  Parameter Values at Different Iterations Principle

      迭代原则 α(1), γ=2 相对误差 迭代次数
      χ2准则 112.6 0.7018 6
      MDP 112.6 0.7096 20

      图  2  正则化参数随迭代次数的变化规律

      Figure 2.  Regularization Parameter Changes with the Number of Iterations

    • 将地下待测空间划分为22×22×10=4 840个单元格,每个单元格均为边长为0.1 m的正方体,地下正演模型如图 3所示。在仅考虑感应磁化的条件下,磁倾角为55°,磁偏角为-16°,3个不同位置、不同形态的异常体的具体位置及参数如表 2。如图 4所示,在反演时加入0.001×std(d, 1)×rand(484, 1)的误差作为实测数据,并利用Bzz分量进行反演。反演过程中设置最高迭代次数为15次,设置物性参数阈值为m(k)>8A/m,所得反演结果的空间分布如图 5所示。从图中结果可以看出:该反演方法能够较好的反映出不同位置、不同物性参数模型的轮廓形态,具有较好的横向分辨率,纵向上随深度增加模型分辨率有所下降,存在边缘效应,这是由于重构模型的分辨率不可避免的随深度增加而下降,但是能够大致分布在相应的位置。

      图  3  理论模型图

      Figure 3.  Theoretical Model Figure

      表 2  异常体参数

      Table 2.  Abnormal Body Parameters

      异常体 X方向位置/m Y方向位置/m Z方向位置/m 磁化强度M(A·m-1)
      异常体1 1.55~1.85 0.35~1.75 0.1~0.4 20
      异常体2 0.25~0.55 0.75~1.75 0.3~0.6 30
      异常体3 1.35~1.75 0.25~0.95 0.2~0.4 40

      图  4  Bzz理论数据与高斯噪声

      Figure 4.  Bzz Theoretical data and Gaussian Noise

      图  5  反演结果平面图与反演结果三维立体图

      Figure 5.  Slice of Inversion Result and 3D Inversion Result

    • 测区位于石家庄某地,测区内放置两个长方体铁块。在地理坐标系下铁块1位于测区东南方向,中心坐标为(1.55, 0.5, 0.16),块体尺寸为0.2 m×0.4 m×0.05 m,铁块2位于测区西北方向,中心坐标为(0.65, 1.2, 0.27),块体尺寸为0.2 m×0.4 m×0.1 m。实验装置为十字形传感器阵列结构,利用Bartington公司的三轴磁通门传感器搭建磁梯度张量探头。图 6(a)展示了测量装置及铁块的实际位置分布,图 6(b)Bzz分量实际测量数据。将待测空间划分为22×22×10=4 840个单元格,每个单元格均为边长为0.1 m×0.1 m×0.1 m的长方体,地面的观测点为22×22=484个网格。分别采用MDP准则和χ2准则对磁异常进行反演,得到的反演结果如图 7所示。图 7(a)7(c)为模型体反演结果平面图,图 7(b)7(d)为物性参数阈值m(k)>5 A/m时的三维反演结果。其中利用χ2准则进行计算需要迭代12次,而利用MDP方法进行迭代需要17次。根据反演结果可以看出,相比MDP准则,χ2准则的反演结果的横向和纵向模型分辨率较好,能够更好地确定异常体的轮廓及位置,进一步证明了本文方法的有效性。

      图  6  实验装置及实测Bzz数据

      Figure 6.  Experimental Devices and Bzz Data

      图  7  反演结果

      Figure 7.  Inversion Result

    • 1) 利用深度加权矩阵和最小支撑矩阵对经典Tikhonov正则化理论框架下的反演模型进行约束,提高了反演解的稳定性;

      2) 采用奇异值分解法(SVD)对核函数矩阵进行分解,利用物性参数上下限函数将反演解转换到合理的值域范围内;

      3) 利用χ2准则对迭代过程中的正则化参数进行选择,实现了正则化参数选择的自适应功能,提高了迭代速度和求解精度;

      4) 虽然本文方法能够使获得的物性参数值分布在相应的位置,但是重构模型的分辨率不可避免地随着深度的增加而下降,存在一定的边缘效应。

参考文献 (17)

目录

    /

    返回文章
    返回