王彥國 吳姿穎 鄧居智 楊亞新 黃 松 羅 瀟

(東華理工大學放射性地質與勘探技術國防重點學科實驗室，江西南昌 330013)

1 引言

位場數(shù)據(jù)處理是重磁資料進行地質—地球物理綜合解釋的基礎，而反演是數(shù)據(jù)處理的重要環(huán)節(jié)。常用的位場反演方法包括基于網(wǎng)格剖分的物性反演法和快速自動反演方法。

基于網(wǎng)格剖分的物性反演方法是基于反演理論、在最小二乘意義下使目標函數(shù)達到極小的(非)線性反演。這類方法需要引入先驗約束條件，需大量的計算機內(nèi)存和較長的計算時間，制約著這類反演方法的應用，尤其在地質勘探程度較低的地區(qū)或無地質約束區(qū)域，難以獲得良好的反演結果[1]。近年來，Li等[2]、姚長利等[3,4]提出的改進措施在一定程度上提高了這類方法的計算效率，改善了反演效果。

常用的快速自動反演方法主要有歐拉反褶積法和歸一化總梯度法，該類方法可以在無地質約束條件下快速計算場源參數(shù)信息。其中，歐拉反褶積是以位場和位的導數(shù)滿足歐拉齊次方程為理論基礎，由Peters[5]首先提出，Thompson[6]進一步完善，Reid等[7]推廣至三維資料處理中的一種反演方法。該方法能自動或半自動地反演出場源質心、上頂或邊界位置，還可以利用反演的構造指數(shù)推斷場源的幾何形狀。由于該方法具有較強的靈活性和實用性，目前已成為了地球物理工作者常用的反演方法和重磁資料解釋工具。但事實上，歐拉反褶積是一個靈敏度較高的方程式，導數(shù)計算結果或計算點位稍有偏差，便會使反演結果產(chǎn)生較大的偏離，從而影響反演精度[8]。近年來，國內(nèi)外學者對歐拉方程進行了各種形式的改進：與解析信號相結合的An-EUL反演法[9]，與Tilt angle相結合的Tilt-Euler法[10],以及其他的改進措施[11-16]，這些改進方法在一定程度上提高了歐拉反褶積法的可靠性。歸一化總梯度法是建立在導數(shù)換算和向下延拓基礎之上的歸一化形式，是由前蘇聯(lián)學者別列茲金提出，該方法對油氣儲集層有著較好的識別效果[17]； 肖一鳴等[18，19]將該方法引進到國內(nèi)，并介紹了初步應用情況； 侯重初等[20]用傅氏變換代替傅氏級數(shù)計算歸一化總梯度，并進行了詳細的模型試算；張鳳旭等[21]利用Hilbert變換計算了重力歸一化總梯度，在理論上提高了異常分辨率；肖鵬飛等[22]利用積分迭代法進行歸一化總梯度中向下延拓的穩(wěn)定計算；張鳳琴等[23]采用了DCT變換計算重力歸一化總梯度；郭燦燦等[24]利用泰勒級數(shù)迭代法代替積分迭代法，進行歸一化總梯度的向下延拓計算；蘇超等[25]采用幾何平均代替算術平均進行空間總梯度的歸一化處理。這些改進措施在一定程度上提高了歸一化總梯度法的計算穩(wěn)定性和反演精度，然而，該方法僅能識別淺部地質體的位置信息，無法反映出深部地質體的深度信息，或計算的深度值與真實值偏差過大，也無法判斷場源的幾何形狀。

針對常規(guī)歸一化總梯度法的缺點，本文提出了基于冪次平均的離散歸一化總梯度法，并通過模型試驗和實際資料處理驗證了方法的可行性、有效性和實用性。

2 基本原理

2.1 傳統(tǒng)歸一化總梯度法

傳統(tǒng)歸一化總梯度法是別列茲金提出的，其表達式為

(1)

(2)

(3)

2.2 迭代濾波

