總強度磁異常各階垂向?qū)?shù)換算新方法

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

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

總強度磁異常各階垂向?qū)?shù)換算新方法

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

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

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

總強度磁異常；垂向?qū)?shù)；空間域和頻率域；拉普拉斯方程；球體磁場模型

1 引 言

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

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

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

2 常規(guī)FFT法

由傅里葉變換，三度磁性目標(biāo)總強度磁異常Bm的一階垂向?qū)?shù)Bmz與其頻譜關(guān)系［8］為

根據(jù)位場理論，總強度磁異常近似滿足拉普拉斯方程，利用外部Dirichlet問題條件，其空間延拓公式可表示為

式中，F(xiàn)Bm（u，v，0）為z＝0平面上總強度磁異常Bm（x，y，0）的頻譜。

將式（2）對z方向求偏導(dǎo)數(shù)，得

由于同一函數(shù)有相同的傅里葉表達式，對比式（1）和式（3）兩端，得

當(dāng)z＝0時，式（4）即為平面上總強度磁異常頻譜向其一階垂向?qū)?shù)頻譜轉(zhuǎn)換的關(guān)系式，同理，還可推導(dǎo)出總強度磁異常沿垂向的第k階導(dǎo)數(shù)頻譜為這里稱為總強度磁異常與其第k階垂向?qū)?shù)間的頻率域轉(zhuǎn)換因子（也稱頻率響應(yīng)函數(shù)）。計算時，首先將平面上總強度磁異常格網(wǎng)數(shù)據(jù)利用快速傅里葉變換轉(zhuǎn)換為頻譜FBm（u，v，0），再乘以垂向?qū)?shù)頻率域轉(zhuǎn)換因子獲取第k階垂向?qū)?shù)的頻譜，最后通過反傅里葉變換即可得到總強度磁異常各階垂向?qū)?shù)。

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

為減弱高頻噪聲對總強度磁異常垂向?qū)?shù)計算的影響，這里提出基于拉普拉斯方程的逐步求解法。在求解過程中，將運用總強度磁異常的某些特殊性質(zhì)，限于篇幅，略去推導(dǎo)過程，給出如下定理。

定理：在無源空間內(nèi)，總強度磁異常沿垂直方向的積分及任意階偏導(dǎo)數(shù)均為準(zhǔn)調(diào)和函數(shù)。

結(jié)合上述定理，總強度磁異常各階垂向?qū)?shù)的求解步驟如下：

（1）對z＝0平面上總強度磁異常數(shù)據(jù)進行傅里葉變換，設(shè)其在頻率域頻譜為FBm（u，v，0），由式（4），F(xiàn)Bm（u，v，0）乘以頻率域轉(zhuǎn)換因子即為總強度磁異常沿垂直方向積分的頻譜，再利用反傅里葉變換得到積分值J（x，y，0）；

（2）計算J（x，y，0）沿x方向和y方向的二階偏導(dǎo)數(shù)Jxx（x，y，0）與Jyy（x，y，0），由定理知，總強度磁異常沿垂直方向的積分是準(zhǔn)調(diào)和函數(shù)，滿足拉普拉斯方程，則總強度磁異常沿垂直方向的一階偏導(dǎo)數(shù)為

（3）計算Bm在水平方向的兩個二階偏導(dǎo)數(shù)與而總強度磁異常自身也為準(zhǔn)調(diào)和函數(shù)，同樣利用拉普拉斯方程，得總強度磁異常沿垂直方向的二階偏導(dǎo)數(shù)為

（4）由于總強度磁異常沿垂直方向的任意階偏導(dǎo)數(shù)都是準(zhǔn)調(diào)和函數(shù)，得第k（k＝3，4，…，n）階垂直導(dǎo)數(shù)為

式中，Bmkz為總強度磁異常沿垂直方向的第k階偏導(dǎo)數(shù)和分別是Bm（k－2）z沿x方向和y方向的二階偏導(dǎo)數(shù)。

4 調(diào)和函數(shù)二階水平方向?qū)?shù)計算方法

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

4.1 三點二階中心差分法

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

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

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

4.2 雙三次樣條曲線函數(shù)法

雙三次樣條曲線函數(shù)是對x方向和y方向分別采用三次樣條曲線函數(shù)進行擬合［9］。以x方向測線為例，設(shè)有觀測數(shù)據(jù) ｛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］，現(xiàn)要求構(gòu)造函數(shù)S（x）滿足以下三個條件：

（1）插值條件

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

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

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

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

5 仿真試驗分析

