总强度磁异常各阶垂向导数换算新方法

翟国君，卞光浪，2，黄谟涛

1.海军海洋测绘研究所，天津300061；2.大连舰艇学院海测工程系，辽宁大连116018

总强度磁异常各阶垂向导数换算新方法

翟国君1，卞光浪1，2，黄谟涛1

1.海军海洋测绘研究所，天津300061；2.大连舰艇学院海测工程系，辽宁大连116018

总强度磁异常垂向导数换算在磁性目标解释推断过程中具有重要意义。分析频率域总强度磁异常各阶垂向导数转换因子的滤波特性，指出在求解总强度磁异常垂向导数时，常规傅里叶变换法会明显放大观测数据中高频噪声的影响，甚至可能淹没掉真实信息。给出一个重要定理：总强度磁异常沿垂直方向的积分及各阶偏导数均为准调和函数，在此基础上，提出基于拉普拉斯方程的逐步求解法，联合空间域和频率域换算总强度磁异常沿垂直方向的各阶导数。同时，采用球体磁场模型验证所提方法的有效性，并推导地磁场方向和磁化强度方向不一致时，球体总强度磁异常沿垂直方向的一阶和二阶导数表达式。结论表明：利用样条法计算的垂向导数结果精度明显优于常规傅里叶变换法计算结果，且具有较强的抗噪能力。

总强度磁异常；垂向导数；空间域和频率域；拉普拉斯方程；球体磁场模型

1 引 言

磁测方法早期主要应用于地球物理和磁力勘探领域。近年来，随着人类开发与利用海洋资源进程的加快，该方法以成本低、快速高效等突出优点，在水下小尺度磁性目标精密探测方面备受关注，特别是沉船、铁锚、海底电缆、未爆炸物以及水下潜航器的定位与识别方面更具明显优势［1－2］。

磁异常导数被广泛应用于磁性目标磁异常解释，是压制区域场、圈定局部场、分离叠加异常以及异常反演中常用的方法［3］。例如，利用解析信号法和欧拉反褶积法确定磁性目标平面位置，或者采用泰勒级数法实现磁场空间稳健延拓，都需要以高精度磁异常水平和垂直方向的一阶、二阶甚至高阶导数作为基础数据。现阶段，获取磁场垂向导数主要有两种方式，其一是将磁力传感器设计为磁阵列模式［4］，利用微分原理直接计算磁场的垂向导数，该模式对载体运行姿态要求较高，海洋上受风、流、浪等因素影响，实现起来较为困难；其二是将采集的剖面或平面磁场信息进行数据转换处理，这种转换既可在空间域内进行，也可在频率域内进行。在空间域换算方面，文献［5］尝试将样条函数插值法引入位场垂向导数解算中，但计算精度有限。目前普遍采用的磁异常垂向导数换算方法是傅里叶变换法［6］（简记为“常规FFT法”），该方法理论较为完善，但其频率响应函数相当于高频放大器，对观测数据中的噪声具有显著的放大作用，当观测数据含有一定误差时，计算精度较低，且存在较强的边界效应。为避免复杂的复数运算，文献［7］提出采用离散余弦变换（DCT）计算磁异常导数，该方法去除了原信号的相关性，在二度磁性目标磁异常垂向导数换算中应用效果较好，关于该方法能否适用于三度磁性目标磁异常垂向导数换算还有待于进一步研究。

针对总强度磁异常沿垂直方向导数的精确求解问题，笔者给出了关于总强度磁异常的一个重要定理，并综合利用空间域和频率域运算，提出一种求解总强度磁异常各阶垂向导数的新方法。文中详细研究了该方法的基本原理和技术措施，通过仿真球体磁场模型检验了方法的有效性，并将其计算结果与常规FFT法进行了分析比较。

2 常规FFT法

由傅里叶变换，三度磁性目标总强度磁异常Bm的一阶垂向导数Bmz与其频谱关系［8］为