由于導數(shù)換算和向下延拓都會不同程度地放大高頻干擾成分，需要通過低通濾波增強計算結果的穩(wěn)定性。傳統(tǒng)的低通濾波算子是別列茲金提出的圓滑濾波因子

(4)

式中：N為諧波數(shù)；k=1,2,…,N。該濾波因子在壓制高頻成分時，低頻成分在一定程度上也有所削弱。為此，本文選擇如下迭代濾波因子提高延拓的計算精度和穩(wěn)定性

(5)

式中：α為正則化因子，α越大，對高頻成分壓制越強，一般情況下α≥1(本文取α=1)；n為迭代次數(shù)。

將式(5)代入式(3)，有

(6)

2.3 冪次平均歸一化總梯度

傳統(tǒng)歸一化函數(shù)MAs(z)為算術平均，本文采用冪次平均代替算術平均進行歸一化處理，即

(7)

式中p為冪次數(shù)，當p=1時，EAs(z)便是算術平均MAs(z)。不同形狀的地質體具有不同的最佳p值，p=1對應水平圓柱體，p=2對應巖脈，p=4則對應臺階(下文模型試驗將予以證實)?；趦绱纹骄臍w一化總梯度表達式為

(8)

2.4 歸一化總梯度的離散化

由于式(1)和(8)只能識別最淺部的地質體信息，較深部的地質體則無法被同時檢測出來，或計算結果與實際情況偏差過大。為了提高方法的實用性，下面對式(8)進行離散化處理。

設觀測面上存在Q個明顯的解析信號As的極大值，將測線離散成Q段，每段包含一個解析信號極大值，兩個極大值之間的極小值點定義為離散點(圖1)。

那么，離散化后第q段的歸一化總梯度為

(9)

圖1 離散化處理示意圖

2.5 方法實現(xiàn)步驟

基于冪次平均的離散歸一化總梯度計算流程如下：

(1)給定初值迭代次數(shù)n(可取1)，采用迭代濾波法計算地下空間的導數(shù)Uxz(x,z)、Uzz(x,z)及As(x,z)；

(2)根據(jù)As(x,0)的異常形態(tài)，利用極大值個數(shù)將測線離散成Q段；

(6)改變區(qū)間段位置，重復上述步驟(3)～(5)，完成所有區(qū)間段的計算；

3 模型試驗

3.1 冪次數(shù)與地質體幾何形狀關系

為了驗證冪次數(shù)與地質體形狀存在相關性，分別設計了臺階、巖脈和無限長水平圓柱體等三種常用二度體模型進行試驗，并通過改變地質體參數(shù)進一步說明這種相關性。

3.1.1 水平圓柱體模型

建立兩個單一水平圓柱體模型，模型1、模型2的質心埋深分別為1km、2km，磁化傾角分別為45°和90°，半徑均為0.5km。圖2給出了這兩個水平圓柱體模型經(jīng)磁異常轉換得到的max(UH)及其對應的深度值與迭代次數(shù)n和冪次數(shù)p的關系曲線。從圖2a、圖2c可以看出，無論p取何值，max(UH)均隨迭代次數(shù)n的增加而先增大后減小，存在一個極值點； 當冪次數(shù)p增大時， max(UH)逐漸減小。從圖2b、圖2d可以看出，max(UH)對應的深度值隨迭代次數(shù)增加先緩慢遞增，再迅速增大，存在一個明顯的轉折點，而該轉折點對應的深度值與max(UH)曲線上極大值對應的深度值一致。由于p=1與p=2時max(UH)曲線上極大值對應的深度值相同，因此可認為圓柱體的最佳冪次數(shù)p為1，此時最佳迭代次數(shù)的max(UH)對應的深度值分別為0.9km和1.8km，均與模型深度1.0km和2.0km存在10%的相對誤差。圖3a、圖3b分別是兩個水平圓柱體模型磁異常歸一化總梯度場，可以看出，UH場關于水平圓柱體質心埋深位置左右呈對稱分布，極大值圈閉位置與圓柱體質心位置基本吻合。

3.1.2 巖脈模型

