斜磁化条件下多磁源目标的边界识别方法

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

(1.陆军工程大学石家庄校区 车辆与电气工程系， 河北 石家庄 050003；2.中国空气动力研究与发展中心 高速空气动力研究所， 四川 绵阳 621000)

0 引言

边界识别是解释磁测数据的重要过程，其主要目的是通过获得的磁测数据对地下或水下的未知磁源的物理边缘进行可视化，在军事侦察(未爆弹、潜艇和水雷等)、水文及工程地质勘探等领域有着广泛的应用[1-5]。

传统的边界识别方法主要包括总水平导数法、倾斜角法以及Theta 图(TM)法等[6-9]，但是这些传统的边界识别方法存在一定局限性：总水平导数法利用磁场沿z轴方向的分量数据Bz在水平方向上梯度的模值对水平边界进行计算，该方法不能同时对强弱磁异常进行显示；倾斜角法计算Bz分量沿z轴方向空间导数获得的Bzz与总水平导数法模值的反正切，获得的弧度值为磁性目标的二维边界结果，该方法能够有效地平衡不同强弱磁异常的振幅，但是磁源的边界没有非常明确的定义；TM法计算总水平导数法获得的模值与Bz分量在3个方向上梯度模值比值的反余弦获得弧度值，利用弧度值结果能够同时描绘不同磁性强弱的磁源，但是计算精度难以满足对磁源边界进行清晰描述的要求。在倾斜磁化条件下，利用传统边界识别方法获得的磁源边界将会偏离磁源的真实位置[10-12]。此外，在噪声影响下，磁源的识别精度会受到一定影响。

尽管目前已经开始将磁梯度张量信息引入到磁源的二维边界反演计算中，但是多局限于对垂直磁化条件下的磁源进行计算[13]，而基于磁梯度张量斜磁化条件下的多磁源目标二维边界反演技术仍然处于起步阶段。且在多磁源目标边界计算过程中，测量得到的磁分量信号受环境磁场、系统误差以及人为操作不当等因素会对磁测数据引入噪声，严重影响计算结果。

因此，针对多磁源目标二维边界反演中存在的问题，本文提出基于磁梯度张量的多磁源边界识别方法。该方法首先基于改进K奇异值分解(K-SVD)的方法对磁测信号进行降噪处理，以实现对磁测数据的预处理，然后利用互相关理论对不同计算区域内的单目标磁源磁化方向进行估计，以实现将磁测数据转化为垂直磁化条件下的数据，最后计算磁梯度张量矩阵特征值的绝对值之和与垂直磁化条件下的磁梯度张量Bzz分量的比值，获得磁源的水平边界。

1 磁梯度张量理论

在笛卡尔坐标系下，磁梯度张量数据是磁场B=[Bx,By,Bz]T沿着x轴、y轴、z轴3个方向的空间导数,Bx、By和Bz为磁场三分量数据。假设x轴指向地理北向，y轴指向地理东向，z轴指向垂直向下的方向，在观测平面(x,y,z)的磁梯度张量数据可表示为

(1)

式中：Γ为磁梯度张量矩阵；Bij(i,j=x,y,z)为磁梯度张量分量，磁梯度张量矩阵为对称矩阵，其矩阵的迹为0. 因此，磁梯度张量具有Bxx、Bxy、Bxz、Byy、Byz这5个独立分量，其他分量都可以从这5个独立分量推导出来。磁梯度张量矩阵的特征值分解可以表示为

(2)

式中：λ1、λ2、λ3为磁梯度张量矩阵的特征值，按照单调递减的顺序依次排列为|λ1|≥|λ2|≥|λ3|≥0，v1、v2、v3分别为特征值λ1、λ2、λ3对应的特征向量.

为了获得磁梯度张量数据，设计了一个由4个三轴磁通门组成的十字形磁梯度张量系统，系统如图1所示。两个磁力计在同一方向上的基线距离为2d，d为传感器测量轴到磁梯度张量系统中心之间的距离。十字形磁梯度张量系统的测量数据可以定义为

(3)

式中：Bij(i=1,2,3,4,j=x,y,z)为第i个传感器j轴的分量数据。

图1 十字形磁梯度张量系统Fig.1 Cross magnetic gradient tensor system

该测量系统的结构误差[14]可以表示为

(4)

根据(4)式可以看出，构造的磁梯度张量系统测量值与理论值间的误差较小，可以在实际测量中忽略不计。

2 多磁源二维边界反演方法

2.1 基于改进K-SVD的磁测数据预处理方法

在对磁源的探测与识别过程中，常常需要考虑噪声的影响，因此在计算过程中要降低磁场数据中的噪声。在降噪过程中，为了保证不同方向上磁场分量之间的相互关系，将传感器上不同的磁场分量利用矢量合成原理转化为总幅值异常以及方向余弦进行降噪处理：

(5)

式中：B为磁总场幅值异常；l、m、n分别为沿坐标轴3个方向的方向余弦。

K-SVD降噪模型[15]可以表示为