根据位场理论，总强度磁异常近似满足拉普拉斯方程，利用外部Dirichlet问题条件，其空间延拓公式可表示为

式中，FBm（u，v，0）为z＝0平面上总强度磁异常Bm（x，y，0）的频谱。

将式（2）对z方向求偏导数，得

由于同一函数有相同的傅里叶表达式，对比式（1）和式（3）两端，得

当z＝0时，式（4）即为平面上总强度磁异常频谱向其一阶垂向导数频谱转换的关系式，同理，还可推导出总强度磁异常沿垂向的第k阶导数频谱为这里称为总强度磁异常与其第k阶垂向导数间的频率域转换因子（也称频率响应函数）。计算时，首先将平面上总强度磁异常格网数据利用快速傅里叶变换转换为频谱FBm（u，v，0），再乘以垂向导数频率域转换因子获取第k阶垂向导数的频谱，最后通过反傅里叶变换即可得到总强度磁异常各阶垂向导数。

3 基于拉普拉斯方程的逐步求解法

为减弱高频噪声对总强度磁异常垂向导数计算的影响，这里提出基于拉普拉斯方程的逐步求解法。在求解过程中，将运用总强度磁异常的某些特殊性质，限于篇幅，略去推导过程，给出如下定理。

定理：在无源空间内，总强度磁异常沿垂直方向的积分及任意阶偏导数均为准调和函数。

结合上述定理，总强度磁异常各阶垂向导数的求解步骤如下：

（1）对z＝0平面上总强度磁异常数据进行傅里叶变换，设其在频率域频谱为FBm（u，v，0），由式（4），FBm（u，v，0）乘以频率域转换因子即为总强度磁异常沿垂直方向积分的频谱，再利用反傅里叶变换得到积分值J（x，y，0）；

（2）计算J（x，y，0）沿x方向和y方向的二阶偏导数Jxx（x，y，0）与Jyy（x，y，0），由定理知，总强度磁异常沿垂直方向的积分是准调和函数，满足拉普拉斯方程，则总强度磁异常沿垂直方向的一阶偏导数为

（3）计算Bm在水平方向的两个二阶偏导数与而总强度磁异常自身也为准调和函数，同样利用拉普拉斯方程，得总强度磁异常沿垂直方向的二阶偏导数为

（4）由于总强度磁异常沿垂直方向的任意阶偏导数都是准调和函数，得第k（k＝3，4，…，n）阶垂直导数为

式中，Bmkz为总强度磁异常沿垂直方向的第k阶偏导数和分别是Bm（k－2）z沿x方向和y方向的二阶偏导数。

4 调和函数二阶水平方向导数计算方法

目前，计算调和函数沿水平方向的二阶偏导数通常在频率域内进行转换（计算步骤与常规FFT法计算垂向导数类似），但该方法存在的不足是会显著放大观测数据中噪声的影响，且具有较强的边界效应。在数值计算中，三点二阶中心差分法和三次样条曲线函数法经常被用于计算函数在一维方向的二阶导数，本文根据平面磁场数据的二维空间分布特点，将两者引入并应用到调和函数二阶水平方向偏导数计算中。下面以总强度磁异常为例，给出三点二阶中心差分法和双三次样条曲线函数法计算其沿水平方向二阶偏导数的具体方法原理。

4.1 三点二阶中心差分法

假设x1，x2，…，xM为格网化数据沿x方向的等间距节点，且x1＜x2＜…＜xM，xi－xi－1＝Δx。则在（xi，y）节点处总强度磁异常沿x方向的二阶偏导数可表示为

类似地，设y1，y2，…，yN为格网化数据沿y方向的等间距节点，且y1＜y2＜…＜yN，yiyi－1＝Δy。则在（x，yi）节点处总强度磁异常沿y方向的二阶偏导数可表示为

