AFCM聚类和张量不变量用于磁源多目标定位

李青竹，李志宁，石志勇，范红波

（陆军工程大学 车辆与电气工程系，河北 石家庄 050003）

1 引言

多个磁目标同步探测是极具挑战的工作。作为磁法勘探领域的研究热点，目前多目标磁定位研究已在山脉矿产勘探［1］、脑磁图定位［2］、航天器 磁 特 征 检 测［3］和 未 爆 弹（Unexploded Ordnance，UXO）探测［4］等领域中具有重要应用。一般地，目标定位和识别（Positioning and Recognition，PR）即是利用测量磁场数据对区域内磁性目标的位置和磁矩进行估算的问题。当磁化均匀且源中心与观测点间距大于目标最大尺寸的2.5倍时，磁源通常被认为是磁偶极子［5］。磁偶极子物性参数可用六个自由度进行描述，即空间位置（x，y，z）和磁矩（mx，my，mz），故多个紧凑源PR检测即对多个磁偶极子的位置和磁矩共六个参数进行估计。磁梯度张量（Magnetic Gradient Tensor，MGT）［6］是磁场在空间正交方向上的变化率，也是磁标势的二次偏微分，其信息度更高，受磁化方向影响小，已成为磁探测领域研究前沿。此外，MGT特征值和不变量等信息能减少系统旋转误差和背景磁场干扰，对于磁源PR更具优势。目前基于MGT测量的磁偶极子单点定位算法已较成熟，如张量衍生不变关系定位法［7］和二阶磁张量欧拉反褶积定位法［8］，通过磁梯度张量系统（Magnetic Gradient Tensor System，MGTS）［5］进行航线测量，可对测区内孤立磁源具有较高的单点定位精度。然而，当多个目标出现在不同埋深，且其分布和磁矩均不相同时，弱磁性目标常因多个磁异常场相互叠加被淹没，这使得仅利用单点方法进行磁源PR探测陷入瓶颈。

针对多磁源PR问题，1993年Hansen等［9］提出利用Werner反褶积方法对多源磁目标进行二维反演；此 后Hansen［10］又 将 多 源Werner反 褶 积方法拓展到了三维并进行了盆地复杂数据测量。这些算法的一个主要缺点是需要在实际计算之前定义磁性目标的数量，可实际应用中目标数目是很难提前知晓的。王林飞等人［11］通过归一化磁场垂向导数与z方向解析信号的比值得到MGT解析信号来识别目标最佳范围和深度，李金朋等人［12］也提出利用归一化磁源强度和同质性方程对测区平面进行边界识别以估计平面下方目标的位置和埋深，但这些方法需要人为对等高线图进行鉴别，以确定目标的水平分布和数目，探测效率低、虚警率高且结果受噪声干扰较大。

由此可见，对测区内目标数目进行自适应识别以适用于更成熟的定位方法并提高位置估计精度是必要的。Ugalde等［13］曾使用核密度分布算法过滤不相关目标位置解，通过模糊c-means聚类（Fuzzy c-means Clustering，FCM）［14］算法的聚类质心来估计异常源的地质分布和走向，这种解决方法可能为磁偶极子的多目标检测和识别提供思路。此外，Davis［15］开发了一种基于希尔伯特变换的扩展欧拉反褶积方法用于多个UXO目标的自动检测，Davis使用不同滑动窗口内的反演结果来估计目标位置。参考划分方窗区域的思路，我们尝试对目标分布区域进行预先识别，以提高多目标PR检测精度。然而，经典FCM算法需要预设初始聚类质心数目，但仅靠经验法预设数目无法保证FCM收敛到最优解，也不能为多目标的自适应PR检测提供帮助。目前已有一些自适应模糊c-means（Adaptive Fuzzy c-means Clustering，AFCM）聚类算法［16-17］，通过评价聚类效果的有效性函数（Clustering Validity Function，CVF）和试错迭代来确定最优聚类数。本文中，我们设计了适用于目标初始位置点云的AFCM聚类有效性函数，为实现多目标自适应PR检测提供可能。

