翟國君，卞光浪，2，黃謨濤

1.海軍海洋測繪研究所，天津300061；2.大連艦艇學院海測工程系，遼寧大連116018

總強度磁異常各階垂向導數換算新方法

翟國君1，卞光浪1，2，黃謨濤1

1.海軍海洋測繪研究所，天津300061；2.大連艦艇學院海測工程系，遼寧大連116018

總強度磁異常垂向導數換算在磁性目標解釋推斷過程中具有重要意義。分析頻率域總強度磁異常各階垂向導數轉換因子的濾波特性，指出在求解總強度磁異常垂向導數時，常規(guī)傅里葉變換法會明顯放大觀測數據中高頻噪聲的影響，甚至可能淹沒掉真實信息。給出一個重要定理：總強度磁異常沿垂直方向的積分及各階偏導數均為準調和函數，在此基礎上，提出基于拉普拉斯方程的逐步求解法，聯合空間域和頻率域換算總強度磁異常沿垂直方向的各階導數。同時，采用球體磁場模型驗證所提方法的有效性，并推導地磁場方向和磁化強度方向不一致時，球體總強度磁異常沿垂直方向的一階和二階導數表達式。結論表明：利用樣條法計算的垂向導數結果精度明顯優(yōu)于常規(guī)傅里葉變換法計算結果，且具有較強的抗噪能力。

總強度磁異常；垂向導數；空間域和頻率域；拉普拉斯方程；球體磁場模型

1 引 言

磁測方法早期主要應用于地球物理和磁力勘探領域。近年來，隨著人類開發(fā)與利用海洋資源進程的加快，該方法以成本低、快速高效等突出優(yōu)點，在水下小尺度磁性目標精密探測方面?zhèn)涫荜P注，特別是沉船、鐵錨、海底電纜、未爆炸物以及水下潛航器的定位與識別方面更具明顯優(yōu)勢［1－2］。

磁異常導數被廣泛應用于磁性目標磁異常解釋，是壓制區(qū)域場、圈定局部場、分離疊加異常以及異常反演中常用的方法［3］。例如，利用解析信號法和歐拉反褶積法確定磁性目標平面位置，或者采用泰勒級數法實現磁場空間穩(wěn)健延拓，都需要以高精度磁異常水平和垂直方向的一階、二階甚至高階導數作為基礎數據?，F階段，獲取磁場垂向導數主要有兩種方式，其一是將磁力傳感器設計為磁陣列模式［4］，利用微分原理直接計算磁場的垂向導數，該模式對載體運行姿態(tài)要求較高，海洋上受風、流、浪等因素影響，實現起來較為困難；其二是將采集的剖面或平面磁場信息進行數據轉換處理，這種轉換既可在空間域內進行，也可在頻率域內進行。在空間域換算方面，文獻［5］嘗試將樣條函數插值法引入位場垂向導數解算中，但計算精度有限。目前普遍采用的磁異常垂向導數換算方法是傅里葉變換法［6］（簡記為“常規(guī)FFT法”），該方法理論較為完善，但其頻率響應函數相當于高頻放大器，對觀測數據中的噪聲具有顯著的放大作用，當觀測數據含有一定誤差時，計算精度較低，且存在較強的邊界效應。為避免復雜的復數運算，文獻［7］提出采用離散余弦變換（DCT）計算磁異常導數，該方法去除了原信號的相關性，在二度磁性目標磁異常垂向導數換算中應用效果較好，關于該方法能否適用于三度磁性目標磁異常垂向導數換算還有待于進一步研究。

針對總強度磁異常沿垂直方向導數的精確求解問題，筆者給出了關于總強度磁異常的一個重要定理，并綜合利用空間域和頻率域運算，提出一種求解總強度磁異常各階垂向導數的新方法。文中詳細研究了該方法的基本原理和技術措施，通過仿真球體磁場模型檢驗了方法的有效性，并將其計算結果與常規(guī)FFT法進行了分析比較。

2 常規(guī)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 調和函數二階水平方向導數計算方法

目前，計算調和函數沿水平方向的二階偏導數通常在頻率域內進行轉換（計算步驟與常規(guī)FFT法計算垂向導數類似），但該方法存在的不足是會顯著放大觀測數據中噪聲的影響，且具有較強的邊界效應。在數值計算中，三點二階中心差分法和三次樣條曲線函數法經常被用于計算函數在一維方向的二階導數，本文根據平面磁場數據的二維空間分布特點，將兩者引入并應用到調和函數二階水平方向偏導數計算中。下面以總強度磁異常為例，給出三點二階中心差分法和雙三次樣條曲線函數法計算其沿水平方向二階偏導數的具體方法原理。

4.1 三點二階中心差分法

假設x1，x2，…，xM為格網化數據沿x方向的等間距節(jié)點，且x1＜x2＜…＜xM，xi－xi－1＝Δx。則在（xi，y）節(jié)點處總強度磁異常沿x方向的二階偏導數可表示為

類似地，設y1，y2，…，yN為格網化數據沿y方向的等間距節(jié)點，且y1＜y2＜…＜yN，yiyi－1＝Δy。則在（x，yi）節(jié)點處總強度磁異常沿y方向的二階偏導數可表示為

值得注意的是，利用式（8）和式（9）求解總強度磁異常沿水平方向二階偏導數時，邊界處的二階導數需要附加條件，本文仿真計算時，考慮到僅對比梯度變化較大區(qū)域的計算精度，所以未計算邊界部分的二階偏導數。

4.2 雙三次樣條曲線函數法