值得注意的是，利用式（8）和式（9）求解总强度磁异常沿水平方向二阶偏导数时，边界处的二阶导数需要附加条件，本文仿真计算时，考虑到仅对比梯度变化较大区域的计算精度，所以未计算边界部分的二阶偏导数。

4.2 双三次样条曲线函数法

双三次样条曲线函数是对x方向和y方向分别采用三次样条曲线函数进行拟合［9］。以x方向测线为例，设有观测数据 ｛x1，Bm（x1）｝，｛x2，Bm（x2）｝，…，｛xM，Bm（xM）｝，且x1＜x2＜…＜xM，整个区间［x1，xM］被节点分成M－1个小区间［x1，x2］、［x2，x3］、…、［xM－1，xM］，现要求构造函数S（x）满足以下三个条件：

（1）插值条件

（2）在每个小区间［xi，xi＋1］（i＝1，2，…，M－1）上，S（x）是x的三次多项式，即

（3）在整个区间［x1，xM］上，S（x）具有连续的二阶导数，在区间的每个节点xi（i＝1，2，…，M）上，函数值、一阶导数和二阶导数连续，即

式中，i＝2，3，…，M－1。

利用双三次样条曲线函数计算水平方向二阶偏导数时，需要给定两个边界约束条件，本文选择自然三次样条插值函数，即Sxx（x1）＝Sxx（xM）＝0，结合方程组（12）得到的Sxx（xi）即为xi（i＝2，3，…，M－1）节点处总强度磁异常沿x方向的二阶偏导数。按照相同的方法，可进一步得到各测点处总强度磁异常沿y方向的二阶偏导数。

5 仿真试验分析

为检验提出的总强度磁异常各阶垂向导数换算方法的有效性，仿真三个磁化球体磁异常数据进行分析。

仿真测区面积为295m×295m，测线和测点间距均为5m。地磁场倾角和偏角分别为45°与5°，为检验方法的普遍适用性，这里设球体的磁化强度方向和地磁场方向不一致，球体空间位置及磁性参数分别如表1所示。

表1 仿真球体模型相关参数Tab.1 Parameters of synthetic sphere model

为定量分析总强度磁异常的垂向一阶导数、二阶导数计算结果精度，在文献［8］给出的球体总强度磁异常表达式基础上，笔者进一步推导磁化强度方向和地磁场方向不一致时，磁化球体总强度磁异常垂向一阶导数Bmz和二阶导数Bmzz解析表达式，分别如式（13）与式（14）所示

式中，I0和A0分别为地磁场倾角和偏角；Bxz、Byz、Bzz和Bxzz、Byzz、Bzzz分别是Bv（v＝x，y，z）沿垂直方向的一阶偏导数与二阶偏导数，其中，

根据球体总强度磁异常以及式（13）和式（14），仿真得到z＝0m平面上总强度磁异常及其垂向一阶和二阶偏导数等值线如图1所示。

图1 总强度磁异常及其垂向一阶和二阶偏导数理论等值线图Fig.1 Contour map of theoretical TMA as well as its 1st and 2ed vertical derivatives

图1（a）中总强度磁异常最大值为48.96nT，最小值为－20.69nT，平均值为1.05nT；图1（b）中总强度磁异常垂向一阶导数最大值为5.30nT／m，最小值为－2.18nT／m，平均值为0.001 8nT／m；图1（c）中总强度磁异常垂向二阶导数最大值为0.70nT／m2，最小值为－0.30nT／m2，平均值为－0.000 1nT／m2。

本文在前面给出了两种垂向导数计算方法，分别是常规FFT法和基于拉普拉斯方程的逐步求解法。在利用后者求解垂向导数时，提供了两种求解水平方向二阶偏导数的算法，简称为“差分法”和“样条法”。为充分验证提出方法的适用性，这里分无噪声和有噪声两种情况分别进行分析。

5.1 无噪声情况