本文方法为基于AFCM聚类和张量不变量的目标多源定位方法。首先，利用MGT不变量中对磁化方向相对不敏感的归一化磁源强度和张量缩并［18］来对MGTS网格内测量节点进行目标预识别和区域划分。接着利用张量衍生不变关系定位方法对识别区域内的网格节点进行单点定位计算，获得一定数量的目标初始位置解集。最后，利用AFCM算法对这些解集形成的稠密点云进行三维聚类并自动估计簇质心数目和质心位置。目标正上方处的MGT矩阵和位置坐标可用以计算磁矩矢量。该方法具备如下优势：（1）分布区域的预识别能极大提高AFCM的聚类自适应识别精度；（2）由于MGT随目标距离四次方衰减［18］，处于源正上方区域的节点处定位结果更加可靠，且该过程能排除大部分虚假解；（3）自适应算法的目标识别过程不需要人为对测区等高线特征进行鉴别和判断。

2 张量不变量用于目标区域识别

2.1 MGT的不变量

无源静磁场中，MGT矩阵可对角化且仅五个张量分量相对独立［5］。以G表 示MGT矩 阵，λ1、λ2和λ3表示G的特征值，v1、v2、v3分别为特征值λ1、λ2、λ3对 应 的 特 征 向 量，且 满 足λ3-I0λ2+I1λ-I2=0 及 |λ1|≥|λ3|、|λ2|≥|λ3|、λ2≥λ3≥λ1，λ3称中间特征值，λ1和λ2分别称最小和最大特征值。若gij（i，j=x，y，z）为MGT分量，则G及其部分MGT不变量可由以下方程表示：

其中：不变量I0和I2分别为G的迹和行列式；u即为归一化磁源强度（Normalized Source Strength，NSS）［18］；c为G的Frobenius范数，又称张量缩并（Contraction of Tensor，CT）。

2.2 不变量改进倾斜角的目标区域识别

平面网格测量的优势是测区内观测点多，相比单点或航线测量信息量更大。网格测量可利用磁场特征量计算，获取测区下方目标的水平分布信息，并极大提高多磁性目标定位的效率和精度。由于待测物体间磁异常场相互叠加，弱磁物体常被强磁物体磁场所掩盖，一些常规特征量，如磁异常总场强度（Total Magnetic Intensity，TMI）［6］、MGT分量、NSS、CT等，很难通过设置区域阈值等高线来圈定目标区域。当阈值设置过大对弱磁性和过埋深目标发现率低，而阈值设置过小会使得目标虚警率高，这为磁异常的自动检测带来困难。

倾斜角（θTilt）［19］是位场标势的总水平导数和垂向导数之比的反正切函数的角度值，表达式为：

其中：U为位场标势，∂U∂x、∂U∂y和∂U∂z分别表示U在三轴正交方向上的偏导数。在磁场中，U为磁标势，其偏导数即为磁场矢量。磁场倾斜角对源埋深和磁化强度不敏感，且为固定阈值的无量纲特征量，可同时描述区域内不同深度和强度的磁源，特别适合用于多目标数目和分布区域的自动检测。

倾斜角的缺点是容易受到磁化方向的影响，因此倾斜角能一定程度识别测区中埋深和磁矩不同的多个目标的数目和大致分布，但是分布区域受磁化方向影响将出现偏差。一般而言，MGT不变量中NSS和CT对磁化方向最不敏感［18］。在此基础上，本文将NSS、CT引入倾斜角公式，得到两个由不变量改进的倾斜角，表达式为：

其中：u为NSS，c为CT，其水平导数可利用空间域差分法估计得到；∂u∂z和∂c∂z不服从Laplace方程，故其不是NSS和CT的垂向导数［18］，但可利用频率域Fourier变换计算得到；θNSSTilt即为NSS不变量改进倾斜角，θCTTilt即为CT不变量改进倾斜角。θNSSTilt和θCTTilt无量纲，其取值范围为-π/4到π/4，致密源正上方处取值为π/4，且当观测点远离源上方区域后将迅速衰减至-π/4。θNSSTilt和θCTTilt对磁化方向、源距离均不敏感，当选定一固定等高线阈值时，该值能在区域MGTS检测极限范围内［5］快速圈定目标数目和分布位置，实现多目标预识别。

3 AFCM聚类算法用于多目标位置估计