為檢驗提出的總強度磁異常各階垂向?qū)?shù)換算方法的有效性，仿真三個磁化球體磁異常數(shù)據(jù)進行分析。

仿真測區(qū)面積為295m×295m，測線和測點間距均為5m。地磁場傾角和偏角分別為45°與5°，為檢驗方法的普遍適用性，這里設(shè)球體的磁化強度方向和地磁場方向不一致，球體空間位置及磁性參數(shù)分別如表1所示。

表1 仿真球體模型相關(guān)參數(shù)Tab.1 Parameters of synthetic sphere model

為定量分析總強度磁異常的垂向一階導(dǎo)數(shù)、二階導(dǎo)數(shù)計算結(jié)果精度，在文獻［8］給出的球體總強度磁異常表達式基礎(chǔ)上，筆者進一步推導(dǎo)磁化強度方向和地磁場方向不一致時，磁化球體總強度磁異常垂向一階導(dǎo)數(shù)Bmz和二階導(dǎo)數(shù)Bmzz解析表達式，分別如式（13）與式（14）所示

式中，I0和A0分別為地磁場傾角和偏角；Bxz、Byz、Bzz和Bxzz、Byzz、Bzzz分別是Bv（v＝x，y，z）沿垂直方向的一階偏導(dǎo)數(shù)與二階偏導(dǎo)數(shù)，其中，

根據(jù)球體總強度磁異常以及式（13）和式（14），仿真得到z＝0m平面上總強度磁異常及其垂向一階和二階偏導(dǎo)數(shù)等值線如圖1所示。

圖1 總強度磁異常及其垂向一階和二階偏導(dǎo)數(shù)理論等值線圖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）中總強度磁異常垂向一階導(dǎo)數(shù)最大值為5.30nT／m，最小值為－2.18nT／m，平均值為0.001 8nT／m；圖1（c）中總強度磁異常垂向二階導(dǎo)數(shù)最大值為0.70nT／m2，最小值為－0.30nT／m2，平均值為－0.000 1nT／m2。

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

5.1 無噪聲情況

圖2和圖3是根據(jù)圖1（a）中的總強度磁異常數(shù)據(jù)，利用常規(guī)FFT法、差分法和樣條法換算得到總強度磁異常沿垂直方向的一階偏導(dǎo)數(shù)和二階偏導(dǎo)數(shù)等值線圖。

圖2 三種方法計算得到的總強度磁異常沿垂直方向的一階偏導(dǎo)數(shù)等值線圖（等值線間隔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 三種方法計算得到的總強度磁異常沿垂直方向的二階偏導(dǎo)數(shù)等值線圖（等值線間隔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可以看出，無噪聲情況下，三種換算方法得到的垂向一階偏導(dǎo)數(shù)和二階偏導(dǎo)數(shù)等值線圖與相應(yīng)的理論等值線圖實現(xiàn)了很好的吻合，但常規(guī)FFT法隨著導(dǎo)數(shù)階數(shù)的增加，邊界效應(yīng)越來越明顯。為定量計算總強度磁異常各階垂向?qū)?shù)換算精度，并排除邊界效應(yīng)對統(tǒng)計結(jié)果的影響，選取上述各圖中虛線框D域內(nèi)梯度變化較大的格網(wǎng)點數(shù)據(jù)進行分析，即D＝｛（x，y）｜25m≤x≤270m，25m≤y≤270m｝，采用的精度評價指標(biāo)公式為式（15）。

式中，＾Bmkz（i，j）和Bmkz（i，j）（k為正整數(shù)）分別表示格網(wǎng)點上總強度磁異常第k階垂向?qū)?shù)計算值與理論值；Nx和Ny分別為D域內(nèi)x方向和y方向采樣點數(shù)；ρk為第k階垂向?qū)?shù)計算結(jié)果中誤差。

三種方法換算結(jié)果統(tǒng)計情況如表2所示。

表2 三種方法換算的垂向一階和二階導(dǎo)數(shù)結(jié)果統(tǒng)計Tab.2 Statistics of the 1st and 2ed vertical derivative results with three methods

從表2可以看到，三種垂向?qū)?shù)換算方法計算結(jié)果均具有較高的精度，其中，樣條法和常規(guī)FFT法計算結(jié)果更加接近于理論值。由于式（15）給出的精度評價指標(biāo)屬于絕對誤差，其大小與各階垂向?qū)?shù)的振幅緊密相關(guān)，垂向?qū)?shù)的振幅越大，絕對誤差越大；反之，振幅越小，絕對誤差則越小。因此，該指標(biāo)在一定程度上無法全面反映換算結(jié)果的優(yōu)劣，故進一步采用中誤差與理論值振幅之比作為評價指標(biāo)

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

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

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是用含有高斯噪聲的總強度磁異常數(shù)據(jù)換算得到其沿垂直方向的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)等值線圖。

