改進的局部波數法及其在磁場數據解釋中的應用

李麗麗，杜曉娟，馬國慶

吉林大學地球探測科學與技術學院，長春 130021

改進的局部波數法及其在磁場數據解釋中的應用

李麗麗，杜曉娟，馬國慶

吉林大學地球探測科學與技術學院，長春 130021

局部波數法是一種進行磁場數據解釋的常用方法。現有的局部波數法在進行反演時往往需要計算局部波數的導數，會明顯地增大噪聲的干擾，為解釋結果帶來誤差。提出3種利用磁異常在不同位置或不同高度上局部波數的簡單組合來進行場源體深度及構造指數的計算方法，不需要計算局部波數的導數，降低了噪聲的干擾，增強了反演結果的穩定性。通過理論模型試驗，證明改過的局部波數法在有無噪聲的情況下均能很好地完成異常的反演，其反演結果與理論值之間的差距小于理論值的5%。將其應用于四川某地區磁異常的反演中，其結果與解析信號的歐拉反褶積法的反演結果吻合。

磁場；深度；構造指數；導數；局部波數

0 前言

磁數據解釋的主要任務是確定場源體的深度及類型（構造指數）。為此人們提出了很多方法來完成這一任務，有維納反褶積法［1］、解析信號法［2］、歐拉反褶積法［3－4］、最小二乘法［5－6］、場源參數成像［7］及人工神經網絡法［8］等，其中場源參數法是近期應用最為廣泛的一種方法。該方法利用局部波數快速地進行場源體反演；但需人為給定研究區場源體的構造指數，而一個地區場源體的構造指數往往是未知的。為了解決這一問題，人們對該方法進行了改進［9－11］，使其能夠同時完成深度及構造指數的反演。然而，改進后的方法需要計算局部波數的一階甚至更高階導數，對數據精度要求較高，且會增大噪聲干擾；當測量點距較大時，高階導數的計算非常不穩定。Pilkington等［12］在2006年證明局部波數為解析信號［13］的導數，因此局部波數在二維情況下具有不受磁化方向干擾的優點。Salem等［14－15］利用原始局部波數和相位轉換后的局部波數完成了異常體位置的反演，并于2008年試驗了該方法在三維情況下的應用效果。Keating［16］利用局部波數及其垂直導數對異常體的深度及構造指數進行反演，但是該方法會增大噪聲的干擾。

筆者對局部波數法進行改進，采用不同位置和不同高度磁異常的局部波數直接計算場源體的深度和構造指數，不需要計算局部波數的導數，且降低了噪聲干擾。

1 改進的局部波數法

對于磁場T而言，其局部波數kx被定義為局部相位θ的水平方向的導數［17］：

局部相位θ被定義為

其中，?T／?x和?T／?z分別表示異常T的水平與垂直導數。對式（1）進行微分后得到

為異常T的解析信號。Smith等［11］推導出水平位置為x0、埋深為z0的磁性體在點（x，z）處引起的局部波數kx的表達式為

其中，n為與場源體類型對應的構造指數。Stavrev等［18］推導了磁源體類型與構造指數之間的對應關系，其中接觸帶的構造指數為0，巖墻為1，水平圓柱體為2，點源為3。

異常體的水平位置x0可以通過解析信號和局部波數的極大值來確定。局部波數kx在水平位置x＝x0、垂直位置z＝0處的值為

局部波數kx在水平位置x＝x0＋c，垂直位置z＝0處的值為

將式（6）與式（7）相除后可以得到

利用方程（8）可獲得場源體的深度，然后利用方程（6）計算場源體的構造指數：

經過以上的推導可以看出，利用局部波數上2個不同水平位置的點依據方程（8）和（9）就可以獲得異常體的深度和構造指數。在反演過程中，計算點應選取在場源體中心位置附近；因為在該范圍內異常的信噪比較大，得到的反演結果更準確。筆者給出了一種利用不同高度的局部波數kx進行磁異常反演的方法，為了使輸出結果更加穩定，采用上延的方法來獲得不同高度上的局部波數。局部波數kx在水平位置x＝x0、垂直位置z＝－d處的值為

將（6）式與（10）式相除后得到

分別計算方程（6）與（10）的倒數然后相減得到

從上述公式中可以看出，利用2個不同高度的局部波數kx，由（11）和（12）式可以得到異常體的深度和構造指數。除上述2種方法外，筆者還提出了另外一種利用相位轉換后的局部波數進行磁異常解釋的方法。

Salem等［13］定義相位轉換后的局部波數kz的表達式：

并推導出水平位置在x0、埋深為z0磁性體的局部波數kz在點（x，z）處的表達式為

局部波數kz在水平位置x＝x0＋b、垂直位置z＝0處的表達式為

局部波數kz在水平位置x＝x0＋c、垂直位置z＝0處的表達式為

將方程（15）與（16）相除后得到

利用（17）式可以得到場源體的深度，然后根據方程（15）計算異常體的構造指數：

利用方程（17）和（18）就可以計算出場源體的深度和構造指數，考慮反演結果的穩定性，b一般選取為1。