構建了兩個巖脈模型(巖脈1、巖脈2)，巖脈的傾斜度分別為45°和90°，上頂埋深分別為1.0km、2.0km，磁化傾角分別為45°和90°。圖4是兩個巖脈模型經(jīng)磁異常轉換得到的max(UH)及其對應的深度值與迭代次數(shù)n和冪次數(shù)p的關系曲線。從圖4a、圖4c可以看出，max(UH)仍隨迭代次數(shù)n的增加而先增大后減小，存在一個極值點；當冪次數(shù)p增大時， max(UH)同樣逐漸減小。從圖4b、圖4d可以看出， max(UH)對應的深度值也隨迭代次數(shù)增加而先緩慢遞增，然后迅速增大，存在一個明顯的轉折點，而該轉折點對應的深度值與max(UH)曲線上極大值對應的深度值基本一致。由于p=2與p=3時max(UH)曲線上極大值對應的深度值相同，因此巖脈的最佳冪次數(shù)p=2，此時最佳迭代次數(shù)下的max(UH)所對應的深度值分別為1.0km和1.9km，與模型深度1.0km和2.0km基本吻合。圖5是兩個巖脈模型的磁異常歸一化總梯度場，可見UH場關于巖脈上頂埋深位置左右呈對稱分布，極大值圈閉位置與模型上頂埋深位置基本重合。

圖2 兩個水平圓柱體模型的max(UH)及其對應的深度隨p和n的變化曲線(a)模型1 max(UH)隨參數(shù)p和n的變化曲線； (b)模型1 max(UH)對應的深度值隨參數(shù)p和n的變化曲線；(c)模型2 max(UH)隨參數(shù)p和n的變化曲線； (d)模型2 max(UH)對應的深度值隨參數(shù)p和n的變化曲線

圖3 水平圓柱體模型1(a)及模型2(b)磁異常冪次平均歸一化總梯度場

圖4 巖脈模型1 max(UH)隨p和n的變化曲線(a)、 巖脈模型1 max(UH)所對應的深度隨p和n的變化曲線(b)、 巖脈模型2 max(UH)隨p和n的變化曲線(c)及巖脈模型2 max(UH)所對應的深度隨p和n的變化曲線(d)

圖5 巖脈模型1(a)、巖脈模型2(b)磁異常冪次平均歸一化總梯度場

3.1.3 臺階模型

臺階模型1、模型2的傾斜度分別為45°和60°，上頂埋深分別為1.0km和2.0km，磁化傾角分別為30°和60°。圖6為兩個臺階模型的磁異常max(UH)及其對應的深度與迭代次數(shù)n和冪次數(shù)p的關系曲線。從圖6a、圖6c可以看出，max(UH)均隨n的增加而先增大后減小，存在極值點； 當冪次數(shù)p增大時，max(UH)則逐漸減小。從圖6b、圖6d可以看出，除p=1外，在其他p值時，max(UH)所對應的深度隨迭代次數(shù)的增加先逐漸增大，當達到最佳迭代次數(shù)后迅速增大，因此max(UH)所對應的深度曲線在最佳迭代次數(shù)時存在一個明顯的拐點； 最佳迭代次數(shù)時max(UH)所對應的深度隨p值的增加逐漸接近模型埋深，由于p=4與p=5時的最佳迭代次數(shù)下的max(UH)所對應的深度相同，因此認為p=4是臺階模型的最佳冪次數(shù)，此時max(UH)所對應的深度分別為1.1km和2.1km，均略大于模型深度1.0km和2.0km。圖7是上述兩個臺階模型磁異常歸一化總梯度場等值線圖，可以看出，UH關于臺階上頂埋深位置左右呈對稱分布， max(UH)所對應的位置與臺階上頂埋深位置基本吻合。

上述三組模型試驗結果表明，基于冪次平均的歸一化總梯度法是可行的，可以有效地識別不同形狀地質體的深度參數(shù)； 不同幾何形狀的地質體對應的最佳冪次數(shù)是不同的，可利用該特點估計地質體類型；歸一化總梯度場的分辨率隨著地質體埋深增加而有所降低。

