CN113466954B - 一种自反馈类正则化反演方法 - Google Patents

Abstract

本发明提供了一种自反馈类正则化反演方法，步骤S1：将平面磁异常数据进行化极并成图；步骤S2：根据化极后的磁异常平面图进行圈定，并对圈定数据进行插值；步骤S3：寻找极大值的平面位置，将此平面位置作为反演的初始位置；步骤S4：开始反演迭代过程，首先给初始模型赋值，然后对模型参数进行迭代更新，同时计算拟合差；迭代完毕后，得到输出反演的结果；步骤S5：得到反演模型参数以及吻合度；吻合度满足要求，则输出结果；吻合度不满足，则返回步骤S2，重新圈定数据，直到输出吻合度满足要求为止。本发明的新方法有效避免观测数据中噪声的影响，使反演过程更加稳健，更容易跳出局部极值，得出的反演结果更加接近实际位置。

Description

一种自反馈类正则化反演方法

技术领域

本发明涉及磁法勘探技术领域，具体涉及一种自反馈类正则化反演方法。

背景技术

磁法勘探在矿产勘查、地质勘探、深部构造、地下金属勘探中具有重要的作用。其中，磁偶极子通常可以用来描述多种物体的磁异常特征，因此，针对磁偶极子位置的反演具有重要的实际应用意义。

在磁偶极子反演中，由于磁测干扰往往比较强烈，反演容易受到噪声的影响变得不稳定，容易陷入局部极值。同时由于噪声存在，往往有很多模型能拟合观测数据。如果只是运用观测数据与正演响应数据的拟合差作为反演中最优化的目标，那么反演过程中，观测数据本身的噪声会对反演过程产生很大的影响。因为在反演过程中，如果只是不断减小正演模型数据与观测数据的拟合差，磁干扰强时反演结果往往不可靠。因此，如何改善偶极子反演中的这一问题，减少观测数据中误差对反演结果的影响，在偶极子反演应用中具有重要意义。

在最新的稳健高斯牛顿法反演中，一般采用下式进行迭代更新模型参数：

其中m^(k)，m^(k+1)分别为第k次和k+1次迭代的模型参数，J_k为第k次迭代的雅各比矩阵，即正演模型数据对反演模型的偏导数矩阵，r_k为第k次迭代时的拟合差，

表示J_k的转置矩阵。

这种高斯牛顿法需要要求J_k矩阵列满秩，同时优化过程中仅与拟合差有关，容易受到数据噪声的影响；同时存在容易陷入局部最小值，对数据过拟合的问题。

综上所述，急需一种自反馈类正则化反演方法以解决现有技术中存在的问题。

发明内容

本发明目的在于提供一种自反馈类正则化反演方法，旨在解决现有反演方法中存在容易受数据噪声影响，同时存在容易陷入局部最小值，对数据过拟合的问题，具体技术方案如下：

一种自反馈类正则化反演方法，包括以下步骤：

步骤S1：将采集的平面磁异常数据进行化极并成图；

步骤S2：根据化极后的磁异常平面图进行圈定，并对圈定数据进行插值；

步骤S3：在化极并插值后的数据中寻找极大值的平面位置，并将此平面位置作为反演的初始位置；

步骤S4：开始反演迭代过程，首先给初始模型赋值，然后对模型参数进行迭代更新，同时计算拟合差；迭代完毕后，得到输出反演的结果；

步骤S5：反演迭代结束，得到反演模型参数以及吻合度；如果吻合度满足要求，则输出结果；如果吻合度不满足，则返回步骤S2，重新圈定数据，直到输出吻合度满足要求为止。

以上技术方案中优选的，所述步骤S4中，采用式(1)和式(2)对模型参数进行迭代更新：

其中，m^(k)、m^(k+1)分别为第k次和k+1次迭代的模型；J_k为第k次迭代的雅各比矩阵，即正演模型数据对反演模型的偏导数矩阵；r_k为第k次迭代时的拟合差；

表示J_k的转置矩阵；I为单位矩阵；λ_k为第k次迭代的调节因子，为了与传统的超参数化反演对应称之为正则化因子。

以上技术方案中优选的，在偶极子反演中，选取的反演目标函数如式(3)：

其中，m为偶极子需要反演的模型；Φ_d(m)为正演模型数据与观测数据之间的拟合差的二范数，即：