經過以上的推導，筆者給出進行磁性體深度和構造指數計算的3種方法，這3種方法不需要局部波數的導數參與計算，有效地提高了結果的穩定性。在利用上述方法進行反演時，為了提高結果的準確性，計算點應選取在局部波數kz兩個極值點之間；因為兩個極值點對應于異常體的邊界，當測點在異常體上方時其反演結果是最準確的［4］。

2 模型試驗

為了試驗方法的適用性，設計了如下模型：在長為100km的剖面中點位置處存在一個頂面埋深為4km的巖脈，頂面寬度為10km，磁化強度為20A／km，地磁傾角為60°，磁偏角為0°，測量點距為1 km。利用筆者提出的3種方法對磁異常進行反演，其反演結果為2個計算點之間的水平距離（c）和上延高度（d）的函數（圖1），并給出了反演結果與理論值之間的標準偏差（sd）。

從反演結果中可以看出：當水平距離小于14 km、上延高度小于7km時，本文3種方法反演得到的構造指數均與理論值之間的偏差小于0.01，反演得到的深度與真實深度之差小于0.1km，其誤差小于理論值的2%，所以本文方法能很好地完成異常的反演，但隨著水平距離與上延高度的增大，標準偏差增大，反演結果的精度降低。

在進行實際數據處理時，噪聲是必須考慮的一個因素。為了試驗方法存在噪聲時的應用效果，在圖1a的原始異常中加入均值為0、均方差為1nT的隨機噪聲，采用筆者提出的3種方法對異常進行反演（圖2）。

從圖2中可以看出，在存在噪聲時本文方法的反演結果有一定的波動。筆者采用如下準則來獲得最終結果：1）選取水平距離為0～10km的所有深度解，計算其平均值；2）將與平均值差異大于平均值25%（根據實際情況設定）的解去除，重新計算平均值；3）以新平均值為最后深度，計算構造指數。依據該準則獲得的反演結果如表1所示。

表1 不同方法的反演結果Table 1 Inversion results of different methods

從表1中可以看出，在存在噪聲的情況下本文方法的反演結果與理論值之間差距也較小。

為了獲得更加準確的結果，可采用向上延拓技術先對異常進行處理，然后再進行反演。筆者對異常進行向上延拓處理，當延拓高度為2km時局部波數已足夠圓滑，即局部波數kx為單峰值（圖3）。

從反演結果（圖3）可以看出，與圖2相比，向上延拓后異常反演結果的不規則波動消失，計算反演結果的平均值，式（8）和（9）反演得到的深度為4.01 km，構造指數為1.01；式（11）和式（12）反演得到的深度為4.02km，構造指數為1.01；式（17）和式（18）反演得到的深度為4.01km，構造指數為1.01，因此對異常進行向上延拓處理后反演結果的精度得到提高，其值更加準確。因此，實際數據處理前對數據進行向上延拓處理后會提高反演的精度。

3 實際數據處理

利用筆者提出的方法解釋四川某段實測航磁異常（圖4）。數據采樣間距為100m。將異常向上延拓200m，根據解析信號和局部波數的最大值得出異常體的水平位置在距起始點2 950m處。

從反演結果可以看出：利用不同方法反演得到的構造指數均接近1，表明地質體為一巖墻，反演得到深度的平均值為543m。

圖1 模型理論異常及不同方法的反演結果Fig.1 Theoretical anomaly and inversion results of different methods

該地區沒有其他的數據來驗證反演結果的正確性，將筆者的反演結果與解析信號歐拉反褶積法［19］的反演結果進行對比。解析信號歐拉反褶積法的反演結果顯示異常體的位置在（2 945±7）m，深度為（554.6±18.2）m，構造指數大小為0.98±0.03。這與本文方法的反演結果相一致，表明本文方法反演結果是準確的。

圖2 加噪后異常及不同方法的反演結果Fig.2 Noise－corrupted anomaly and the inversion results of different methods

圖3 濾波后異常及不同方法的反演結果Fig.3 Smoothed anomaly and the inversion results of the methods

圖4 實際航磁異常的反演結果Fig.4 Inversion results of real aeromagnetic anomaly

4 結論

1）筆者提出3種改進的局部波數法進行磁性體反演方法，這些方法僅利用不同水平位置和不同高度的局部波數來進行深度和構造指數的計算，不需要計算局部波數的導數，降低了對數據精度的要求。

2）通過理論模型試驗證明，新方法在有無噪聲情況下，均能準確地計算出地質體的深度和構造指數。

3）在有噪聲情況下，采用向上延拓技術處理數據后可獲得更加精確的結果，為實際數據解釋提供了理論依據。

4）將新方法應用于四川實測剖面磁數據的解釋中，其反演結果與解析信號的歐拉反褶積法的結果相吻合，表明本文方法在實際中也有較好的應用效果。

（References）：

［1］ Hartman R R，Teskey D J，Friedberg J L.A System for Rapid Digital Aeromagnetic Interpretation［J］.Geophysics，1971，36：891－918.

