低緯度磁異常波數域化極方法機理分析與功能定位

柴玉璞 萬海珍

（中國石油東方地球物理公司綜合物化探處，河北涿州072751）

0 引言

非垂直向下的磁化強度（矢量）和非垂直向下的磁異常場強度（矢量），致使磁異常形態復雜化，從而使磁異常與其場源的關系變得很不直觀。隨著磁傾角變小，此問題越來越嚴重，給磁異常解釋帶來許多困難。Baranov［1］提出的磁異常化極概念（非垂直向下的磁化強度和磁異常場強度均轉換為垂直向下，就像將磁源體搬到北極所產生的磁異常場一樣，故稱化極），奠定了解決這一問題的理論基礎；Bhattacharyya［2］提出的位場波數域轉換理論和方法，使磁異常化極得以實現。然而，隨著磁傾角的減小，化極過程不穩定現象越發嚴重，使低緯度磁異常化極效果很不理想。于是各種低緯度波數域化極方法相繼出現。這些方法可分為三類：①波數域壓制型化極法，如偽傾角化極法［3］、平方余弦壓制化極法［4］、維納濾波化極法［5］及雙曲正弦壓制化極法［6］等；②波數域反演化極法［7］；③波數域偏移抽樣化極法［8-10］。

本文主要討論三個問題：①從理論上證明壓制型化極法肯定不能用于低緯度磁異常定量解釋，但可用于定性解釋；②從理論上揭示，采用等效源法實現的波數域反演化極和采用理論化極算子實現的波數域反演化極，前者效果總是優于后者的根本原因，給出該方法難以用于低緯度磁測資料定量解釋的功能定位；③闡明波數域偏移抽樣化極法的基本思想和方法原理，指出該方法相當苛刻的應用條件。

1 可用于定性解釋的低緯度化極方法——壓制型化極法

從位場波數域轉換算法誤差方程［8-10］出發，導出了實用含噪磁異常波數域化極算法誤差方程［9］（附錄A）

式中：Re（·）表示取實部；DFT（·）、IDFT（·）分別表示正、反離散傅里葉變換；R（·）、R′（·）分別為理論化極算子和實用化極算子；g′1（nΔ）、g2（nΔ）分別為含噪待化極磁異常場和理論化極磁異常場；n（n1，n2）、m（m1，m2）分別為空間域和波數域坐標矢量，n1＝1，2，…，N1表示空間域x方向的離散坐標點，n2＝1，2，…，N2表示空間域y方向的離散坐標點，m1＝1，2，…，N1表示波數域u方向的離散坐標點，m2＝1，2，…，N2表示波數域v方向的離散坐標點；是空間域采樣間隔對角矩陣，Δ1、Δ2分別為空間域x、y方向的采樣間隔；是波數域采樣間隔對角矩陣，d1、d2分別為波數域u、v方向的采樣間隔；表示波數域和空間域數據規模對角矩陣；e（nΔ）為噪聲，也稱觀測誤差。

利用式（1），可分析各種壓制型化極法的基本功能。式中左端是以R′（md）為化極算子的實用化極算法表達式；右端第一項是理論化極異常場，第二項是化極損失項，第三項是觀測誤差。在低緯度區，理論化極算子R（md）實際是一放大作用極強的扇通濾波器，其中心線方位角θ0＝90°＋D（D為磁偏角）；實用化極算子R′（md）是R（md）經壓制的結果，也是一放大作用很強的扇通濾波器，中心線與R（md）的中心線重合。數學推導和計算可證明

也是一種濾波器。不過，式（2）對偽傾角化極法［3］和平方余弦壓制化極法［4］而言是扇通濾波器（圖1），對維納濾波化極法［5］和雙曲正弦壓制化極法［6］而言是準扇通濾波器。這四個濾波器的中心線都與R（md）、R′（md）的中心線重合，但濾波強度弱很多。

筆者在研究磁異常化極中的誤差規律時發現，像理論化極算子R（md）這類扇通濾波器，可分解為一個方向低通濾波器（中心線方位角為θ0＝90°＋D）和兩個高通濾波器（其中濾波作用較強的濾波器中心方位角為θ0＝90°＋D，濾波作用較弱的濾波器中心方位角為θ0＝D）［10］。其實實用化極算子R′（md）也可做同樣的分解。扇通濾波器可做同樣的分解，只是中心線方位角為θ0＝D的濾波器消失，因為該扇通濾波器的極小值為零，而R（md）、R′（md）的極小值都不為零。

圖1 扇通濾波器濾波等值線圖

下面對分解出的各個濾波器的功能及其對化極結果的影響作分析。

