重磁異常解釋的歸一化局部波數(shù)法

馬國(guó)慶，黃大年，李麗麗，于平

吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130021

1 引言

重磁異常解釋的主要任務(wù)之一是確定異常體的空間賦存狀態(tài).局部波數(shù)法是近幾年應(yīng)用較為廣泛的一種自動(dòng)解釋方法，該方法利用局部波數(shù)可快速地完成場(chǎng)源體的反演（Thurston and Smith，1997），但該方法在進(jìn)行場(chǎng)源體深度反演時(shí)需要已知場(chǎng)源體的構(gòu)造指數(shù)，但未知地區(qū)場(chǎng)源體的構(gòu)造指數(shù)往往是未知的，構(gòu)造指數(shù)的錯(cuò)誤選取會(huì)為反演結(jié)果帶來(lái)較大的誤差.后來(lái)人們對(duì)該方法進(jìn)行改進(jìn)，使其能夠同時(shí)完成深度及構(gòu)造指數(shù)的反演（Smith etal.，1998；Salem and Smith，2005），然而改進(jìn)后方法需要計(jì)算異常的三階甚至更高階導(dǎo)數(shù)，對(duì)數(shù)據(jù)精度要求較高，且當(dāng)測(cè)量點(diǎn)距較大時(shí)，高階導(dǎo)數(shù)的計(jì)算是非常不穩(wěn)定的.Salem等（2005）利用原始局部波數(shù)和相位轉(zhuǎn)換后局部波數(shù)進(jìn)行異常體的反演，并于2008年（Salem etal.，2008）試驗(yàn)了該方法在三維情況下的應(yīng)用效果.Keating（2009）利用局部波數(shù)及其垂直導(dǎo)數(shù)對(duì)異常體的深度及構(gòu)造指數(shù)進(jìn)行反演，但是該方法會(huì)增大噪聲的干擾.Ma（2013）利用不同水平位置與不同高度局部波數(shù)曲線(xiàn)的組合進(jìn)行異常的解釋?zhuān)@得了良好的效果.歸一化總梯度法是一種利用高精度重力異常來(lái)確定異常源分布的方法，該方法不需要附加條件，計(jì)算簡(jiǎn)便，并能半定量地確定地質(zhì)體的位置而被廣泛用于研究含油氣構(gòu)造或檢測(cè)油氣（肖一鳴和張林詳，1984；沈慶夏等，2010），但該方法的結(jié)果由展開(kāi)式的項(xiàng)次來(lái)決定，為異常的解釋工作帶來(lái)了不確定性（孟平等，2003），為此提出采用Fourier變換和Hilbert變換來(lái)完成歸一化總梯度的計(jì)算，避免了展開(kāi)項(xiàng)次的約束（曾華霖等，1999；肖鵬飛等，2006）.

本文提出歸一化局部波數(shù)法進(jìn)行重磁異常的解釋?zhuān)⑻岢龆喾N不同的歸一化函數(shù).通過(guò)理論模型和實(shí)際數(shù)據(jù)驗(yàn)證方法的可行性和有效性，并給出了不同歸一化方式的應(yīng)用效果.

2 歸一化局部波數(shù)法

對(duì)于異常T而言，局部波數(shù)kx被定義為

局部相位θ的水平方向的導(dǎo)數(shù)：

局部相位θ被定義為其中，?T／?x和?T／?z分別表示異常T的水平與垂直導(dǎo)數(shù).對(duì)式（1）進(jìn)行微分后得到

Salem等（2005）定義相位轉(zhuǎn)換后局部波數(shù)kz的表達(dá)式

歸一化局部波數(shù)的表達(dá)式為：

其中，N（z）為歸一化函數(shù)，主要包含以下幾種方式

其中：

局部波數(shù)為異常的二階導(dǎo)數(shù)，因此歸一化局部波數(shù)法相對(duì)歸一化總梯度法受背景異常干擾更小.為了降低下延過(guò)程的不穩(wěn)定性及噪聲的干擾，采用ISVD（Integrated Second Vertical Derivative）法來(lái)完成局部波數(shù)的延拓工作（Fedi and Florio，2002；Rapolla etal.，2002），具體表達(dá)式為：