雙三次樣條曲線函數是對x方向和y方向分別采用三次樣條曲線函數進行擬合［9］。以x方向測線為例，設有觀測數據 ｛x1，Bm（x1）｝，｛x2，Bm（x2）｝，…，｛xM，Bm（xM）｝，且x1＜x2＜…＜xM，整個區(qū)間［x1，xM］被節(jié)點分成M－1個小區(qū)間［x1，x2］、［x2，x3］、…、［xM－1，xM］，現要求構造函數S（x）滿足以下三個條件：

（1）插值條件

（2）在每個小區(qū)間［xi，xi＋1］（i＝1，2，…，M－1）上，S（x）是x的三次多項式，即

（3）在整個區(qū)間［x1，xM］上，S（x）具有連續(xù)的二階導數，在區(qū)間的每個節(jié)點xi（i＝1，2，…，M）上，函數值、一階導數和二階導數連續(xù)，即

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

利用雙三次樣條曲線函數計算水平方向二階偏導數時，需要給定兩個邊界約束條件，本文選擇自然三次樣條插值函數，即Sxx（x1）＝Sxx（xM）＝0，結合方程組（12）得到的Sxx（xi）即為xi（i＝2，3，…，M－1）節(jié)點處總強度磁異常沿x方向的二階偏導數。按照相同的方法，可進一步得到各測點處總強度磁異常沿y方向的二階偏導數。

5 仿真試驗分析

為檢驗提出的總強度磁異常各階垂向導數換算方法的有效性，仿真三個磁化球體磁異常數據進行分析。

仿真測區(qū)面積為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。

本文在前面給出了兩種垂向導數計算方法，分別是常規(guī)FFT法和基于拉普拉斯方程的逐步求解法。在利用后者求解垂向導數時，提供了兩種求解水平方向二階偏導數的算法，簡稱為“差分法”和“樣條法”。為充分驗證提出方法的適用性，這里分無噪聲和有噪聲兩種情況分別進行分析。

5.1 無噪聲情況

圖2和圖3是根據圖1（a）中的總強度磁異常數據，利用常規(guī)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可以看出，無噪聲情況下，三種換算方法得到的垂向一階偏導數和二階偏導數等值線圖與相應的理論等值線圖實現了很好的吻合，但常規(guī)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可以看到，三種垂向導數換算方法計算結果均具有較高的精度，其中，樣條法和常規(guī)FFT法計算結果更加接近于理論值。由于式（15）給出的精度評價指標屬于絕對誤差，其大小與各階垂向導數的振幅緊密相關，垂向導數的振幅越大，絕對誤差越大；反之，振幅越小，絕對誤差則越小。因此，該指標在一定程度上無法全面反映換算結果的優(yōu)劣，故進一步采用中誤差與理論值振幅之比作為評價指標

以樣條法為例，總強度磁異常垂向一階導數和二階導數理論值振幅分別為7.48nT／m和1nT／m2，計算得到的垂向一階、二階導數中誤差和振幅之比分別為0.11%與0.17%，驗證了該方法具有很高的垂向導數換算精度，且二階垂向導數換算結果的相對精度要稍低于一階垂向導數換算結果。因此，當實測資料中無噪聲或噪聲微小時，既可采用常規(guī)FFT法，也可采用樣條法計算總強度磁異常的垂向一階導數和二階導數。

考慮到利用泰勒級數法進行總強度磁異常空間延拓時，往往還需要獲取總強度磁異常的高階垂向導數，本文進一步對比分析了常規(guī)FFT法和樣條法的垂向三階、四階和五階導數換算結果。結果表明：常規(guī)FFT法計算結果存在較強的邊界效應，且隨著垂向導數階數的增加，其邊界效應更加明顯，到第五階垂向導數時，誤差已非常大；而樣條法雖然也存在一定的邊界效應，但其換算效果要明顯優(yōu)于常規(guī)FFT法結果。因此，實踐中建議采用樣條法計算總強度磁異常高階垂向導數。

5.2 含有噪聲情況

為更接近實際情況，進一步檢驗提出方法的穩(wěn)健性，在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%；常規(guī)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可以看出，當總強度磁異常觀測數據中含有噪聲時，同一噪聲水平情況下，樣條法計算結果精度最高，差分法次之，常規(guī)FFT法精度最低；而對于同一種方法，垂向二階導數計算結果相對精度要低于一階導數計算結果相對精度。同時，隨著噪聲水平的逐漸增大，各階垂向導數換算結果精度也不斷降低。

表3 不同噪聲水平下三種方法換算得到的垂向導數中誤差Tab.3 Mean square error of three methods with different noise level

6 結 論

（1）分析了常規(guī)FFT法在換算總強度磁異常各階垂向導數過程中存在的不足；給出了一個重要定理：在無源空間內，總強度磁異常沿垂直方向的積分和各階偏導數同樣為準調和函數，在此基礎上，提出了一種新的總強度磁異常各階垂向導數換算方法。

（2）圍繞調和函數的二階水平方向偏導數求解問題，給出了三點二階中心差分法和雙三次樣條曲線函數法。通過仿真試驗對三點二階中心差分法和雙三次樣條曲線函數法應用效果進行了比較分析，結果表明，利用樣條法換算得到的各階垂向導數精度較高，且具有較強的抗噪聲能力。

（3）當總強度磁異常觀測數據不含噪聲時，常規(guī)FFT法、差分法和樣條法換算的總強度磁異常垂向一階導數與二階導數結果精度大體相當，但隨著其垂向導數階數的增加，常規(guī)FFT法邊界效應現象較其他兩種方法嚴重得多。當總強度磁異常觀測數據含有噪聲時，隨著噪聲水平的不斷提高，垂向導數換算結果精度逐漸降低。但總體而言，樣條法垂向導數換算結果精度要明顯優(yōu)于常規(guī)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.（姚長利，黃衛(wèi)寧，管志寧.綜合利用位場及其垂直梯度的快速樣條曲化平方法［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