如圖5和圖6所示，當(dāng)總強度磁異常數(shù)據(jù)中加入2%高斯白噪聲后，垂向一階導(dǎo)數(shù)和二階導(dǎo)數(shù)計算結(jié)果精度要明顯低于無噪聲情況下的結(jié)果，表明噪聲對垂向?qū)?shù)換算結(jié)果存在較大影響。在三種方法中，樣條法計算結(jié)果精度最高，根據(jù)D域內(nèi)格網(wǎng)點數(shù)據(jù)，利用式（15）和式（16）計算得到的ρ1＝0.07nT／m，η1＝0.93%，ρ2＝0.05nT／m2，η2＝5%；差分法計算結(jié)果精度次之，對應(yīng)ρ1＝0.13nT／m，η1＝1.7%，ρ2＝0.11nT／m2，η2＝11%；常規(guī)FFT法計算結(jié)果精度最低，其中ρ1＝0.21nT／m，η1＝2.8%，ρ2＝0.13nT／m2，η2＝13%。從相對精度指標(biāo)變化情況來看，垂向二階導(dǎo)數(shù)較一階導(dǎo)數(shù)受噪聲的影響更加嚴(yán)重。

圖5 三種方法計算得到的總強度磁異常沿垂直方向的一階偏導(dǎo)數(shù)等值線圖（等值線間隔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%等不同噪聲水平情況下，三種垂向?qū)?shù)方法計算得到的垂向一階、二階導(dǎo)數(shù)結(jié)果精度。

圖6 三種方法計算得到的總強度磁異常沿垂直方向的二階偏導(dǎo)數(shù)等值線圖（等值線間隔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可以看出，當(dāng)總強度磁異常觀測數(shù)據(jù)中含有噪聲時，同一噪聲水平情況下，樣條法計算結(jié)果精度最高，差分法次之，常規(guī)FFT法精度最低；而對于同一種方法，垂向二階導(dǎo)數(shù)計算結(jié)果相對精度要低于一階導(dǎo)數(shù)計算結(jié)果相對精度。同時，隨著噪聲水平的逐漸增大，各階垂向?qū)?shù)換算結(jié)果精度也不斷降低。

表3 不同噪聲水平下三種方法換算得到的垂向?qū)?shù)中誤差Tab.3 Mean square error of three methods with different noise level

6 結(jié) 論

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

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

（3）當(dāng)總強度磁異常觀測數(shù)據(jù)不含噪聲時，常規(guī)FFT法、差分法和樣條法換算的總強度磁異常垂向一階導(dǎo)數(shù)與二階導(dǎo)數(shù)結(jié)果精度大體相當(dāng)，但隨著其垂向?qū)?shù)階數(shù)的增加，常規(guī)FFT法邊界效應(yīng)現(xiàn)象較其他兩種方法嚴(yán)重得多。當(dāng)總強度磁異常觀測數(shù)據(jù)含有噪聲時，隨著噪聲水平的不斷提高，垂向?qū)?shù)換算結(jié)果精度逐漸降低。但總體而言，樣條法垂向?qū)?shù)換算結(jié)果精度要明顯優(yōu)于常規(guī)FFT法換算結(jié)果。

［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.（于波，翟國君，劉雁春，等.利用航磁數(shù)據(jù)向下延拓得到海面磁場的方法［J］.測繪學(xué)報，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.（汪炳柱.用樣條函數(shù)法求重力異常垂向二階導(dǎo)數(shù)和向上延拓計算［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.（張鳳旭，張鳳琴，孟令順，等.基于離散余弦變換的磁位譜分析及磁異常導(dǎo)數(shù)計算方法［J］.地球物理學(xué)報，2007，50（1）：297－304.）

［8］ GUAN Zhining.Magnetic Field and Magnetic Exploration［M］.Beijing：Geological Publishing House，2005.（管志寧.地磁場與磁力勘探［M］.北京：地質(zhì)出版社，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

國家自然科學(xué)基金（40671161）

雷秀麗）

2010－12－28

2011－03－01

翟國君（1961—），男，博士，高級工程師，博士生導(dǎo)師，主要從事海底地形測量和海洋大地測量技術(shù)與應(yīng)用研究。

E－mail：zhaigj＠163.com