其中，m為泰勒展開(kāi)式的階數(shù).為了減小噪聲的干擾，引入Laplace方程來(lái)完成二階垂直導(dǎo)數(shù)的計(jì)算：對(duì)于其他階垂直導(dǎo)數(shù)的計(jì)算采用如下的方式來(lái)進(jìn)行（圖1）.

垂直導(dǎo)數(shù)的計(jì)算采用兩個(gè)水平導(dǎo)數(shù)之和來(lái)完成，水平導(dǎo)數(shù)采用空間域方法來(lái)計(jì)算，具有受噪聲干擾小和計(jì)算穩(wěn)定的優(yōu)勢(shì)，因此該方法能穩(wěn)定地完成異常的延拓工作.由于導(dǎo)數(shù)隨埋深衰減速度更快，因此足夠多的階數(shù)參與計(jì)算就會(huì)得到準(zhǔn)確的結(jié)果，一般展開(kāi)到7～10階就可以滿(mǎn)足精度要求.

圖1 不同階垂直導(dǎo)數(shù)的計(jì)算方法Fig.1 Computation method of vertical derivatives with different orders

3 理論模型試驗(yàn)

水平位置為100m，埋深為20m，半徑為5m，與圍巖密度差為1g·cm-3的圓柱體所引起的重力異常如圖2a所示.分別利用Fourier變換法和ISVD法將異常向下延拓4m，ISVD法計(jì)算時(shí)階數(shù)m為7.

從圖2b中可以看出，ISVD方法能較穩(wěn)定地完成異常的下延工作，且與理論異常之間的差距較小.Fourier變換計(jì)算結(jié)果穩(wěn)定性較差，異常出現(xiàn)不規(guī)則地波動(dòng)，為此需對(duì)數(shù)據(jù)進(jìn)行低通濾波處理，濾波后數(shù)據(jù)較圓滑，但與原始數(shù)據(jù)的差距較大.

利用歸一化局部波數(shù)法對(duì)該異常進(jìn)行反演（圖3）.

從圖3中可以看出，歸一化局部波數(shù)法能準(zhǔn)確地獲得地質(zhì)體的位置與深度信息，且中值、幾何平均和調(diào)和平均歸一化局部波數(shù)相對(duì)算術(shù)平均歸一化局部波數(shù)具有更高的分辨率，反演得到的結(jié)果更加收斂.

先給出歸一化局部波數(shù)法對(duì)于多個(gè)地質(zhì)體產(chǎn)生異常的應(yīng)用效果.地下水平位置50、100和150m存在埋深均為10m的垂直薄板狀體，其寬度為3m，磁化強(qiáng)度為20A／m，在磁傾角為70°時(shí)板狀體引起的磁異常如圖4a所示，利用歸一化局部波數(shù)法對(duì)磁異常進(jìn)行反演.

圖2 （a）圓柱體重力異常；（b）下延4m后重力異常Fig.2 （a）Synthetic gravity anomalies of a horizontal cylinder；（b）Gravity anomalies after downward continuation by 4m

圖3 （a）單個(gè)圓柱體產(chǎn)生的重力異常；（b）異常的局部波數(shù)；（c）算術(shù)平均歸一化局部波數(shù)；（d）中值歸一化局部波數(shù)；（e）幾何平均歸一化局部波數(shù)；（f）調(diào)和平均歸一化局部波數(shù)Fig.3 （a）Gravity anomaly generated by a horizontal cylinder；（b）Local wavenumber of the data in Fig.3a；（c）Normalized local wavenumber based on arithmetic mean；（d）Normalized local wavenumber based on median value；（e）Normalized local wavenumber based on geometric mean；（f）Normalized local wavenumber based on harmonic mean

從圖4中可以看出，歸一化局部波數(shù)法對(duì)于磁異常也有較好的應(yīng)用效果，其歸一化局部波數(shù)的最大值處埋深為9.8m，與異常體的實(shí)際深度相一致，具有不受磁化方向干擾的優(yōu)勢(shì).

