张坚，林春生，周建军

(1.海军青岛雷达声纳修理厂，山东 青岛 266100；2.海军工程大学兵器工程系，湖北 武汉 430033)

航空磁性探测广泛应用于考古、矿物勘探和水中磁性目标探测等领域[1].但是在飞机上对地磁异常场进行测量时，飞机干扰磁场叠加在地磁异常场上，掩盖了地磁异常场原有的特性，所以必须对飞机干扰磁场进行补偿，利用测量数据求解飞机干扰磁场模型系数，使用所求系数计算干扰磁场，最后从测量的总磁场中减去干扰磁场，完成飞机磁补偿.因此，飞机干扰磁场模型的求解是飞机磁补偿中的关键技术.

美国学者Tolles和Lawson在上世纪40年代最早对飞机干扰磁场的数学模型进行了研究，建立了经典的Tolles和Lawson模型，但没有给出模型的求解方法[2].其后，Leliak和Bickel均对该模型进行了改进[3-4]，但其求解精度仍然不高.飞机干扰磁场模型属于线性模型，并且存在较强的复共线性.线性模型求解的常用方法有最小二乘估计LSE和岭估计RE等[5-6].最小二乘估计运算量小，但在模型存在复共线性时均方误差很大；岭估计可以在一定程度上克服模型复共线性，但在实际探测中，飞机上电气设备(如继电器、电机等)产生较强的与飞机姿态变化无关的干扰磁场，这种测量噪声叠加在观测数据上无法进行补偿，并会对岭估计的求解效果产生较大影响.针对飞机干扰磁场模型求解时均方误差较大的问题，本文中提出小波阈值去噪结合岭估计(WTRE)的模型参数求解方法，并通过参数估计和飞机磁补偿的仿真实验对方法进行了验证.

1 飞机干扰磁场建模

图1 飞机干扰磁场模型坐标系

飞机干扰磁场包括剩余磁场、感应磁场和涡流磁场，本文中采用经典的Tolles和Lawson模型进行飞机干扰磁场的建模，并对模型进行简化，模型坐标系如图1所示.图中，机载磁传感器为坐标原点，地磁场矢量He与X，Y，Z轴的夹角为(X0，Y0，Z0)，其方向余弦定义为(cosX0，cosY0，cosZ0).

1)剩余磁场Hp是由飞机上铁磁性物质的剩余磁化产生的，在短时间内可以视为稳定不变，建模时将其分解为与X，Y，Z轴平行的三分量Hpx，Hpy，Hpz，即Hp=(HpxHpyHpz)T,

则剩余磁场系数有3项.

其中，K为感应磁场系数矩阵，有9项系数.

其中，L为涡流磁场系数矩阵，有9项系数.

上述模型是传统的飞机干扰磁场模型，共有21项模型系数.利用感应磁场系数矩阵K的对称性[3]，以及方向余弦cos2X0+cos2Y0+cos2Z0=1的性质，可将3种磁场相加后再进行合并与化简，最终得到总的飞机干扰磁场ΔH为ΔH=ABTHe+kzzHe

(1)

A=(Hpx，Hpy，Hpz，kxx-kzz，kxy+kyx，kzx+kxz，kyz+kzy，kyy-kzz，lxx-lzz，lyx，lzx，lxy，lyy-lzz，lzy，lxz，lyz)

B= (cosX0/He，cosY0/He，cosZ0/He，cos2X0，cosX0cosY0，cosX0cosZ0，cosY0cosZ0，cos2Y0，

(1)式是简化后的飞机干扰磁场模型，共有16项系数.飞机干扰磁场的补偿过程是，利用测量数据估计线性方程组(1)中的16项未知系数，再根据估计出的系数计算飞机干扰磁场，从磁场测量值中减去飞机干扰磁场，得到补偿后的磁场值.

2 小波阈值去噪结合岭估计的模型求解方法

(2)

式中，Y为n×1观测向量；X为n×m列满秩设计矩阵，为已知量；β为m×1待求的模型参数；e为n×1随机误差向量；In为n阶单位阵；σ2为方差参数.

对于飞机干扰磁场模型，自变量X由方向余弦、方向余弦的导数及其乘积项构成，Y由磁场总强度的测量值构成，e为测量误差向量.因为在一定范围内地磁偏角和倾角可看做不变，显然Y是飞机姿态角变化的函数.用LS估计求解模型时均方误差很大，原因在于设计矩阵XTX的奇异性.该矩阵产生奇异性的原因是，在保证飞机稳定飞行的情况下，飞机的机动角度不可能很大，以致XTX每一行的对应元素值非常接近，每行之间有比较强的相关性，导致飞机干扰磁场模型存在严重的复共线性.

2.3小波阈值去噪与岭估计的结合小波阈值去噪是Donoho提出的一种简洁有效的去噪方法，通过对小波系数进行非线性阈值处理，恢复噪声中的信号.算法的基本思想如下：属于Besov空间的信号在小波域的能量主要集中在有限的几个小波系数中，而白噪声的能量却分布在整个小波域内.因此对wj，k，可以找到一个合适的阈值λ，当wj，k<λ时，认为此时的wj，k主要由噪声引起，将其置为零；当wj，k>λ时，认为此时的wj，k主要由信号引起，将这部分wj，k直接保留下来(硬阈值法)或者按某一固定量向零收缩(软阈值法).最后将处理过的小波系数进行小波重构而得到去噪后的信号.

在航空磁性探测过程中，飞机上电气设备(如继电器、电机等)会产生较强的干扰磁场，这种测量噪声与飞行姿态无关，叠加在观测数据上将对模型求解产生较大的影响.由于测量噪声在有限频带内可以近似看作白噪声，因此本文中采用小波阈值去噪作为岭估计前的预处理手段，对测量噪声进行滤除.小波阈值去噪结合岭估计(WTRE)的实施过程如下：求解模型参数之前，采用小波阈值去噪对观测数据进行预处理，降低电气磁场等测量噪声对模型求解的影响，阈值函数选取软阈值；然后利用岭估计求解干扰磁场模型系数，岭参数的选取采用L曲线法.

3 算法验证

3.1参数估计仿真实验采用Matlab7.0软件仿真验证WTRE在飞机干扰磁场模型求解中的有效性.设地磁场大小为5×104nT，地磁场方向余弦可以用地磁场方位角θ和俯仰角φ的函数来表示，且假设dφ/dt=1，dθ/dt=cos(t/2)，φ∈[0，2π]，θ∈[0，π].观测数据采样点数为400，即N=400.设定模型系数的真值如下：剩余磁场系数[10，35，90]，5个感应磁场系数均为0.5，8个涡流磁场系数均为0.2.利用设定的系数真值，仿真产生信噪比分别为SNR1=21.47 dB和SNR2=9.13 dB时的飞机干扰磁场信号，如图2所示.

图2 不同信噪比下飞机干扰磁场信号仿真

衡量矩阵病态程度的方法之一是求它的条件数.计算得到法矩阵XTX的条件数为1.35×1015，呈严重病态性.分别采用LS估计(LSE)、岭估计(RE)和WTRE 3种方法进行求解，计算结果如表1所示，

表1 LS估计、岭估计和WTRE的计算结果

岭估计采用L曲线法确定k值.采用均方误差(MSE)作为求解效果的评价标准：

由表1的计算结果可以得出以下结论：

1)LS估计的均方误差远大于岭估计，原因是法矩阵XTX的条件数很大，模型复共线性严重.