［2］ Nabighian M N.The Analytic Signal of Two Dimensional Magnetic Bodies with Polygonal Cross Section：Its Properties and Use for Automated Anomaly Interpretation［J］.Geophysics，1972，37：507－517.

［3］ Thompson D T.“EULDPH”：A New Technique for Making Computer－Assisted Depth Estimates fromMagnetic Data［J］.Geophysics，1982，47：31－37.

［4］ Reid A B，Allsop J M，Granser H，et al.Magnetic Interpretation in Three Dimensions Using Euler Deconvolution［J］.Geophysics，1990，55：80－91.

［5］ Abdelrahman E M.Discussion on“A Least－Squares Approach to Depth Determination from Gravity Data”by O P Gupta［J］.Geophysics，1990，55：376－378.

［6］ Salem A，Ravat D，Mushayandebvu M F，et al.Linearized Least－Squares Method for Interpretation of Potential Field Data from Sources of Simple Geometry［J］.Geophysics，2004，69：783－788.

［7］ Thurston J B，Smith R S.Automatic Conversion of Magnetic Data to Depth，Dip，and Susceptibility Contrast Using the SPI Method［J］.Geophysics，1997，62：807－813.

［8］ Elawadi E，Salem A，Ushijima K.Detection of Cavities from Gravity Data Using A Neural Network［J］.Exploration Geophysics，2001，32：75－79.

［9］ Thurston J B，Smith R S，Guillon J.A Multi－Model Method for Depth Estimation from Magnetic Data［J］.Geophysics，2002，67：555－561.

［10］ Salem A，Smith R S.Depth and Structural Index from the Normalized Local Wavenumber of 2DMagnetic Anomalies［J］.Geophysical Prospecting，2005，51：83－89.

［11］ Smith R S，Thurston J B，Dai T，et al.ISPI：The Improved Source Parameter Imaging Method［J］.Geophysical Prospecting，1998，46：141－151.

［12］ Pilkington M，Keating P.The Relationship Between Local Wavenumber and Analytic Signal in Magnetic Interpretation［J］.Geophysics，2006，71（1）：L1－L3.

［13］ 馬國慶，杜曉娟，李麗麗.利用水平與垂直導數的相關系數進行位場數據的邊界識別［J］.吉林大學學報：地球科學版，2011，41（增刊）：345－348.

Ma Guo－qing，Du Xiao－juan，Li Li－li.Edge Detection of Potential Field Using Correlation Coefficient of Horizontal and Vertical Derivatives［J］.Journal of Jilin University：Earth Science Edition，2011，41（Sup.）：345－348.

［14］ Salem A，Ravat D，Smith R，et al.Interpretation of Magnetic Data Using an Enhanced Local Wavenumber（ELW）Method［J］.Geophysics，2005，70（2）：L7－L12.

［15］ Salem A，Williams S，Fairhead D.Interpretation of Magnetic Data Using Tilt－Angle Derivatives［J］.Geophysics，2008，73：L1－L10.

［16］ Keating P.Improved Use of the Local Wavenumber in Potential－Field Interpretation［J］.Geophysics，2009，74（6）：L75－L85.

［17］ Bracewell R N.The Fourier Transform and Its Application［M］.New York：McGraw Hill Book Co，1965.

［18］ Stavrev P，Reid A.Degrees of Homogeneity of Potential Fields and Structural Indicates of Euler Deconvolution［J］.Geophysics，2007，72（1）：L1－L12.

［19］ Keating P，Pilkington M.Euler Deconvolution of the Analytic Signal and Its Application to Magnetic Interpretation［J］.Geophysical Prospecting，2004，52：165－182.

Improved Local Wavenumber Methods in the Interpretation of Magnetic Fields

Li Li－li，Du Xiao－juan，Ma Guo－qing

College of GeoExploration Science and Technology，Jilin University，Changchun 130026，China

Local wavenumber method is a widely used method in the interpretation of magnetic fields，the previous local wavenumber methods need to compute the derivatives of local wavenumber in the inversion，which increase the effect of noise obviously，and bring the errors to inversion results.The authors proposed using the simple combinations of different locations and heights of local wavenumber of magnetic anomaly to compute the depth and structural index of the causative source，and these methods do not need to compute the derivatives of the local wavenumber，which reduce the effect of the noise，and increase the stability of the inversion results.The new methods are demonstrated on theoretical magnetic field data without and with random noise，and the difference between the inversion results and theoretical values is less than 5%of the theoretical values.The practical application of the new methods is tested on real magnetic anomaly from Sichuan basin，Southwest China，the inversion results are consistent with the results computed by the Euler deconvolution of analytic signal.The proposed local wavenumber methods can finish the inversion of magnetic data effectively.

magnetic fields；depth；structural index；derivatives；local wavenumber

book=2012,ebook=558

P631.2

A

1671－5888（2012） 04－1179－07

2011－10－23

國家深部探測技術與實驗研究專項項目（SinoProbe－09－01）

李麗麗（1983－），女，博士研究生，主要從事地球物理數據處理與解釋方面的研究，E－mail：lilili10＠mails.jlu.edu.cn。