具体实施方式

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

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

实施例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实测数据人为添加了一个磁偶极子干扰大大增加了反演难度。但从结果来看，在实测数据上本实施例的新反演算法和流程还是能很准确地反演出异常体的真实位置。可见，本实施例的新方法可以十分有效地在各种实际情况下反演出偶极子的真实位置。

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