航线测量中，目标单点定位算法可对孤立磁源进行连续追踪定位，连续计算的位置解将在目标空间区域附近形成稠密点云。将航线测量离散化，并将数条航线并列排开，即形成网格测量。相较于航线测量，网格测量观测点虽不连续，但在目标区域预识别基础上，仍能利用单点定位算法计算目标位置正上方识别区域内的节点的位置点云，这为多目标解集的AFCM聚类估计提供依据，如图1所示。

图1 识别区域内磁偶极子定位结果形成的稠密点云Fig.1 Dense point cloud formed by magnetic dipole positioning results in the recognition area

3.1 张量衍生不变关系的目标单点定位

尹等［20］提出了MGT几何不变量的概念，李青竹等［7］针对MGTS测量过程将这种量归纳为MGT衍生不变关系。若磁偶极子用六个未知量描述：位置矢量r=（x，y，z）T和磁矩矢量m=（mx，my，mz）T，其中r为偶极子到MGTS观测点间的相对坐标，则有如下关系。

3.1.1 磁矩夹角不变关系

其中，‖·‖是2范数，即m与r间夹角θ是不变的，且等于中间特征值λ3与归一化磁源强度u之比的反余弦函数。

3.1.2 特征向量不变关系

即中间特征值对应特征向量v3与m、r所在平面是垂直的，而最大、最小特征值对应特征向量v1、v2与m、r所在平面是共面的。

利用衍生不变关系可对区域内孤立磁偶极子进行定位，该方法对每个观测点提供了四种可能的坐标解，这些位置解关于观测点空间象限对称，且只有一个是真实解［7］：

3.2 磁源定位初始位置点云的AFCM聚类分析

当检测区域存在多个磁源时，衍生不变关系定位方法估计位置将近似聚集在磁源附近，图1中每个点云簇的质心可用于估计源的位置。聚类分析在解决许多问题中一直发挥着重要作用，而聚类数［16］是最重要的参数。计算的磁偶极子初始位置向量在目标真实位置附近聚集程度很高，而零散的虚假解和失败解可能极度偏离目标区域，这使得根据经验法则随机选择聚类质心时会出现虚假质心和位置偏移的可能。如果预定义的聚类数不正确，聚类结果可能会偏离给定数据集的真实结构。因此，实际应用中多目标自适应PR检测的关键即是利用算法自动找到的最优聚类数。

首先给出FCM聚类的计算过程。假定R=｛r1，r2，…，rn｝为总数为n的有限目标初始位置向量数据集，其中ri为第i个三维坐标向量；V=｛v1，v2，…，vc｝为总数为c的聚类簇质心，其中vi为第i个三维聚类中心；U=（uik）（n×c）为模糊划分矩阵，uik表示第i个样本属于第k簇的隶属（紧凑）程度。uik和vk的初始化及更新过程［17］遵循：

其中：ri-vk表示第i个向量到第k个簇质心之间的相对距离向量，m（m∈［1，∞））是模糊加权指数，控制成员的模糊性。m值越大由此产生的成员隶属便更加模糊，Pal等人［21］建议m应该在1.5和2.5之间，若无特殊要求通常m=2。FCM算法通过更新U和V，实现如式（9）所示目标函数的最小化迭代过程：

FCM聚类的具体步骤如下：

（1）初始化聚类的簇质心数c、模糊指数m、最大迭代次数itermax和阈值ξ；

（2）根据隶属度的约束随机初始化模糊分区U（0）；

（3）在第t步迭代，根据式（8）计算全部c簇质心V（t）；

（5）根据式（8）计算U（t+1），并返回步骤（3）；

（6）将当前对象按照最大隶属度原则排列进簇，并完成迭代。

FCM聚类需要预设聚类数c，AFCM算法可通过设置聚类数的搜索范围运行FCM生成不同聚类数的聚类结果，选择合适的聚类有效性函数评价聚类结果，最终确定FCM最优聚类数。1991年，Xie和Beni设计了一种基于模糊紧凑度和分 离度的CVF［22］。随后，Bensaid对其进行 改进，设计了对集群范围不敏感的CVF［23］，这种方法十分适用于磁源定位中初始位置点云集群分布不规律且集群容积大小不一的情形。然而，Bensaid方法在聚类数与真实对象数量接近时，聚类有效性函数值单调递减且接近于零，使簇质心数寻优过程会过早失去鲁棒性和决策功能。因此，我们在Bensaid方法基础上，通过引入惩罚函数消除簇数寻优过程的单调递减趋势，设计了一种新的基于模糊紧凑度和分离度的模糊聚类有效性函数实现AFCM聚类：