其中，F[m]为模型正演数据矢量；

为模型正演数据矢量的第i个数据；d^obs为观测数据矢量；

为观测数据矢量的第i个数据；W_d代表的是观测误差的对角矩阵，

是观测误差的对角矩阵的第i行第i列的数据，N表示观测数据的总量；

根据式(1)和式(2)对模型参数不断修改迭代，直至完成设定的迭代次数。

以上技术方案中优选的，正演模型数据通过偶极子正演得到，具体是通过式(5)得到磁总场：

ΔT＝H_axL₀+H_ayM₀+Z_aN₀ (5)，

其中，L₀＝cosI₀cosA₀，M₀＝cosI₀sinA₀，N₀＝sinI₀；I₀为地磁场的磁偏角，A₀为地磁场的磁倾角，ΔT为磁总场，H_ax、H_ay和Z_a分别为磁场的x、y和z分量上的大小。

以上技术方案中优选的，其中H_ax、H_ay和Z_a的计算方式如下：

设总磁化的磁偏角为I、磁倾角为A，磁偶极子的中心为O(x₀，y₀，b)，观测点坐标为(x，y，z)，磁化强度为M；利用泊松方程，计算出磁场各分量的表达式：

其中，M_x、M_y、M_z分别表示磁化强度的x、y和z方向的分量，其中M_x＝McosIcosA，M_y＝McosIsinA，M_z＝MsinI；u₀是磁导率；G是万有引力常数；σ是物体密度；π为圆周率；V_xx，V_yx，V_zx，V_xy，V_yy，V_zy，V_xz，V_yz，V_zz分别是引力位在对应方向上的二阶偏导数；

引力位V由公式(7)计算得到：

其中，v是物体体积，r为物体中心到观测点的距离；

根据磁矩m的定义得m＝Mv，令α＝cosIcosA,β＝cosIsinA,γ＝sinI，将引力位的二阶偏导数代入式(6)消去引力位，直接求得磁场各分量的公式：

以上技术方案中优选的，给初始模型赋值包括6个参数，分别为x、y和z三个位置参数、磁化倾角和偏角两个角度参数以及磁偶极子极距。

以上技术方案中优选的，x和y的初值采用步骤S3中化极插值后数据极大值的位置，模型深度z的初值根据实际情况赋值。

以上技术方案中优选的，步骤S5中，根据式(9)计算吻合度：

其中，mean代表的是求平均值。

以上技术方案中优选的，步骤S5中，吻合度大于0.9即满足要求。

应用本发明的技术方案，具有以下有益效果：

本发明提出了一种自反馈类正则化反演方法，提出了一种基于参数本身变化的自反馈技术，避免了反演时参数更新受数据噪声影响大而变得不稳定的问题。这种新的反演方法在反演的迭代过程中添加了类正则化项，这个类正则化项的权重基于参数本身变化进行自反馈变化，即根据参数本身变化量来确定正则化因子，从而可以得出更加精确的反演结果。与传统方法相比，本发明的新方法可以更有效避免观测数据中噪声的影响，同时使反演过程更加稳健，从而可以更容易跳出局部极值，得出的反演结果更加接近实际位置。

本发明将传统的三个位置参数(即x，y和z)反演变成了新的六参数反演，这样做能更好地拟合实际数据，把剩余磁化强度也包含在反演中；改成6参数反演明显要更符合实际情况，可以得到更加准确的反演结果。同时，本发明采用了人工圈定异常以及通过化极来确定磁偶极子异常体的初始中心位置，让反演的初始模型接近真实位置，从而让反演更快收敛并更准确。

另外，本发明通过对实测数据进行插值，有利于提升反演中的数据信息量，从而提升反演精度。此外，本发明添加了吻合度计算，可以通过吻合度知道本次反演的可靠程度，从而有利于判断本次反演的位置是否准确。

除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。

附图说明

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1是本发明自反馈类正则化反演方法的流程图；

图2是磁偶极子产生磁场示意图；

图3a是化极前原始磁异常数据示意图(白色十字代表异常体真实水平中心位置)；

图3b是化极后磁异常数据示意图(白框代表人工圈定的区域)：

图4是本发明反演算法与国际最新的稳健高斯牛顿法反演位置参数x对比；

图5是本发明反演算法与国际最新的稳健高斯牛顿法反演位置参数y对比；