圖6 臺階模型1 max(UH)隨p和n的變化曲線(a)、 臺階模型1 max(UH)所對應的深度隨p和n的變化曲線(b)、臺階模型2 max(UH)隨p和n的變化曲線(c)及臺階模型2 max(UH)所對應的深度隨p和n的變化曲線(d)

3.2 疊加模型試驗

為了進一步驗證新方法的有效性，設計了一個疊加組合模型，組合模型包含了三個不同類型的地質體。地質體1為臺階模型，上頂坐標為(10km，2km)，傾斜角度為60°，磁化傾角為30°，磁化強度為0.1A/m；地質體2為垂直巖脈模型，上頂坐標為(20km，0.5km)，磁化傾角為45°，磁化強度為1A/m；地質體3為水平圓柱體，質心坐標為(30km，1km)，磁化傾角為60°，磁化強度為1A/m。圖8是該組合模型產(chǎn)生的磁異常及磁異常的解析信號，從解析信號可以看出地下存在三個明顯的磁異常體。圖9是整條剖面及離散化剖面的磁異常max(UH)及其對應深度與迭代次數(shù)n和冪次數(shù)p的關系曲線，從圖9a和圖9b中可以看出，冪次數(shù)p=2與p=3時最佳迭代次數(shù)下的max(UH)對應的深度一致，即最佳冪次數(shù)p為2時，此時max(UH)對應的坐標為(20.0km，0.5km)，與巖脈上頂埋深位置完全一致。

圖7 臺階模型1(a)和模型2(b)磁異常冪次平均歸一化總梯度場

圖8 疊加模型磁異常及其解析信號(a)和模型示意圖(b)

圖10a是整條測線的最佳冪次數(shù)、最佳迭代次數(shù)時的磁異常歸一化總梯度場，可以看出，總梯度場僅能反映巖脈的位置信息，無法有效識別臺階和水平圓柱體。依據(jù)地面磁異常解析信號(圖8a)特征，將測線離散化為3段：0～16.0km、16.0～25.5km和25.5～40.0km。由于未離散化的磁異常歸一化總梯度可以較好地識別出16.0～25.5km之間的磁異常體，因此該段無需再處理。圖9c～圖9f分別是第一段(0～16.0km)和第三段(25.5～40.0km)磁異常max(UH)及其深度值隨迭代次數(shù)n和冪次數(shù)p的關系曲線?？梢钥闯觯谝欢蝺绱螖?shù)p=4與p=5時，最佳迭代次數(shù)下的max(UH) 對應的深度值一致，即該段最佳冪次數(shù)為4，與上述臺階模型試驗對應的最佳冪次數(shù)一致， 對應的坐標為(10.0km，2.1km)，與臺階上頂埋深位置基本一致；第三段最佳冪次數(shù)p=1，為水平圓柱體對應的冪次數(shù)，最佳冪次數(shù)時最佳迭代次數(shù)的max(UH)對應的坐標為(30.0km，0.9km)，與理論質心坐標(30.0km，1.0km)較為接近。圖10b是離散歸一化總梯度場，由圖可見，可以根據(jù)極值圈閉識別三個地質體的位置，有效提高了歸一化總梯度法的實用性。

圖9 疊加模型整條剖面max(UH)隨p和n的變化曲線(a)、 疊加模型整條剖面max(UH)對應的深度隨p和n的變化曲線(b)、 疊加模型0～16.0km的max(UH)隨p和n的變化曲線(c)、 疊加模型0～16.0km的max(UH)對應的深度隨p和n的變化曲線(d)、 疊加模型25.5～40.0km的max(UH)隨p和n的變化曲線(e)及疊加模型25.5～40.0km的max(UH)對應的深度隨p和n的變化曲線(f)

圖10 疊加模型整條剖面(a)及離散剖面拼貼(b)的磁異常冪次平均歸一化總梯度場

4 實例應用