图2和图3是根据图1（a）中的总强度磁异常数据，利用常规FFT法、差分法和样条法换算得到总强度磁异常沿垂直方向的一阶偏导数和二阶偏导数等值线图。

图2 三种方法计算得到的总强度磁异常沿垂直方向的一阶偏导数等值线图（等值线间隔0.2nT／m）Fig.2 Contour map of calculated the 1st vertical derivative of TMA with three methods（contour interval is 0.2nT／m）

图3 三种方法计算得到的总强度磁异常沿垂直方向的二阶偏导数等值线图（等值线间隔0.05nT／m2）Fig.3 Contour map of calculated the 2ed vertical derivative of TMA with three methods（contour interval is 0.05nT／m2）

从图2和图3可以看出，无噪声情况下，三种换算方法得到的垂向一阶偏导数和二阶偏导数等值线图与相应的理论等值线图实现了很好的吻合，但常规FFT法随着导数阶数的增加，边界效应越来越明显。为定量计算总强度磁异常各阶垂向导数换算精度，并排除边界效应对统计结果的影响，选取上述各图中虚线框D域内梯度变化较大的格网点数据进行分析，即D＝｛（x，y）｜25m≤x≤270m，25m≤y≤270m｝，采用的精度评价指标公式为式（15）。

式中，＾Bmkz（i，j）和Bmkz（i，j）（k为正整数）分别表示格网点上总强度磁异常第k阶垂向导数计算值与理论值；Nx和Ny分别为D域内x方向和y方向采样点数；ρk为第k阶垂向导数计算结果中误差。

三种方法换算结果统计情况如表2所示。

表2 三种方法换算的垂向一阶和二阶导数结果统计Tab.2 Statistics of the 1st and 2ed vertical derivative results with three methods

从表2可以看到，三种垂向导数换算方法计算结果均具有较高的精度，其中，样条法和常规FFT法计算结果更加接近于理论值。由于式（15）给出的精度评价指标属于绝对误差，其大小与各阶垂向导数的振幅紧密相关，垂向导数的振幅越大，绝对误差越大；反之，振幅越小，绝对误差则越小。因此，该指标在一定程度上无法全面反映换算结果的优劣，故进一步采用中误差与理论值振幅之比作为评价指标

以样条法为例，总强度磁异常垂向一阶导数和二阶导数理论值振幅分别为7.48nT／m和1nT／m2，计算得到的垂向一阶、二阶导数中误差和振幅之比分别为0.11%与0.17%，验证了该方法具有很高的垂向导数换算精度，且二阶垂向导数换算结果的相对精度要稍低于一阶垂向导数换算结果。因此，当实测资料中无噪声或噪声微小时，既可采用常规FFT法，也可采用样条法计算总强度磁异常的垂向一阶导数和二阶导数。

考虑到利用泰勒级数法进行总强度磁异常空间延拓时，往往还需要获取总强度磁异常的高阶垂向导数，本文进一步对比分析了常规FFT法和样条法的垂向三阶、四阶和五阶导数换算结果。结果表明：常规FFT法计算结果存在较强的边界效应，且随着垂向导数阶数的增加，其边界效应更加明显，到第五阶垂向导数时，误差已非常大；而样条法虽然也存在一定的边界效应，但其换算效果要明显优于常规FFT法结果。因此，实践中建议采用样条法计算总强度磁异常高阶垂向导数。

5.2 含有噪声情况

为更接近实际情况，进一步检验提出方法的稳健性，在z＝0m平面上总强度磁异常理论值中加入2%高斯白噪声（振幅为1.39nT），仿真得到含噪声的总强度磁异常等值线如图4所示。

图4 含有高斯噪声的总强度磁异常等值线间隔（等值线间隔2nT）Fig.4 Contour map of TMA corrupted with Gauss noise（contour interval is 2nT）

图5和图6是用含有高斯噪声的总强度磁异常数据换算得到其沿垂直方向的一阶导数和二阶导数等值线图。