式（10）中的分子代表了第k簇的紧凑性，而分母代表了第k簇质心到其他簇质心的平均距离，用于分离性的度量。则此时分子分母之比代表了第k簇的聚类效果，L（U，V）即为全部聚类簇的聚类效果（比率）之和。显然L值随c的增长先变小后增大，而L取最小时对应的c值即为最优聚类数，这种决策模式能极大程度地排除磁偶极 子 定 位 的 虚 假 解 和 偏 差 解。Bensaid［23］和Sun［16］描述的迭代试错方法可以用于估计AFCM的最优聚类数量。

3.3 多目标磁矩矢量估计

AFCM聚类算法可以同步估计区域内多个磁偶极子的三维位置坐标。由于MGT分量随源的相对距离四次方衰减［18］，对于一切致密源上方的张量场受区域其他源和地磁场的影响可以忽略，即各磁偶极子正上方处测得MGT可认为仅由该源产生。设估计的源坐标rc=（xc，yc，zc），则源正上方处MGTS与源间相对矢量r0=（0，0，-zc），MGTS测得该点处张量矩阵为Gc，则根据磁偶极子正演方程和欧兰反褶积公式［24］可得

综上，基于AFCM聚类的多目标自适应PR检测方法技术路线如图2所示。

图2 磁偶极子多目标自适应PR检测技术路线Fig.2 Technical route of multi-target adaptive PR detection for magnetic dipoles

4 磁偶极子多目标的网格测量仿真

仿真环境为Matlab。在10 m×10 m的测区中，总共包含了20个磁偶极子，编号分别为1至20，其水平分布随机，埋深范围为-0.40 m至-0.95 m、磁矩强度为6 Am2至20 Am2，磁化倾角和偏角分别在［-90°，90°］和［-180°，180°］随机取值。图3显示了20个磁偶极子的分布情况，预设的物性参数列于表1中。

表1 预设和估计的20个磁偶极子目标的物性参数Tab.1 Preset and estimated physical parameters of the 20 magnetic dipoles

图3 多源磁偶极子目标示意图Fig.3 Schematic diagram of multi-target magnetic dipoles

其中r0=（x0，y0，z0）为位置坐标，M、Im和Dm分别为磁矩模、磁化倾角和偏角，磁矩m=（Mcos（Im）cos（Dm），Mcos（Im）sin（Dm），Msin（Im））。设置地磁场TMI为57 000 nT，设置测量平面z=0，坐标轴x正向向东，y正向向北，z正向向上，此时地磁倾角和偏角分别为-55°和-5°。利用平面十字形MGTS［7］在网格各节点处差分测量MGT分量，系统基线距离为0.3 m，网格间隔为0.1 m。该平 面MGT分量、MMA、NSS和CT如图4所示，其中（a）至（f）为测得MGT分量：（g）至（i）分别为计算的MMA、NSS、CT值。可见，由MMA、NSS和CT仅能发现部分强磁异常区域；磁梯度张量具有更高的分辨率，其中gzz分量识别效果最好但受磁化方向影响较大。此外，上述各参量无法找到有效等高线阈值判断磁性体是否存在。

图4 10 m×10 m测区平面内MGT、MMA、NSS和CT的识别结果Fig.4 Recognition results of MGT，MMA，NSS and CT in the 10 m×10 m survey area

提出的两个不变量改进倾斜角θNSSTilt和θCTTilt识别结果分别如图5（a）、（b）所示，图5（c）、（d）为设置等高线阈值为0时圈定的目标识别区域，图中用白圈标注了磁偶极子的真实位置。由图5可知，θNSSTilt和θCTTilt有效增强了测区平面内的较弱磁异常和深层磁异常目标。识别区域刚好处于磁源位置正上方处，且区域内节点数目充分，这为提高目标定位精度和实现位置点云AFCM聚类提供了可能。

图5 不变量改进倾斜角θNSSTilt和θCTTilt圈定的磁偶极子识别区域（等高线阈值0）Fig.5 Magnetic dipole identification area delineated by improved tilt angles θNSSTilt and θCTTilt（contour threshold 0）

