磁力勘探是利用地表的磁测数据来评估地下未知磁异常体的轮廓形态^[1]。在传统反演方法中，常假设总磁化方向与感应磁化方向是一致的。然而，大多数情况下，由于剩余磁化的存在，导致总的磁化方向产生明显的改变，与感应磁化的方向并不一致^[2]。同时，由于地下待求未知参数远远多于观测数据点数^[3]，会导致线性方程求解的模糊性，难以得到方程最优解。为了降低测量误差，常采用L₂或者Tikhonov正则化思想^[4-5]。确定正则化参数的主要方法有L(LC)曲线、广义交叉验证(GCV)^[6-8]，其中LC曲线法计算量较大，需要计算多个正则化参数来确定最优解; GCV方法难以保证在每次迭代过程中都能获得GCV功能矩阵的极小值。

针对上述问题，本文提出基于归一化磁源强度的自适应正则化聚焦反演方法。首先，利用归一化磁源强度数据进行反演，减弱剩余磁化对反演结果的影响；然后，引入深度加权矩阵和最小支撑矩阵，避免由于反演参数多于采集点数而导致的欠定问题，并有效解决核函数随深度增大而快速衰减的问题；最后，根据Morozov偏差原则自适应地确定目标函数在迭代过程中的正则化参数，提高了迭代速度和求解精度。

1 Tikhonov正则化聚焦原理

磁测数据正演公式如下：

(1)

式中，核函数矩阵A的元素A_ij为第i(i=1, 2, 3, …, p)个观测点第j(1, 2, 3, …, n)个单元的响应, 向量d ∈为观测数据矩阵，向量m ∈为待测模型参数矩阵。利用Tikhonov正则化公式将目标函数矩阵转化为最小二乘公式：

(2)

式中，W_d=diag(1/η₁, …, 1/η_p)为数据加权矩阵，η_i为第i个基准噪声的标准偏差，D为正则化矩阵，m_apr为模型参数先验信息，α为正则化参数。将式(2)转为如下形式：

(3)

(4)

式中，m_apr)。

正则化矩阵可以表示为D^(k)= W_e^(k)W_depth。Li等^[6]给出深度加权矩阵W_depth=diag(1/(z_j+ξ)^β)，其中z_j为地下第j个剖分单元的中心埋深，ξ>0是为避免出现奇点而引入的数值，β为可以控制加权函数大小的正数。Vatankhah给出最小支撑矩阵W_e^(k)=diag((m^(k)－m^(k－1))+ε²)^－1/2, k>0，W⁽⁰⁾= I, m⁽⁰⁾=m_apr，ε>0为聚焦参数^[8]。

2 基于奇异值分解的归一化磁源强度自适应反演方法 2.1 迭代奇异值分解

设置迭代过程为：

(5)

(6)

在求解J(α)时，对矩阵进行奇异值分解(SVD)，并在迭代过程中设置物性参数上下限约束函数[m_min, m_max]。矩阵Ã可以表示为Ã=UΣV^T，式中，U、V分别为p×p和n×n正交矩阵，Σ=diag([δ₁, δ₂, …δ_p])，其中奇异值δ₁≥δ₂≥…≥δ_p>0。为了防止高频观测误差被放大，引入正则化滤波因子。因此，对式(5)进行奇异值分解得：

(7)

式中，u_i和v_i分别为矩阵U和V的第i列。

2.2 归一化磁源强度数据计算

磁梯度张量可以由一个二阶矩阵表示：

(8)

式中，G为磁梯度张量，B_i(i=x，y，z)为磁场三分量，B_αβ(α、β=x, y, z)为磁梯度张量的9个分量。

归一化磁源强度是由磁偶极子磁梯度张量矩阵推导的旋转不变量，分布于场源中心。假定磁偶极子在观测面(x，y，z)上引起的磁梯度张量矩阵的特征值按单调递减顺序为λ₁≥λ₂≥λ₃，则归一化磁源强度u(x、y、z)可以表示为：