在實(shí)際數(shù)據(jù)解釋中噪聲和區(qū)域場(chǎng)的干擾是不可避免的.為了驗(yàn)證方法的穩(wěn)定性，在圖4模型中部加入一埋深為20m，半徑為3m的圓柱體，磁化強(qiáng)度為20A／m，并加入均值為0，方差為1nT的隨機(jī)噪聲，磁異常如圖5a所示.采用圖1所示方法計(jì)算得到的局部波數(shù)如圖5b所示，可以看出局部波數(shù)能很好地消除區(qū)域異常的影響.利用歸一化局部波數(shù)法進(jìn)行含噪磁異常的解釋?zhuān)▓D5c-5f）.從反演結(jié)果中可以看出，本文方法受噪聲和區(qū)域場(chǎng)的干擾小，能有效地完成異常的解釋工作.

下面給出歸一化局部波數(shù)法在較為復(fù)雜情況下的應(yīng)用效果.地下水平位置40m、80m、120m、160m和220m存在埋深分別為7m、9m、12m、15m和15m的垂直薄板狀體，其寬度為1m，磁化強(qiáng)度均為20A／m，磁傾角為70°時(shí)產(chǎn)生的磁異常如圖6a所示.

反演結(jié)果得到異常體的水平位置分別為40m、80m、120m、160m和220m，根據(jù)結(jié)果最大值判斷出地質(zhì)體的埋深分別為6.9m、8.8m、11.5m、13.3m和14.4m.可以看出當(dāng)異常體之間距離與異常體埋深的比值較小時(shí)，受其他鄰近異常體的干擾，反演結(jié)果的準(zhǔn)確性較低.當(dāng)異常體之間的水平距離較大時(shí)，埋深較大的異常體也可以得到準(zhǔn)確的結(jié)果，所以反演結(jié)果的準(zhǔn)確性依賴(lài)于異常體埋深與異常體之間水平距離的比值.

圖4 （a）三個(gè)板狀體產(chǎn)生的磁異常；（b）磁異常局部波數(shù)；（c）算術(shù)平均歸一化局部波數(shù)；（d）中值歸一化局部波數(shù)；（e）幾何平均歸一化局部波數(shù)；（f）調(diào)和平均歸一化局部波數(shù)Fig.4 （a）Original magnetic anomaly generated by three dikes；（b）Local wavenumber of the data；（c）Normalized local wavenumber based on arithmetic mean；（d）Normalized local wavenumber based on median value；（e）Normalized local wavenumber based on geometric mean；（f）Normalized local wavenumber based on harmonic mean

4 實(shí)際數(shù)據(jù)應(yīng)用

將本文提出的歸一化局部波數(shù)法應(yīng)用于實(shí)際數(shù)據(jù)的解釋.圖7a為埃及紅海西部邊緣Hamrawien地區(qū)一條實(shí)測(cè)磁異常剖面（Salem etal.，1999），點(diǎn)距為10m.現(xiàn)已驗(yàn)證地下存在兩個(gè)板狀體，其埋深分別為555.7m和441.2m.圖7b為磁異常的局部波數(shù).圖7c—7f為不同方法獲得的歸一化局部波數(shù)，反演結(jié)果顯示地質(zhì)體中心的水平位置分別為4530m和14860m，其埋深分別為546m和447m，歸一化局部波數(shù)法的反演結(jié)果與實(shí)際埋深相接近，說(shuō)明該方法具有良好的實(shí)際應(yīng)用效果.

從反演結(jié)果中可以看出較深地質(zhì)體的歸一化局部波數(shù)幅值相對(duì)較小，反映不是很清晰.為了清晰地獲得不同埋藏深度地質(zhì)體的位置，對(duì)歸一化函數(shù)進(jìn)行簡(jiǎn)單的變形，采用分段函數(shù)來(lái)表示：

不同的地質(zhì)體采用異常體范圍內(nèi)的歸一化值，可有效地凸顯出不同深度異常體的效應(yīng)，能更清晰地顯示異常體的分布，因此在應(yīng)用本文方法進(jìn)行復(fù)雜異常解釋時(shí)具有一定的局限性.利用改進(jìn)的歸一化局部波數(shù)法對(duì)異常進(jìn)行反演（圖8）.