為了驗證基于冪次平均的離散歸一化總梯度法的實用性，選取中國內(nèi)蒙古某地區(qū)的地面磁異常場進行試驗，研究區(qū)地表全部被沉積層所覆蓋。圖11是網(wǎng)格化為250m×50m后的磁異常圖，可以看出，研究區(qū)南側存在一條近東西走向的條帶狀磁異常，在西北側還存在一個長軸近東西走向的橢圓狀磁異常。為了查明這兩個磁異常對應的磁源參數(shù)信息，在磁異常圖中選取了兩條南北走向的試驗剖面AA′和BB′，對應的磁異常及解析信號曲線見圖12?？梢钥闯?，AA′剖面上有兩個明顯的解析信號極值，推斷沿這條剖面存在兩個磁性體；BB′剖面上只有一個明顯的解析信號極值，認為地下主要存在一個磁源體。

首先采用基于圓滑濾波因子qm(式(4))的常規(guī)歸一化總梯度法對兩條剖面進行處理。圖13是常規(guī)歸一化總梯度場極值及其深度值隨諧波數(shù)N的變化曲線?？梢钥闯?，在AA′、BB′剖面上，最佳諧波數(shù)N時歸一化總梯度場極值對應的深度分別為100m和400m。圖14是最佳諧波數(shù)時AA′和BB'剖面磁異常的常規(guī)歸一化總梯度斷面圖，可以看出，AA′剖面常規(guī)歸一化總梯度在坐標(4000m，100m)處存在明顯的極大值，分別在坐標(1800m，100m)、(3450m，100m)和(3650m，100m)處也存在極值，但極值圈閉范圍較小，可能是磁性體引起的，但也可能是噪聲壓制不徹底導致的；BB′剖面的常規(guī)歸一化總梯度在坐標(1450m，400m)處存在最大值，還在坐標(2200m，400m)、(4250m，1800m)兩處存在明顯的極值圈閉。

圖11 中國內(nèi)蒙古某地區(qū)磁異常等值線圖

圖12 圖11中AA′剖面(a)和BB′剖面 (b)的磁異常及其解析信號

圖13 常規(guī)歸一化總梯度場極值及其對應的深度隨諧波數(shù)N的變化曲線

(a)AA′剖面常規(guī)歸一化總梯度場極大值與諧波數(shù)N的關系曲線； (b)AA′剖面常規(guī)歸一化總梯度場極大值對應的深度與諧波數(shù)N的關系曲線； (c)BB′剖面常規(guī)歸一化總梯度場極大值與諧波數(shù)N的關系曲線； (d)BB′剖面常規(guī)歸一化總梯度場極大值對應的深度與諧波數(shù)N的關系曲線

再采用基于冪次平均的離散歸一化總梯度法對兩條剖面磁異常進行處理。根據(jù)AA′剖面的解析信號(圖12a)，可把水平位置2850m作為離散點將測線分為兩段。圖15a、圖15b分別是AA′剖面磁異常max(UH)及其對應的深度與冪次數(shù)和迭代次數(shù)的關系曲線；圖15c、圖15d分別是該剖面離散化后0～2850m的磁異常max(UH)及其對應的深度與冪次數(shù)和迭代次數(shù)的關系曲線。可以看出，水平位置4000m處的磁性體最佳冪次數(shù)p=1，推斷地質體類似為水平圓柱體，該最佳冪次數(shù)下最佳迭代次數(shù)時的max(UH)對應的深度為350m；在水平位置為1800m處的磁性體最佳冪次數(shù)p=2，即可認為磁性體近似于巖脈，最佳冪次數(shù)下最佳迭代次數(shù)時的max(UH)對應的深度為550m。從BB′剖面的解析信號圖(圖12b)可以看出，該剖面解析信號異常主要有一個極值點，因此該剖面無需離散化。圖16a、圖16b分別為BB′剖面磁異常max(UH)及其對應的深度隨參數(shù)p和n的變化曲線，從中可知，在水平位置1400m處地下存在一個磁性體，該磁性體的最佳冪次數(shù)p=2，即磁源接近于巖脈形狀，最佳冪次數(shù)下最佳迭代次數(shù)時的max(UH)對應的深度為450m。圖17a、圖17b分別是AA′、BB′剖面的冪次平均離散歸一化總梯度斷面圖，可以看出，AA′剖面的磁異常冪次平均離散歸一化總梯度場在坐標(1800m，550m)和(4000m，350m)處清晰地展示出兩個明顯的極值圈閉，雖然常規(guī)歸一化總梯度(圖14a)也在水平位置1800m和4000m處存在兩個極值，但極值對應的深度(均為100m)與改進方法獲得的深度存在較大偏差。BB′剖面的磁異常冪次平均歸一化總梯度場在坐標(1400m，450m)處存在一個明顯的極值圈閉，在坐標(2200m，400m)也存在一個極值圈閉，可能在該處也存在一個磁性體，只是地表磁異常解析信號分辨率不夠高導致未被識別；常規(guī)歸一化總梯度(圖14b)在水平位置1400m和2200m處識別的磁性體埋深均為400m，與改進方法計算結果基本一致。