选择θCTTilt识别结果作为目标区域，式（7）可计算出识别区域内各网格点对应的目标四组位置解。利用θCTTilt区域水平坐标和埋深坐标正负号筛选这些解，获得的初始位置解集的点云如图6所示，图6（a）为估计的磁性目标初始位置示意图，6（b）为其俯视图。可见，筛选后的解集二、三悬浮在源上方构成了大部分被遗漏虚假解，而筛选后的解集一、四更倾向于真实解，并位于在目标真实位置附近。但实际测量中目标真实位置是未知的，我们无法再继续筛选解集了，这导致了目标埋深的估计值可能总体性地偏小。

图6 识别区域内目标初始位置解集形成的稠密点云Fig.6 Dense point cloud formed by the solution set of the initial position of the target in the recognition area

提出的AFCM聚类算法对初始位置点云的聚类结果如图7所示，图中给出了估计的簇质心位置和目标真实位置。算法精确地估算了磁偶极子的数目，继而获得了各簇质心处估计的目标位置。表1中对比了估计的各目标三维坐标、磁矩强度和磁化方向，并统计了坐标值和磁矩强度的总体估计精度（Total Estimation Accuracy，PTEA）作为性能指标，表示参数真实值减去估计偏差量后与自身百分比的最小值：

其中：Pr（Pr≠0）为参数真实值；Pei为该参数第i（i=1，2，…，N）个估计值，N为目标数量。为更真实模拟实际测量场景，在相同条件下给每个MGT分量加入了相互独立的均值为0 nT/m、方差分别为2 nT/m和5 nT/m的高斯白噪声。两次测试的TEA分别记为TEAd1、TEAd2，见表1，同时给出了计算的目标磁矩和坐标估计偏差。结果显示在题设噪声环境下的目标数目识别仍然准确，该聚类结果仍具有较好的稳定性和鲁棒性。

5 多目标小型磁铁的网格测量实验

5.1 2.1 m×2.1 m测区五磁体定位

在中国石家庄某地平坦的地面进行实验，测量区域为2.1 m×2.1 m，x轴正向向北，y轴正向向东，z轴正向向上。实验时间为夏季某日18：00，环境温度28℃，系统为基线距离0.5 m的平面十字形MGTS［5-7］，传感器型号Mag-03（英国Bartington），采集卡NI-6281，采样率500 Hz。实验如图8所示，在测区内共放置序号为1～5的五个小型磁源，其磁矩强度和方向未知，采样间隔为0.1 m，各网格点采样时间不少于2 s。由于设备电磁噪声频段主要集中在50 Hz、200 Hz等频段（国内交流电场频率），为尽量避免电磁场影响和保留真实磁异常信号，利用零相位低通滤波器进行降噪，并将截止频率保守地设置为5 Hz［5］。预先测得该区域内地磁场TMI均值为53 394.36 nT，地磁倾角和偏角分别为56.33°（向下）和-20.17°（自北偏西）。

图8 实验一：五个磁铁的网格定位实验Fig.8 Experiment 1：grid positioning of five magnets

测得该平面处的MGT分量和常规特征量如图9所示。图9（a）至（f）为测得MGT分量；图9（g）至（i）分别为计算的MMA、NSS、CT值。图10为θNSSTilt和θCTTilt的区域识别效果，白圈标注了磁体的真实位置。可见，五号磁铁相对于其他目标有更大的磁矩，导致1～4号的异常场几乎被悉数淹没。MGT和常规特征量已无法对全部磁铁进行有效识别，而设定阈值为0的两个改进倾斜角能准确圈定出五个磁体的分布位置，其中θCTTilt识别效果更佳。

图9 五个小型磁铁的MGT、MMA、NSS和CT识别结果Fig.9 MGT，MMA，NSS and CT recognition results of five small magnets

图10 θNSSTilt和θCTTilt圈定的5个磁铁目标的识别区域（等高线阈值0）Fig.10 Recognition area of 5 magnet objects delineated by θNSSTilt and θCTTilt（contour threshold 0）