(9)

式中，C_m=10^-7H/m，m为磁偶极子磁矩，r为异常源到观测点的距离。由式(9)可见，归一化磁源强度与距离的四次方成反比，与磁矩成正比。

2.3 基于Morozov偏差原则的自适应正则化参数选择

徐新禹等^[9]讨论了在Tikhonov正则化方法中如何有效选择正则化参数。Zhdanov等^[10]提出一种自适应正则化参数的选择方法，定义α^(k)=α⁽¹⁾q^(k)，q为衰减因子。但是这种方法不能保证每次迭代得到的正则化参数为最佳解，影响求解精度。本文提出基于Morozov偏差原则(MDP)的正则化参数选择方法。

Mead^[11]指出，MDP是基于后验策略的偏差原理，利用先验误差信息，使得实测数据与预测数据的残差为最小值的正则化参数即为最优正则化参数。MDP遵循自由度为p-n的卡方分布：

(10)

利用奇异值分解法求解最终的正则化参数：

(11)

3 仿真分析

测试模型见图 1，黑色表示待测异常体，其中格点代表观测点矩阵，各部分参数见表 1，其中地磁场倾角为65°，地磁场偏角为－25°。将地下待测空间划分为22×22×10=4 840个单元格，每个单元格均为边长为0.1 m的正方体，地面观测点为22×22=484个网格。反演时，在模型理论归一化磁源强度数据中加入0.03×std(d, 1)×rand(484, 1)的高斯噪声作为实测数据，理论模型正演结果见图 2。可以看出，归一化磁源强度数据与异常体理论模型位置基本重合，进一步论证了归一化磁源强度数据是弱敏感于磁化方向的。反演结果见图 3，可以看出，该方法能够准确还原不同位置、不同形状的模型的轮廓形态，具有较好的横向分辨率。从图 3(b)可以看出，纵向分辨率有所下降，存在一定的边缘效应，这是由于磁总场模量数据仅包含异常体的幅值信息，不包含相位信息。

4 实验验证

测试装置为十字形传感器阵列结构，利用Bartington公司的三轴磁通门传感器搭建的磁梯度张量探头，该传感器的分辨率为0.1 nT，实测仪器包括十字形探头、数字采集模块以及软件操作终端。实验中，将探头固定在无磁台架上，在台架的横向和纵向滑道上标注刻度值，通过扫描对实测区域的每一点磁梯度张量值进行测量(图 4)。测区位于石家庄某地，在地理坐标系下，对水平圆柱体铁筒和长方体铁块的组合磁异常进行测试。

长方体铁块位于测区东南方向，中心坐标为(1.5, 0.4, 0.3)，尺寸为0.23 m×0.36 m×0.2 m。铁筒位于测区西北方向，中心坐标为(0.6, 1.2, 0.2)，圆筒直径为0.1 m，长0.46 m。

实测数据见图 5。利用归一化磁源强度数据进行反演，获得的三维反演结果见图 6。由反演结果可以看出，本文方法均能很好地重构地下磁性体的轮廓形态，而且反演更加灵活，能够对任意数量的磁性体进行反演。其缺点是在反演过程中损失了相位信息。

5 结语

1) 利用归一化磁源强度数据对磁性目标进行反演，反演结果不受剩余磁化方向的影响，能够很好地校正剩磁效应带来的偏差。

2) 引入深度加权矩阵和最小支撑矩阵对反演模型进行约束，避免了由于反演参数多于采集点数而导致的多解性问题，有效解决了核函数随深度增大而快速衰减的问题。

3) 在迭代过程中，通过奇异值分解法对核函数矩阵进行奇异值分解，利用物性参数上下限函数将反演解转换到合理值域；采用Morozov偏差原则对迭代过程中的正则化参数进行选择，实现了正则化参数的自适应调整。

虽然本文方法能够重构异常体的轮廓形态，但是归一化磁源强度数据仅具有幅值信息，不包含相位信息，反演结果存在一定的边缘效应。