圖14 AA′(a)、BB′(b)剖面磁異常常規(guī)歸一化總梯度場

圖15 AA′剖面max(UH)隨p和n的變化曲線(a)、 AA′剖面max(UH)對應的深度隨p和n的變化曲線(b)、 AA′剖面0～2850m磁異常max(UH)隨p和n的變化曲線(c)及AA′剖面0～2850m磁異常max(UH)對應的深度隨p和n的變化曲線(d)

圖16 BB′剖面max(UH) 隨參數(shù)p和n的變化曲線(a)及max(UH)對應的深度隨參數(shù)p和n的變化曲線(b)

圖17 AA′(a)、BB′(b)剖面基于冪次平均的磁異常離散歸一化總梯度

由于沒有其他地質、鉆孔等已知資料驗證上述結果，這里采用了磁法數(shù)據(jù)處理中常用的AN-EUL反褶積法對AA′、BB′兩條剖面進行相應處理，結果分別見圖18和圖19。從圖18可以看出，利用AN-EUL反褶積反演AA′剖面磁數(shù)據(jù)得到了兩個磁性體，坐標分別為(1727.1m，543.8m)、(3839.2m，350.1m)，與冪次平均離散歸一化總梯度法反演的位置基本一致；AN-EUL反褶積在上述兩處反演得到的構造指數(shù)分別為1.1和1.7，即兩處的磁性體分別接近于巖脈(理論構造指數(shù)為1)和水平圓柱體(理論構造指數(shù)為2)，同樣與冪次平均離散歸一化總梯度法估計的地質體類型相一致。從圖19中可以看出，利用AN-EUL反褶積反演BB′剖面，在坐標(1524.6m，480.3m)處反演得到一個類似巖脈的磁性體(反演的構造指數(shù)為1.1)，反演結果與冪次平均歸一化總梯度得到的磁性體參數(shù)同樣吻合較好，這進一步證實了冪次平均離散歸一化總梯度法的有效性和實用性。

圖18 AA′剖面磁異常AN-EUL反演結果(a)位置解; (b)構造指數(shù)解

圖19 BB′剖面磁異常AN-EUL反演結果(a)位置解; (b)構造指數(shù)解

5 結論

本文在分析常規(guī)歸一化總梯度法的優(yōu)缺點基礎上，利用迭代濾波法代替圓滑濾波實現(xiàn)穩(wěn)定向下延拓，采用冪次平均替換算術平均進行歸一化計算，并對歸一化總梯度場進行離散化處理，由此推導出了基于冪次平均的離散歸一化總梯度法。模型試驗證實，基于冪次平均的離散歸一化總梯度法可以用于處理疊加異常場，能夠有效地識別地質體的位置及幾何形狀參數(shù)。實例應用表明，相對于常規(guī)歸一化總梯度法，新方法不僅可以有效地反映疊加異常的場源分布，獲得更高的反演精度，還可以推斷地質體的幾何形狀特征，有效提高歸一化總梯度法的地質解釋能力。