從圖8中可以看出，通過(guò)分段歸一化函數(shù)可以使改進(jìn)的局部波數(shù)法能更清晰地反映地下地質(zhì)體的分布，能更好地完成異常的解釋工作.

圖5 （a）含噪磁異常；（b）磁異常局部波數(shù)；（c）算術(shù)平均歸一化局部波數(shù)；（d）中值歸一化局部波數(shù)；（e）幾何平均歸一化局部波數(shù)；（f）調(diào)和平均歸一化局部波數(shù)Fig.5 （a）Noise-bearing magnetic anomalies；（b）Local wavenumber of the data；（c）Normalized local wavenumber based on arithmetic mean；（d）Normalized local wavenumber based on median value；（e）Normalized local wavenumber based on geometric mean；（f）Normalized local wavenumber based on harmonic mean

利用歸一化局部波數(shù)法對(duì)朱日和地區(qū)磁異常進(jìn)行解釋?zhuān)湓即女惓Ｈ鐖D9所示.

計(jì)算不同深度（0m，10m，30m，50m，100m，150m）磁異常的幾何平均歸一化局部波數(shù)（圖10）.

根據(jù)圖10所示不同深度歸一化局部波數(shù)可以看出，大部分礦體埋藏深度在50～150m范圍內(nèi)，且異常有向東延伸的趨勢(shì).

5 結(jié)論

本文提出歸一化局部波數(shù)法進(jìn)行重磁異常的解釋工作，并提出了多種不同的歸一化方式.通過(guò)理論模型試驗(yàn)證明本文方法能有效地完成異常的反演工作，且通過(guò)對(duì)比發(fā)現(xiàn)算術(shù)平均歸一化局部波數(shù)所獲得結(jié)果的分辨率較低.將該方法應(yīng)用于實(shí)際數(shù)據(jù)的解釋?zhuān)@得了地下地質(zhì)體的深度，與已知深度相一致，并針對(duì)不能清晰地反映較深地質(zhì)體的缺點(diǎn)采用分段歸一化的方式來(lái)使異常體得到更加清晰地顯示.

圖6 （a）原始磁力異常；（b）磁異常局部波數(shù)；（c）算術(shù)平均歸一化局部波數(shù)；（d）中值歸一化局部波數(shù)；（e）幾何平均歸一化局部波數(shù)；（f）調(diào)和平均歸一化局部波數(shù)Fig.6 （a）Synthetic magnetic anomalies generated by multiple sources；（b）Local wavenumber of the data；（c）Normalized local wavenumber based on arithmetic mean；（d）Normalized local wavenumber based on median value；（e）Normalized local wavenumber based on geometric mean；（f）Normalized local wavenumber based on harmonic mean

圖7 （a）加拿大安大略省北部實(shí)測(cè)磁異常；（b）磁異常局部波數(shù)；（c）算術(shù)平均歸一化局部波數(shù)；（d）中值歸一化局部波數(shù)；（e）幾何平均歸一化局部波數(shù)；（f）調(diào)和平均歸一化局部波數(shù)Fig.7 （a）Measured magnetic anomaly；（b）Local wavenumber of the data；（c）Normalized local wavenumber based on arithmetic mean；（d）Normalized local wavenumber based on median value；（e）Normalized local wavenumber based on geometric mean；（f）Normalized local wavenumber based on harmonic mean

圖8 （a）算術(shù)平均歸一化局部波數(shù)；（b）中值歸一化局部波數(shù)；（c）幾何平均歸一化局部波數(shù)；（d）調(diào)和平均歸一化局部波數(shù)Fig.8 （a）Normalized local wavenumber based on arithmetic mean；（b）Normalized local wavenumber based on median value；（c）Normalized local wavenumber based on geometric mean；（d）Normalized local wavenumber based on harmonic mean

圖9 實(shí)測(cè)磁異常Fig.9 Measured magnetic anomaly of Zhurihe area

圖10 不同深度歸一化局部波數(shù)Fig.10 Normalized local wavenumber of different depths