图6是本发明反演算法与国际最新的稳健高斯牛顿法反演位置参数z对比；

图7是本发明反演算法与国际最新的稳健高斯牛顿法反演拟合差对比；

图8a是异常体在(2m，5m)位置斜磁化添加高斯白噪声合成数据结果；

图8b是异常体在(5m，5m)位置斜磁化添加高斯白噪声合成数据结果；

图8c是异常体在(5m，2m)位置斜磁化添加高斯白噪声合成数据结果；

图8d是异常体在(2m，2m)位置斜磁化添加高斯白噪声合成数据结果；

图9a是异常体在(4.0m，3.01m)位置斜磁化实测数据结果；

图9b是异常体在(3.56m，3.5m)位置斜磁化实测数据结果；

图9c是异常体在(3.96m，3.75m)位置斜磁化实测数据结果；

图9d是异常体在(2.56m，3.7m)位置并伴有另外一个偶极子干扰的斜磁化实测数据结果。

具体实施方式

为了便于理解本发明，下面将对本发明进行更全面的描述，并给出了本发明的较佳实施例。但是，本发明可以以许多不同的形式来实现，并不限于本文所描述的实施例。相反地，提供这些实施例的目的是使对本发明的公开内容的理解更加透彻全面。

除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是旨在于限制本发明。

实施例1：

参见图1，一种自反馈类正则化反演方法，其特征在于，包括以下步骤：

如图3a和图3b所示，步骤S1：将采集的平面磁异常数据进行化极并成图；化极可以消除斜磁化影响，简化磁场形态。因此，化极后的磁异常图相当于是地下的磁异常体投影，因此可以很容易地实现后续的人工圈定磁异常。

优选本实施例中采用的是频率域化极技术，化极后的数据不仅形态变得更加简单，同时数据的极值基本覆盖在地下异常体真实位置，因此按照化极后的磁异常图像可以很容易地进行人工圈定数据进行反演。

步骤S2：根据化极后的磁异常平面图进行人工圈定，并对圈定数据进行插值；本实施例中优选采用克里金插值方法，通过插值增加数据数量，并将插值后的数据作为后续反演的数据；

步骤S3：在化极并插值后的数据中寻找极大值的平面位置，并将此平面位置作为反演(即类正则化反演)的初始位置；

步骤S4：开始反演迭代过程，首先给初始模型赋值，然后对模型参数进行迭代更新，同时计算拟合差；迭代完毕后，得到输出反演的结果；

优选的，步骤S4中，给初始模型赋值包括6个参数，分别为x、y和z三个位置参数、磁化倾角和偏角两个角度参数以及磁偶极子极距。其中，x和y的初值采用步骤S3中化极插值后数据极大值的位置，模型深度z的初值根据实际情况赋值，优选模型深度z的初值赋值为1m，因为在磁偶极子的实际应用中，一般异常体埋深通常在0.5到2米，故初始值设为1米方便迭代收敛。

步骤S5：反演迭代结束，得到反演模型参数以及吻合度；如果吻合度满足要求，则输出结果；如果吻合度不满足，则返回步骤S2，重新人工圈定数据，直到输出吻合度满足要求为止。优选的，吻合度大于0.9即满足要求。

以上为本实施例中自反馈类正则化反演方法处理流程，接着继续对本实施例中步骤S4中的反演(即类正则化反演)过程进行详细说明。

本实施例中所述步骤S4中，采用式(1)和式(2)对模型参数进行迭代更新：

其中，m^(k)、m^(k+1)分别为第k次和k+1次迭代的模型；J_k为第k次迭代的雅各比矩阵，即正演模型数据对反演模型的偏导数矩阵；r_k为第k次迭代时的拟合差；

表示J_k的转置矩阵；I为单位矩阵；λ_k为第k次迭代的调节因子，为了与传统的超参数化反演对应称之为正则化因子；由式(2)可知，调节因子为上一次模型更新量的绝对值之和。

本实施例在偶极子反演中，选取的反演目标函数如式(3)：

其中，m为偶极子需要反演的模型；Φ_d(m)为正演模型数据与观测数据之间的拟合差的二范数，即：

其中，F[m]为模型正演数据矢量；d_i ^pred为模型正演数据矢量的第i个数据；d^obs为观测数据矢量；

为观测数据矢量的第i个数据；W_d代表的是观测误差的对角矩阵，