（3）濾波器R′（md）中較強的方向低通濾波器（濾波器3），會導致化極結果圖中出現沿磁偏角方向、正負相間的平行線狀等值線。這是濾波器3對隨機誤差所形成的、遍布全圖的正負等值線小圈沿磁偏角方向極度拉長的結果。因為隨機誤差所形成的小圈，高頻成分較弱，所以R′（md）中的兩個高通濾波器表現不明顯。

上述三個濾波器的濾波強度都與壓制因子的強度相關。壓制因子增強，則濾波器3的作用減弱，濾波器1和濾波器2的作用增強。也就是說，壓制因子增強，會削弱隨機誤差對化極結果的干擾，但同時會使化極損失效應增強。據此可得如下結論。

（1）當增強壓制強度，使濾波器3所產生的平行線異常減弱到可忽略不計的情況下，濾波器2將大幅度降低化極結果中正異常幅度，而小幅度地增強磁偏角方向上的邊部負異常。

（2）此種情況下，濾波器1對化極結果在垂直于磁偏角方向上有一定程度的圓滑作用。

（3）在上述條件下，濾波器1、濾波器2將一定程度地影響化極異常的形態，但不影響化極異常的中心，也基本不影響化極結果中正異常的分布范圍。

為了驗證上述理論的正確性，設計一方柱體組合模型進行試驗，模型參數見文獻［11］。模型的磁傾角和磁偏角均為0°，觀測范圍為50km×50km，點、線距均為1km。圖2是該模型的總磁異常（ΔT），圖3是磁傾角為90°的理論化極異常。圖4是采用偽傾角法［3］計算的化極結果，其余三種壓制型化極方法的化極結果與偽傾角法化極結果大同小異，故本文未展示。粗看上去，圖3與圖4差異非常大：圖4中化極結果的正異常幅值比圖3小得多，而磁偏角方向上的邊部負異常則比圖3強；此外，垂直于磁偏角方向上的異常等值線出現一定程度的圓滑，而且變得稀疏。然而仔細觀察可以發現，圖4所示的化極結果中正異常的中心位置和分布范圍與圖3基本一致。

圖2 方柱體組合模型磁源的總磁異常（ΔT）分布圖

理論分析和模型計算結果表明：各種波數域壓制型化極法不能用于低緯度磁異常數據的定量解釋，因為這些方法的化極結果與理論化極結果，無論異常形態還是量值，差異都相當大；但是這些方法完全可用于定性解釋，因為化極結果中正異常的中心位置和分布范圍與理論化極異常的基本一致。

圖5是秘魯A地區實測磁異常（ΔT）。該地區的地磁傾角和磁偏角分別為－2°和0°，觀測點、線距均為2km，資料覆蓋范圍為100km×288km。圖6是圖5采用偽傾角法的化極結果，該圖為秘魯A地區的地質勘探提供了重要信息。

圖3 圖2的理論化極異常（磁傾角為90°）

圖4 采用偽傾角法對圖2磁異常的化極結果

圖5 秘魯A地區實測磁異常圖

圖6 利用偽傾角法對圖5的化極異常圖

2 有可能用于定量解釋的低緯度化極方法：波數域反演化極法和波數域偏移抽樣化極法

2.1 波數域反演化極法

從原理上看，波數域反演化極法［7］可用于低緯度磁異常數據的定量解釋。該方法的核心是求解依據基本關系

建立的方程組

式中：G是以1/R（md）為元素的對角矩陣；m′和d′分別為模型向量和觀測值向量；W為加權函數矩陣；S為加權差分矩陣；上標“H”和“T”分別表示復共軛轉置和轉置；μ表示正則化參數。由理論化極算法誤差方程［10］可知，式（3）嚴格成立的條件是G以T（md）為元素，而不是1/R（md），因為式（3）中的d′和m′分別是待化極異常的有限離散譜和化極目標異常的有限離散譜。T（md）與1/R（md）的關系（推導過程詳見附錄A）為

Li等［7］用1/R（md）作為G的元素所獲得的化極結果不理想，改用等效源法計算G的元素，結果大為改善。式（5）從理論上解釋了這一現象產生的原因：采用等效源法得到的G中，至少包含了計算窗外數據的部分影響，而采用1/R（md）計算的G，則完全忽略了計算窗外數據的影響。文獻［7］中的模型算例是對一直立方柱體的斜磁化磁異常場進行化極，所以采用一小直立方柱體為等效源計算G的元素，化極效果相當好。因為小方柱體磁異常的窗外數據與大方柱體磁異常的窗數據有很強的相似性。

由理論化極算法誤差方程導出的式（5），為波數域反演化極過程中采用等效源法計算G的元素，提供了理論依據。然而對實際資料而言，窗外數據難以較精確地模擬。這可能就是該方法難以應用于實際資料處理的根本原因。

2.2 波數域偏移抽樣化極法