Fedi M，F(xiàn)lorio G.2002.A stable downward continuation by using the ISVD method.Geophysical Journal International，151（1）：146-156，doi：10.1046／j.1365-246X.2002.01767.x.

Keating P.2009.Improved use of the local wavenumber in potential-field interpretation.Geophysics，74（6）：L75-L85，doi：10.1190／1.3242270.

Ma G Q.2013.Improved local wavenumber methods in the interpretation of potential field data.Pure and Applied Geophysics，170（4）：633-643，doi：10.1007／s00024-012-0551-z.Meng P，Qin T，Wu Y H.2003.Study of nonuniqueness of the solution to normalized total gradient method.Geophysical Prospecting for Petroleum（in Chinese），42（2）：252-253.

Rapolla A，Cella F，F(xiàn)edi M，etal.2002.Improved techniques in data analysis and interpretation of potential fields：examples of application in volcanic and seismically active areas.Annals of Geophysics，45（6）：733-751，doi：10.4401／ag-3541.

Salem A，Elsirafi A，Ushijima K.1999.Design and application of high-resolution aeromagnetic survey over Gebel Duwi area and its offshore extension.Egypt：Memoirs of the Faculty of Engineering.Kyushu University，59（3）：201-213.

Salem A，Smith R S.2005.Depth and structural index from normalized local wavenumber of 2Dmagnetic anomalies.Geophysical Prospecting，53（1）：83-89，doi：10.1111／j.1365-2478.2005.00435.x.

Salem A，Ravat D，Smith R，etal.2005.Interpretation of magnetic data using an enhanced local wavenumber（ELW）method.Geophysics，70（2）：L7-L12，doi：10.1190／1.1884828.

Salem A，Williams S，F(xiàn)airhead D，etal.2008.Interpretation of magnetic data using tilt-angle derivatives.Geophysics，73（1）：L1-L10，doi：10.1190／1.2799992.

Shen Q X，Wang Z Q，Li R.2010.Simple analysis of the normalized full gradient of gravity anomaly in the exploration of shallow reservoir.Contributions to Geology and Mineral Resources Research（in Chinese），25（1）：72-75.

Smith R S，Thurston J B，Dai T，etal.1998.iSPI—The improved source parameter imaging method.Geophysical Prospecting，46（2）：141-151，doi：10.1046／j.1365-2478.1998.00084.x.

Thurston J B，Smith R S.1997.Automatic conversion of magnetic data to depth，dip，and susceptibility contrast using the SPI method.Geophysics，62（3）：807-813，doi：10.1190／1.1444190.

Xiao P F，Li M，Xu S Z，etal.2006.Stable solution of normalized total gravity gradient.Oil Geophysical prospecting（in Chinese），41（5）：596-598.

Xiao Y M，Zhang L X.1984.The application of normalized gravity gradient in hydrocarbon exploration.Oil Geophysical Prospecting（in Chinese），（3）：247-254.

Zeng H L，Li X M，Yao C L，etal.1999.The modified normalized full gradient of gravity anomalies and its application to Shengli oil field，East China.Petroleum Exploration and Development（in Chinese），26（6）：1-6.

附中文參考文獻(xiàn)

孟平，秦瞳，吳云海.2003.關(guān)于歸一化總梯度異常多解性問(wèn)題的研究.石油物探，42（2）：252-253.

沈慶夏，王志強(qiáng)，李瑞.2010.重力歸一化總梯度法勘探淺層油氣藏淺析.地質(zhì)找礦論叢，25（1）：72-75.

肖鵬飛，李明，徐世浙等.2006.重力歸一化的穩(wěn)定解法.石油地球物理勘探，41（5）：596-598.

肖一鳴，張林詳.1984.重力歸一化總梯度法在尋找油氣中的應(yīng)用.石油地球物理勘探，（3）：247-254.

曾華霖，李小孟，姚長(zhǎng)利等.1999.改進(jìn)的重力歸一化總梯度法及其在勝利油區(qū)油氣藏探測(cè)中的應(yīng)用效果.石油勘探與開(kāi)發(fā)，26（6）：1-6.