是观测误差的对角矩阵的第i行第i列的数据，N表示观测数据的总量；

在反演过程中，要让观测数据与正演模型数据之间的误差，即拟合差达到最小；这个过程需要通过式(1)和式(2)对模型参数不断修改迭代，直至完成设定的迭代次数；迭代次数本领域人员可以根据实际情况进行设定，满足当前后两次反演得到m的值基本不变即可，此时可以认为反演已经收敛。

进一步地，本领域人员知晓，在磁偶极子反演问题中，首先要解决的是磁偶极子的正演问题，对磁偶极子产生的磁场进行正确的模拟才能保证反演的精度，因此接下来将对正演进行详细说明。

本实施例中正演模型数据通过偶极子正演得到，具体是通过式(5)得到磁总场：

ΔT＝H_axL₀+H_ayM₀+Z_aN₀ (5)，

其中，L₀＝cosI₀cosA₀，M₀＝cosI₀sinA₀，N₀＝sinI₀；I₀为地磁场的磁偏角，A₀为地磁场的磁倾角，ΔT为磁总场，H_ax、H_ay和Z_a分别为磁场的x、y和z分量上的大小。

对于一个在地球的带有磁性的物体来说，如果这个物体的形状是球状或者是类椭球状的(如图2所示)，那么这个物体产生的磁场就可以等价为磁偶极子产生的磁场。磁偶极子的磁场可以由泊松公式推导而来，具体地，式(5)中H_ax、H_ay和Z_a可以通过如下计算得到：

如图2所示，x和y轴的坐标系统是位于北部和东部，z轴是垂直向下的；设总磁化的磁偏角为I、磁倾角为A，磁偶极子R的中心为O(x₀，y₀，b)，观测点P坐标为(x，y，z)，磁化强度为M；利用泊松方程，计算出磁场各分量的表达式：

其中，M_x、M_y、M_z分别表示磁化强度的x、y和z方向的分量，其中M_x＝McosIcosA，M_y＝McosIsinA，M_z＝MsinI；u₀是磁导率；G是万有引力常数；σ是物体密度；π为圆周率；V_xx，V_yx，V_zx，V_xy，V_yy，V_zy，V_xz，V_yz，V_zz分别是引力位V在对应方向上的二阶偏导数，如：V_xy表示先对V求x方向偏导再求y方向偏导，

引力位V由公式(7)计算得到：

其中，v是物体体积，r为物体中心到观测点的距离；

根据磁矩m的定义可知m＝Mv，令α＝cosIcosA,β＝cosIsinA,γ＝sinI，将引力位的二阶偏导数代入式(6)消去引力位，可以直接求得磁场各分量的公式：

本实施例中优选的，所述步骤S5中，根据式(9)计算吻合度：

其中，mean代表的是求平均值。

更进一步的，为了测试本实施例方法的有效性，我们通过在理论数据和某地异常体实测数据进行验证。

在理论数据上，通过对数据加高斯白噪声，让数据更加接近真实情况。同时本实施例方法与最新发表的稳健的高斯牛顿法反演进行了对比，如图4，图5，图6和图7所示。在对比中可以发现，本实施例的方法很快就收敛了，并收敛到了真实位置参数。而国际最新发表的稳健的高斯牛顿法反演收敛较慢，同时也无法收敛到真实位置参数。此外，虽然国际最新的稳健高斯牛顿法可以收敛到更低的拟合差，但是是陷入了局部最小值，对数据过拟合，把噪声也进行了拟合，搜索能力较差。本实施例的算法很好地收敛到了真实位置，不会对数据进行过拟合，从而把噪声也拟合了。

如图8a-图8d为异常体的模拟数据，在斜磁化下产生的磁异常，并添加高斯白噪声。从结果可以看到，在理论合成数据中，无论异常体在哪个方位，本实施例的新算法都可以很准确得反演出其真实位置。

此外，本实施例还通过实测异常体数据进行了测试。在实测数据中图9a，图9b，图9c三组实测数据没有设置干扰体，而图9d实测数据人为添加了一个磁偶极子干扰大大增加了反演难度。但从结果来看，在实测数据上本实施例的新反演算法和流程还是能很准确地反演出异常体的真实位置。可见，本实施例的新方法可以十分有效地在各种实际情况下反演出偶极子的真实位置。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种自反馈类正则化反演方法，其特征在于，包括以下步骤：