如图5和图6所示，当总强度磁异常数据中加入2%高斯白噪声后，垂向一阶导数和二阶导数计算结果精度要明显低于无噪声情况下的结果，表明噪声对垂向导数换算结果存在较大影响。在三种方法中，样条法计算结果精度最高，根据D域内格网点数据，利用式（15）和式（16）计算得到的ρ1＝0.07nT／m，η1＝0.93%，ρ2＝0.05nT／m2，η2＝5%；差分法计算结果精度次之，对应ρ1＝0.13nT／m，η1＝1.7%，ρ2＝0.11nT／m2，η2＝11%；常规FFT法计算结果精度最低，其中ρ1＝0.21nT／m，η1＝2.8%，ρ2＝0.13nT／m2，η2＝13%。从相对精度指标变化情况来看，垂向二阶导数较一阶导数受噪声的影响更加严重。

图5 三种方法计算得到的总强度磁异常沿垂直方向的一阶偏导数等值线图（等值线间隔0.2nT／m）Fig.5 Contour map of calculated the 1st vertical derivative of TMA with three methods（contour interval is 0.2nT／m）

这里列出零噪声、1%、2%、3%、5%和10%等不同噪声水平情况下，三种垂向导数方法计算得到的垂向一阶、二阶导数结果精度。

图6 三种方法计算得到的总强度磁异常沿垂直方向的二阶偏导数等值线图（等值线间隔0.05nT／m2）Fig.6 Contour map of calculated the 2ed vertical derivative of TMA with three methods（contour interval is 0.05nT／m2）

从表3可以看出，当总强度磁异常观测数据中含有噪声时，同一噪声水平情况下，样条法计算结果精度最高，差分法次之，常规FFT法精度最低；而对于同一种方法，垂向二阶导数计算结果相对精度要低于一阶导数计算结果相对精度。同时，随着噪声水平的逐渐增大，各阶垂向导数换算结果精度也不断降低。

表3 不同噪声水平下三种方法换算得到的垂向导数中误差Tab.3 Mean square error of three methods with different noise level

6 结 论

（1）分析了常规FFT法在换算总强度磁异常各阶垂向导数过程中存在的不足；给出了一个重要定理：在无源空间内，总强度磁异常沿垂直方向的积分和各阶偏导数同样为准调和函数，在此基础上，提出了一种新的总强度磁异常各阶垂向导数换算方法。

（2）围绕调和函数的二阶水平方向偏导数求解问题，给出了三点二阶中心差分法和双三次样条曲线函数法。通过仿真试验对三点二阶中心差分法和双三次样条曲线函数法应用效果进行了比较分析，结果表明，利用样条法换算得到的各阶垂向导数精度较高，且具有较强的抗噪声能力。

（3）当总强度磁异常观测数据不含噪声时，常规FFT法、差分法和样条法换算的总强度磁异常垂向一阶导数与二阶导数结果精度大体相当，但随着其垂向导数阶数的增加，常规FFT法边界效应现象较其他两种方法严重得多。当总强度磁异常观测数据含有噪声时，随着噪声水平的不断提高，垂向导数换算结果精度逐渐降低。但总体而言，样条法垂向导数换算结果精度要明显优于常规FFT法换算结果。

［1］ GAO Ping，LESLIE M C.A Two－dimensional Generalized Likelihood Ratio Test for Land Mine and Small Unexploded Ordnance Detection［J］.Signal Processing，2000，80（8）：1669－1686.

［2］ YU Bo，ZHAI Guojun，LIU Yanchun，et al.The Downward Continuation Method of Aeromagnetic Data to the Sea Level［J］.Acta Geodaetica et Cartographica Sinica，2009，38（3）：202－209.（于波，翟国君，刘雁春，等.利用航磁数据向下延拓得到海面磁场的方法［J］.测绘学报，2009，38（3）：202－209.）