(6)

式中：αe为最终的稀疏向量解；y是测试数据，y=b+q，b是理论数据，q是均值为0的高斯白噪声；α∈Rs为K-SVD字典矩阵的稀疏向量；D∈Ra×s(s≥a)为字典矩阵，其中s为字典大小，a为测试数据个数；T为容差极限。在实际应用中，噪声为标准高斯白噪声的情况很少，因此，K-SVD降噪方法在应用中存在一定的局限性。因此，本文提出了一种改进K-SVD(IK-SVD)降噪模型来重建磁数据。

(7)

式中：yi是测试数据中的第i列数据；Mi是观测数据去均值化后的第i列数据。

(8)

利用IK-SVD降噪方法对磁总场幅值异常B及方向余弦(l,m,n)进行降噪处理，利用(5)式对信号进行逆向重构[16]：

(9)

式中：B′x、B′y、B′z分别为降噪后的磁测数据三分量；B′为降噪后的磁总场幅值异常；(l′,m′,n′)为降噪后的方向余弦。

2.2 磁梯度张量数据转化理论

在斜磁化条件下，对磁性目标进行边界识别计算时会受到磁化方向的影响而对计算结果产生干扰。为了解决这一问题，本文利用基于归一化磁源强度(NSS)和化极(RTP)互相关的方法对磁化方向进行估计，进而利用获得的磁化方向对降噪后的磁测数据进行RTP计算，最终获得垂直磁化条件下磁源的磁梯度张量数据。NSS是由磁梯度张量矩阵推导出来的一个旋转不变量[17]，可以表示为

(10)

定义实测区域内的磁异常RTP结果ΔTr与NSS数值μ的互相关系数C[18]为

(11)

磁总场异常TMI RTP公式[18]为

(12)

磁场三分量的傅里叶变换[19]可以表示为

(13)

2.3 磁源边界识别理论

在对多目标磁源进行边界计算过程中，利用磁梯度张量矩阵获得的特征值的和(SAE)对磁源的分布区域进行划分,这是因为SAE数据的主正极值与磁源的实际分布位置相一致，受磁化方向影响较小。同时，为了获得磁源的边界，利用特征值边缘检测的方法来进行计算，边界计算结果Es的具体计算公式如下：

(14)

在垂直磁化条件下，磁梯度张量Bzz分量的正负数据分布与磁源的边界一致，因此，磁源的边界位于Es的正负数据交界处，即磁源的边界可以用Es的零轮廓线来进行描述。

3 仿真分析

为了证明本文方法的有效性，在测量区域内放置两个斜磁化条件下的模型进行仿真分析，如图2所示。地磁背景场的磁化方向为Ib=60°、Db=-20°. 模型1为多边形薄板，磁矩为50 A/m，磁化方向为Ia1=60°、Da1=-20°，厚度为0.1 m，中心深度为0.25 m. 模型2为长方形薄板，磁矩为50 A/m，磁化方向为Ia2=45°、Da2=55°，厚度为0.1 m,中心深度为0.35 m.

图2 组合模型Fig.2 Combined model

图3 观测面上磁测数据Fig.3 Magnetic measurement data on the observation surface

分别向4个传感器每个轴上加入信噪比为30%的高斯白噪声，模拟环境因素带来的误差，并向观测数据中随机加入0.5×rand(max(B))的冲击噪声，模拟实际测量中由于操作误差所带入的冲击误差。观测面上获得的磁测数据如图3所示，图3(h)中的白线为模型的理论水平位置，图3(i)中的白线为划分的计算区域。为了证明本文降噪方法的有效性，分别利用形态学滤波器、低通滤波器、小波变换方法、EMD分解方法以及K-SVD方法与本文的IK-SVD方法进行对比，其中形态学滤波器通过构造开- 闭和闭- 开组合的组合滤波器并利用三角形结构元素进行滤波运算，低通滤波器采用巴特沃斯滤波器确定截止频率，EMD滤波器与小波滤波器采用软阈值的滤波方法[16]，小波方法中的小波基为coif5，分解层数为5层。K-SVD计算方法中，块大小为4，训练字典大小为512. 本文方法的窗口大小为4,移动步长为1. 分别利用信噪比SNR以及均方根误差RMSE作为评判降噪能力的指标，其中

表1 不同计算方法的降噪结果Tab.1 Denoised results of different calculation methods nT/m

图4 去噪声后观测面上磁测数据Fig.4 Denoised magnetic measurement data on the observation surface

利用IK-SVD方法降噪后的磁测数据如图4所示。利用(11)式估计两个计算窗口内模型的磁化方向，其中模型1的磁化方向估计值为Is1=56°、Ds1=-18°，模型2的磁化方向估计值为Is2=48°、Ds2=59°. 图5为磁梯度张量转化前后不同计算区域内的Bzz数据。利用(13)式获得垂直条件下的磁梯度张量Bzz分量，如图5(c)和图5(d)所示，白色实线为磁源的理论水平分布区域。与图5(a)和图5(b)对比可以看出：在垂直磁化条件下，实线的位置与Bzz主正极值的分布更加一致；而在斜磁化条件下，Bzz主正极值与实际轮廓位置存在一定的偏移。因此，利用本文的相关性来估计磁化方向的方法可以有效地降低斜磁化的影响。