步骤S1：将采集的平面磁异常数据进行化极并成图；

步骤S2：根据化极后的磁异常平面图进行圈定，并对圈定数据进行插值；

步骤S3：在化极并插值后的数据中寻找极大值的平面位置，并将此平面位置作为反演的初始位置；

步骤S4：开始反演迭代过程，首先给初始模型赋值，然后对模型参数进行迭代更新，同时计算拟合差；迭代完毕后，得到输出反演的结果；

步骤S5：反演迭代结束，得到反演模型参数以及吻合度；如果吻合度满足要求，则输出结果；如果吻合度不满足，则返回步骤S2，重新圈定数据，直到输出吻合度满足要求为止；

所述步骤S4中，采用式(1)和式(2)对模型参数进行迭代更新：

其中，m^(k)、m^(k+1)分别为第k次和k+1次迭代的模型；J_k为第k次迭代的雅各比矩阵，即正演模型数据对反演模型的偏导数矩阵；r_k为第k次迭代时的拟合差；

表示J_k的转置矩阵；I为单位矩阵；λ_k为第k次迭代的调节因子，为了与传统的超参数化反演对应称之为正则化因子。

2.根据权利要求1所述的自反馈类正则化反演方法，其特征在于，在偶极子反演中，选取的反演目标函数如式(3)：

其中，m为偶极子需要反演的模型；Φ_d(m)为正演模型数据与观测数据之间的拟合差的二范数，即：

其中，F[m]为模型正演数据矢量；

为模型正演数据矢量的第i个数据；d^obs为观测数据矢量；

为观测数据矢量的第i个数据；W_d代表的是观测误差的对角矩阵，

是观测误差的对角矩阵的第i行第i列的数据，N表示观测数据的总量；

根据式(1)和式(2)对模型参数不断修改迭代，直至完成设定的迭代次数。

3.根据权利要求2所述的自反馈类正则化反演方法，其特征在于，正演模型数据通过偶极子正演得到，具体是通过式(5)得到磁总场：

ΔT＝H_axL₀+H_ayM₀+Z_aN₀ (5)，

其中，L₀＝cosI₀cosA₀，M₀＝cosI₀sinA₀，N₀＝sinI₀；I₀为地磁场的磁偏角，A₀为地磁场的磁倾角，ΔT为磁总场，H_ax、H_ay和Z_a分别为磁场的x、y和z分量上的大小。

4.根据权利要求3所述的自反馈类正则化反演方法，其特征在于，其中H_ax、H_ay和Z_a的计算方式如下：

设总磁化的磁偏角为I、磁倾角为A，磁偶极子的中心为O(x₀，y₀，b)，观测点坐标为(x，y，z)，磁化强度为M；利用泊松方程，计算出磁场各分量的表达式：

其中，M_x、M_y、M_z分别表示磁化强度的x、y和z方向的分量，其中M_x＝McosIcosA，M_y＝McosIsinA，M_z＝MsinI；u₀是磁导率；G是万有引力常数；σ是物体密度；π为圆周率；V_xx，V_yx，V_zx，V_xy，V_yy，V_zy，V_xz，V_yz，V_zz分别是引力位在对应方向上的二阶偏导数；

引力位V由公式(7)计算得到：

其中，v是物体体积，r为物体中心到观测点的距离；

根据磁矩m的定义得m＝Mv，令α＝cosIcosA,β＝cosIsinA,γ＝sinI，将引力位的二阶偏导数代入式(6)消去引力位，直接求得磁场各分量的公式：

5.根据权利要求1所述的自反馈类正则化反演方法，其特征在于，步骤S4中，给初始模型赋值包括6个参数，分别为x、y和z三个位置参数、磁化倾角和偏角两个角度参数以及磁偶极子极距。

6.根据权利要求5所述的自反馈类正则化反演方法，其特征在于，x和y的初值采用步骤S3中化极插值后数据极大值的位置，模型深度z的初值根据实际情况赋值。

7.根据权利要求2所述的自反馈类正则化反演方法，其特征在于，步骤S5中，根据式(9)计算吻合度：

其中，mean代表的是求平均值。

8.根据权利要求7所述的自反馈类正则化反演方法，其特征在于，步骤S5中，吻合度大于0.9即满足要求。

Images