［3］ WEN B D，SHU K H，YI C Y.A Derivative－based Interpretation Approach to Estimating Source Parameters of Simple 2DMagnetic Sources from Euler Deconvolution，the Analytic－signal Method and Analytical Expressions of the Anomalies［J］.Geophysical Prospecting，2007，55（2）：255－264.

［4］ SALEM A，HAMADA T，ASAHINA J K.Detection of Unexploded Ordnance（UXO）Using Marine Magnetic Gradiometer Data［J］.Exploration Geophysics，2005，58（1）：97－103.

［5］ WANG Bingzhu.Computing the Vertical Second Derivative and Upward Continuation of Gravity Anomaly by Spline Function Method［J］.Oil Geophysical Prospecting，1996，31（3）：415－422.（汪炳柱.用样条函数法求重力异常垂向二阶导数和向上延拓计算［J］.石油地球物理勘探，1996，31（3）：415－422.）

［6］ BLAKELY R J.Potential Theory in Gravity and Magnetic Applications［M］.Cambridge：Cambridge University Press，1995.

［7］ ZHANG Fengxu，ZHANG Fengqin，MENG Lingshun，et al.Magnetic Potential Spectrum Analysis and Calculating Method of Magnetic Anomaly Derivatives Based on Discrete Cosine Transform［J］.Chinese Journal of Geophysics，2007，50（1）：297－304.（张凤旭，张凤琴，孟令顺，等.基于离散余弦变换的磁位谱分析及磁异常导数计算方法［J］.地球物理学报，2007，50（1）：297－304.）

［8］ GUAN Zhining.Magnetic Field and Magnetic Exploration［M］.Beijing：Geological Publishing House，2005.（管志宁.地磁场与磁力勘探［M］.北京：地质出版社，2005.）

［9］ YAO Changli，HUANG Weining，GUAN Zhining.Fast Splines Conversion of Curved－surface Potential Field and Vertical Gradient Data into Horizontal－plane Data［J］.Oil Geophysical Prospecting，1997，32（2）：229－236.（姚长利，黄卫宁，管志宁.综合利用位场及其垂直梯度的快速样条曲化平方法［J］.石油地球物理勘探，1997，32（2）：229－236.）

A New Method to Calculate the Vertical Derivatives of Total Field Magnetic Anoma－

lyZHAI Guojun1，BIAN Guanglang1，2，HUANG Motao1
1.Naval Institute of Hydrographic Surveying and Charting，Tianjin300061，China；2.Department of Hydrography and Cartography，Dalian Naval Academy，Dalian 116018，China

The vertical derivative of the total field magnetic anomaly（TMA）plays an important role in the interpretation process of magnetic objects.Owing to the filter characteristics of vertical derivative operator in frequency domain，the conversion of vertical derivative is inherently unstable and any high－frequency noise presented in the surveying data gets strongly magnified in the transformed map in such a way to mask any useful signal.Based on the harmonic properties of the vertical integral and derivatives of TMA，an algorithm is presented to perform calculation of vertical derivatives using both frequency and space domain transformations.The effectiveness of the suggested techniques has been illustrated by synthetic sphere magnetic model whose TMA’s first and second vertical derivative expressions are deduced when the geomagnetic direction is different from magnetization direction.The conclusion indicates that the presented vertical derivatives calculation method provides better results and allows a lower degrading of the signal－to－noise ratio than the standard Fourier method.

total field magnetic anomaly；vertical derivative；space and frequency domain；Laplace equation；sphere magnetic model

ZHAI Guojun（1961—），male，PhD，senior engineer，PhD supervisor，majors in technology and application of bathymetry and marine geodesy.

1001－1595（2011）06－0671－08

P229

A

国家自然科学基金（40671161）

雷秀丽）

2010－12－28

2011－03－01

翟国君（1961—），男，博士，高级工程师，博士生导师，主要从事海底地形测量和海洋大地测量技术与应用研究。

E－mail：zhaigj＠163.com