图6 不同方法获得的边界识别结果Fig.6 Edge detection results obtained by different methods

图7 降噪及磁梯度张量转化前后的Es轮廓Fig.7 Es values before and after noise reduction and magnetic gradient tensor conversion

4 试验验证

为了验证上述方法的有效性，进行实测试验。测区位于石家庄某地，探测仪器及磁源位置分布如图8所示。试验所用探头是利用英国Bartington公司生产的Mag-03三轴磁通门传感器搭建的十字形磁梯度张量探头，将传感器固定在塑料支架上，两个磁传感器在同一方向上的基线距离为0.4 m，在试验过程中，探头固定在无磁试验台架上，利用扫描的方式对待测区域的每一点进行测量，利用非线性最小二乘方法对传感器安装误差及软磁干扰进行校正及补偿[20]。测区的背景磁化方向为Ib=55°、Db=-16°，测点的间距为0.1 m，在测区内对两个长方形铁块的磁测数据进行测量，获得磁测数据如图9所示。

图8 试验测试仪器及磁源分布Fig.8 Experimental equipment and magnetic source distribution

图9(h)中铁块1尺寸为沿南北方向长度为0.36 m，沿东西方向宽度为0.23 m，中心位于(1.15 m，0.65 m)，顶部深度为0.2 m，厚度为0.1 m. 图9(h)中铁块2尺寸为沿南北方向长度为0.36 m，沿东西方向宽度为0.234 m，中心为(0.45 m，1.5 m)。顶部深度为0.15 m，厚度为0.05 m. 分别利用形态学滤波器、低通滤波器、小波变换方法、EMD分解方法以及K-SVD方法与本文的IK-SVD方法进行对比，计算结果如图10所示。由图10可以看出，通过对获得的磁梯度张量5个独立分量进行对比，可以看出本文降噪方法具有最优的降噪能力。将SAE图划分为覆盖各个铁块的计算区域，如图9(i)中的实线所示。

利用IK-SVD方法降噪后的磁测数据如图11所示。利用(11)式估计两个计算窗口内铁块的磁化方向，其中铁块1的磁化方向估计值为Ie1=12°、De1=-5°，铁块2的磁化方向估计值为Ie2=3°、De2=-13°. 图12为磁梯度张量转化前后不同计算区域内的Bzz数据。利用(13)式获得垂直条件下的磁梯度张量Bzz分量，如图12(c)和图12(d)所示，白色实线为磁源的理论水平分布区域。与图12(a)和图12(b)对比可以看出：在垂直磁化条件下，实线的位置与Bzz主正极值的分布更加一致；而在斜磁化条件下，Bzz主正极值与实际轮廓位置存在一定的偏移。图13为具有噪声的斜磁化条件下不同计算方法识别的边界轮廓。图14为及磁梯度张量转化前后的Es轮廓。由图13和图14可以看出，本文方法能够有效提高多磁源目标在噪声及斜磁化条件下的边界识别能力，进一步证明了本文方法的有效性。

图9 观测面上磁测数据Fig.9 Magnetic measurement data on the observation surface

图10 不同计算方法的降噪结果Fig.10 Denoised results of different calculation methods

图11 去噪声后观测面上磁测数据Fig.11 Denoised magnetic measurement data on the observation surface

图12 磁梯度张量转化前后不同计算区域内的Bzz数据Fig.12 Bzz data in different computing areas before and after magnetic gradient tensor conversion

图13 不同方法获得的边界识别结果Fig.13 Edge detection results obtained by different methods

图14 降噪及磁梯度张量转化前后的Es轮廓Fig.14 Es values before and after noise reduction and magnetic gradient tensor conversion

5 结论

本文针对磁源在二维边界反演过程中存在的问题，提出了斜磁化条件下多磁源目标二维边界识别方法，该方法能有效地在噪声以及斜磁化条件下对多磁源目标的二维边界进行识别。通过仿真和试验证明了本文方法具有如下优势：

1)提出了利用IK-SVD方法对磁测数据总幅值异常以及方向余弦进行预处理，能够在预处理过程中能够有效保证磁测数据三分量之间的内在联系，且具有较强的抗噪性。

2)利用NSS与RTP互相关的方法估计不同计算区域内磁源的磁化方向，将磁测数据转化为垂直磁化条件下的磁测数据，能够有效消除斜磁化对计算结果带来的影响。

3) 通过计算SAE与垂直条件下的磁梯度张量Bzz分量的比值零线位置，能够准确对磁源的水平边界进行求解。

本文方法适用的探测范围可以从米到千米，具有一定的国防和工业应用背景，未来的工作旨在检测剩磁以提高其准确性，该成果可应用于未爆弹探测、矿产勘探等领域。