而后，用张量衍生不变关系方法估计识别区域内各网格节点处的目标单点定位结果，对得到的四组解集进行区域和象限筛选。为了获得更高的聚类精度，在利用θCTTilt区域筛选定位解集时，将θCTTilt阈值上调至0.5，此时各区域内单点定位解集形成的稠密点云如图11（a）和（b）所示，AFCM聚类的结果如图11（c）和（d）。由结果可知，经区域识别、定位筛选的目标解集聚集在真实位置附近，AFCM准确估计了磁体数目为5，簇质心的位置估计列于表2中，并给出了计算的目标坐标估计偏差。该实验对区域内的五个小型磁体的位置坐标估计偏差小于0.061 m，获得了令人满意的结果。

表2 定位实验中对九个磁性体的真实位置和估计的物性参数Tab.2 Real positions and estimated physical parameters of the nine magnets in the positioning experiment

图11 识别区域内的单点定位密集点云与AFCM聚类结果Fig.11 Single-point positioning dense point cloud and AFCM clustering results in the recognition area

5.2 1.2 m×1.2 m测区四磁体定位

在同地域进行第二次实验，测量区域大小为1.2 m×1.2 m，时间为春季某日16：00，环境温度17℃，测量仪器更换为五个HSF-112型（中国Huashun）磁通门传感器搭建的平面十字形MGTS，采集卡为NI-9239，采样率为1 kHz。实验如图12所示，测器内共放置序号6～9的四个小型磁源（6～9号磁铁磁矩强度、方向和形状均不相同［5］），磁矩强度和方向未知，预设坐标列于表3中。采样间隔为0.05 m，各节点采样时间和降噪方法不变。测得该区域内地磁总场强度 均 值53 407.02 nT，磁 倾 角-56.17°，磁 偏角-19.98°。

图12 实验二：四个磁体的网格定位实验Fig.12 Experiment 2：grid positioning of four magnets

为了获得更加清晰分明的区域识别结果，对张量数据网格进行0.01 m间隔的三次样条插值处理，因磁源埋深增加且各目标间磁矩有明显差异，将倾斜角阈值设定为0.3，MGT分量及其特征量如图13所示，θNSSTilt和θCTTilt的区域识别如图14所示。

图13表明6和9号磁铁的磁矩强度明显弱于7和8号，且6和9号的磁张量场已完全被强磁场淹没；7和8号的磁矩方向具有明显差异。此时传统方法已经失去识别目标边界的能力，而图14中的θNSSTilt和θCTTilt仍能清晰地识别出四个磁体的大致分布区域。根据识别结果，将等高线上调至0.4并选择θNSSTilt识别区域筛选定位解集。各区域内单点定位解集的AFCM聚类结果如图15所示。

图15 识别区域内4个磁铁的AFCM聚类结果Fig.15 AFCM clustering results of four magents in the recognition area

结果显示在较弱磁性体磁场被淹没情况下，两种不变量改进倾斜角仍能对磁体进行区域预识别；尽管4个磁体磁矩、形状差距较为明显，经区域识别、定位筛选的目标解集仍聚集在磁体真实位置附近，AFCM准确估计了磁体数目为4，估计的簇质心位置和磁矩列于表2中；磁体8的估计磁矩强度数倍强于磁体6和9，且磁体间磁化方向差距较大，此时实验对区域内的四个小型磁体的位置坐标估计偏差仍小于0.091 m，同样令人满意。

6 结论

本文提出了一种利用AFCM聚类算法和张量不变量改进倾斜角定位多目标磁偶极子的方法。通过引入惩罚函数消除簇数寻优过程的单调递减趋势，设计了一种新的基于模糊紧凑度和分离度的模糊聚类有效性函数实现AFCM聚类。由归一化磁源强度和张量缩并两个张量不变量改进的倾斜角可用于目标分布区域识别和初始位置坐标筛选，而张量衍生不变关系定位方法用于计算识别区域内的各网格节点处对应的四组位置解集，这些解集将在目标真实位置附近形成稠密点云。设计的AFCM算法对这些点云进行三维聚类并自动估计簇质心数目和质心位置。

仿真结果证明了该方法的有效性和鲁棒性；在2.1 m×2.1 m和1.2 m×1.2 m的网格测量实验中，该方法对埋深不超过1 m、磁矩强度和方向差异较大的九个小型磁铁的坐标估计偏差小于0.091 m。然而，由于衍生不变关系单点定位算法的虚假解不能完全筛除，目标估计精度受到了一定的限制，发展更具优势和更高精度的磁偶极子单点定位方法将进一步提高本文方法的定位性能和适用性。