2)利用岭估计进行飞机干扰磁场模型参数的求解可以较好地克服复共线性的影响，且基于L曲线法的岭参数选取方法适用性好、稳定性高，能够满足航空磁性探测的实际需求，估计精度较高.

3)在不同信噪比下，WTRE的均方误差均小于RE，特别是在低信噪比条件下，WTRE对RE的改进更加明显，这是由于WTRE降低了测量噪声对模型求解的影响.

图3 磁偶极子目标信号

式中，μ0=4π×10-7H/m为真空磁导率，M为磁性目标的磁矩矢量，r为目标到磁传感器的距离.通过设定一组参数，可以由上述磁偶极子模型仿真生成目标磁场信号，如图3所示.目标信号与飞机干扰磁场叠加后的混合信号如图4(a)所示.由图4(a)可以看出，飞机干扰磁场的变化幅度超过10 nT，目标信号完全被干扰磁场淹没，无法辨识.分别采用LS估计、岭估计和WTRE所求系数进行磁补偿后的信号如图4(b)～(d)所示.可以看出，LS估计补偿后的信噪比最低，WTRE补偿后的信噪比最高.3种方法中只有WTRE的补偿结果可以准确检测到目标信号，其他2种方法都在不同程度上受到噪声的干扰.这进一步验证了WTRE在飞机磁补偿中的有效性.

(a)目标与飞机干扰磁场的混合信号

(b)LS估计的补偿结果

(c)岭估计的补偿结果

(d)WTRE的补偿结果

4 结论

针对飞机干扰磁场模型求解时均方误差较大的问题，本文中提出小波阈值去噪结合岭估计(WTRE)的模型参数求解方法.首先采用小波阈值去噪对观测数据进行预处理，降低测量噪声对模型求解的影响；然后利用岭估计求解干扰磁场模型系数，岭参数的选取采用L曲线法.仿真结果表明，WTRE的均方误差明显小于LS估计和岭估计，特别是在低信噪比情况下对岭估计的改进尤为明显，采用WTRE所求系数补偿后的信噪比最高，在飞机磁补偿中与其它方法相比具有明显的优越性.

[1] 熊盛青. 我国航空重磁勘探技术现状与发展趋势[J].地球物理学进展， 2009，24(1):113-117.

[2] Tolles W E， Lawson J D. Magnetic compensation of MAD equipped aircraft[R].New York:Airborne Instruments Lab Inc，1950.

[3] Leach B W. Aeromagnetic compensation as a linear regression problem[C]//Information linkage between applied mathematics and industry II. New York:Academic Press Inc， 1980：367-371．

[4] Bickel S H. Small signal compensation of magnetic fields resulting from aircraft maneuvers[J].IEEE Trans on AES，1979， 15(4):515-525.

[2] 陈希孺， 王松桂. 近代回归分析-原理方法及应用[M].合肥:安徽教育出版社， 1987:294-320.

[6] 刘莉. NQD样本下部分线性模型中估计的强相合性[J].湖北大学学报：自然科学版， 2004， 26(4):290-293.

[7] Hansen P C. Analysis of discrete Ill-posed problems by means of the L-curve[J].SIAM Review，1992，34(4):561-580.

[8] 王振杰， 欧吉坤. 用L-曲线法确定岭估计中的岭参数[J].武汉大学学报：信息科学版， 2004， 29(3):235-238.

[9] 林春生， 龚沈光. 舰船物理场[M].北京:兵器工业出版社， 2007.