由含噪實用化極算法方程的推導過程（附錄A）可知：式（1）中，若將DFT替換為移樣離散傅里葉正變換算子DFT0η，將IDFT替換為移樣離散傅里葉反變換算子IDFT0η（下標“0η”表示僅波數域偏移抽樣，偏移參數η為二維行向量），將m替換為m＋η，并令R′［（m＋η）d］＝R［（m＋η）d］，即得波數域偏移抽樣化極算法方程

式中，左端為化極算法表達式，右端第一項為理論化極結果項，第二項為隨機誤差。波數域偏移抽樣化極方法的基本思想是：測線沿磁偏角方向布設，使磁偏角為0°，然后通過沿u軸方向偏移v軸半個采樣間隔（0.5d1）抽樣，大大減小化極算子的量值，以控制隨機誤差對化極結果的影響。

由上述基本思想和式（6）的化極算法表達式可知，波數域偏移抽樣化極法的具體做法包括如下步驟：

（1）對磁偏角為0°的磁測資料做DFT0η變換，偏移抽樣參數η＝（0.5，0）；

（2）將步驟（1）的結果乘以理論化極算子R［（m＋η）d］；

（3）對步驟（2）的結果做IDFT0η變換，偏移參數η＝（0.5，0）；

（4）對步驟（3）的結果取實部，即得偏移抽樣化極結果。

波數域偏移抽樣化極算法方程右端只有兩項，一項是理論化極結果項，一項是觀測誤差項。因此把觀測誤差項的量值減至可忽略的程度，是該方法成功的關鍵。而誤差項量值的減小，是通過沿u軸方向偏移v軸半個采樣間隔（0.5d1）抽樣實現的（d1越大，減小量越大）。這種抽樣方式可使化極算子的量值大大減小，而且隨著N1（d1＝1/N1Δ1）的減小，量值還會進一步減小。然而當沿磁偏角方向上的網格數N1降至64時，化極算子的量值仍然相當大，以至于要求資料相對精度高于0.2%，才能保證觀測誤差項量值小到可忽略不計。

由上述分析可以看出，波數域偏移抽樣化極法可用于低緯度磁異常的定量解釋，但條件相當苛刻：磁偏角方向上網格數受限，且資料精度要求很高。因此在目前資料精度條件下，波數域偏移抽樣化極法不可能用于低緯度區大面積磁測資料的定量解釋。隨著資料精度的不斷提高，將來該方法有可能用于低緯度區小面積高精度磁測資料的定量解釋。

3 結論

（1）理論研究表明，各種波數域壓制型化極法不能應用于低緯度磁測資料的定量解釋，但可以成為低緯度磁測資料定性解釋的有效工具。這一理論成果的獲得，實用化極算法誤差方程的導出是基礎，扇通濾波器定性分解（可分解為一個方向低通濾波器和兩個或一個高通濾波器）觀點的提出和運用是關鍵。這一源于化極問題的原創性概念和觀點，為其他涉及扇通濾波器問題的解決提供了新的思路和方法。

（2）文中的式（5）為波數域反演化極中采用等效源法計算G的元素，提供了嚴格的理論依據。理論上，該方法可用于低緯度磁測資料的定量解釋，然而實際資料的窗外數據難以較精確地模擬。因此，該方法難以應用于實際數據的化極。該方法需求解超大型方程組，可能是難以實用化的另一原因。

（3）波數域偏移抽樣化極法，理論上可用于低緯度磁測資料的定量解釋，但由于應用條件過于苛刻，在目前數據精度條件下不可能用于低緯度、大面積磁測資料的定量解釋。但是隨著未來資料精度的提高，該方法有可能用于低緯度區、小面積、高精度磁測資料的定量解釋。

感謝中國石油集團東方地球物理公司賈繼軍高級工程師在研究初期在編程方面給予的幫助，楊戰軍高級工程師在圖件繪制方面的幫助，以及孫喜明高級工程師在資料準備方面提供的幫助！

附錄A 實用化極算法方程及T（md）表達式的推導

假設g1（x，y）和G1（u，v）分別表示無噪待化極異常及其譜，g2（x，y）和G2（u，v）分別表示化極目標異常及其譜，那么g2（x，y）和G2（u，v）的二維IDFT算法誤差方程［9］為

g1（x，y）和G1（u，v）的二維DFT算法誤差方程［9］可寫作

式（A-2）兩邊乘以R（md），并利用關系

將式（A-4）代入式（A-1），并移項，得理論化極算法誤差方程

在采用實用化極算子的情況下，式（A-6）右端前四項量值很小，可忽略不計。所以實用無噪磁異常化極算法誤差方程簡化為

并將g2（nΔ）由方程左端移到右端，即得正文中的式（1），即含噪實用化極算法方程。

對式（A-5）做如下變換：①將左端的g2（nΔ）移到右端；②對方程兩端作DFT；③方程兩端同除以R（md）DFT［g2（nΔ）］，可得正文中的式（5），即T（md）的表達式。