邊緣特征和深度加權約束的重力三維相關成像反演

安國強,魯寶亮,2,3,高新宇,朱武,3,4,李柏森

(1.長安大學 地質工程與測繪學院,陜西 西安 710054; 2海洋油氣勘探國家工程研究中心,北京100028; 3長安大學 西部礦產資源與地質工程教育部重點實驗室,陜西 西安 710054;4自然資源部 生態地質與災害防控重點實驗室,陜西 西安 710054)

0 引言

相關成像是一種利用相關系數定性解釋地質體空間位置和形態的快速成像方法, 相比于常規重力異常反演方法, 其不需要求解大型方程組, 便能夠快速提供異常體的平面位置、深度信息等。相關成像法最初是由Patella[1]提出, 用于自然電位測量數據進行地質目標體的三維成像與解釋; Mauriello等[2-3]分別將該方法推廣到了大地電磁領域和重磁領域; 此后相關成像方法被廣泛應用于地球物理數據處理中[4-6]。從相關成像的原理出發, 郭良輝等[7-9]提出了重力、重力梯度數據三維相關成像方法以及磁異常的三維相關成像方法; 侯振隆等[10]提出基于泰勒級數的相關成像方法, 提高了成像速度。為了得到真實的物性參數, 孟小紅等[11]和閆浩飛等[12]提出了迭代計算, 反演出了地質體的重磁物性參數。針對直接相關成像結果中出現深部發散的問題, 眾多學者提出了不同的改進方法, 主要包括基于異常分離的相關成像[13], 重磁異常垂直梯度三維相關成像[14-15], 重力數據與重力梯度數據的聯合相關成像[13,16-17], 基于窗口數據的相關成像[18], 基于深度加權的相關成像[16,19]。然而, 目前重力異常相關成像仍然存在橫向分辨率低、深度加權函數超參數過多的問題。

在前人的研究基礎上, 本文提出了基于邊緣特征和深度加權約束的重力三維相關成像反演方法, 該方法較好地克服了相關成像橫向分辨率不高以及深部發散的問題。根據重力相關成像的基本原理, 利用垂向導數(VDR, vertical derivative)和解析信號振幅(ASM, analytical signal amplitude)的邊界均衡方法, 得到了異常體深部的水平邊緣特征信息, 并提出了利用重力異常水平邊緣特征: 均衡垂向導數和均衡解析信號振幅來改善重力異常相關成像橫向分辨率的方法; 針對深度加權函數超參數過多的問題,提出了利用一種更為簡潔的深度加權函數來改善重力異常相關成像縱向分辨率的方法。最后通過模型試驗以及對澳大利亞Olympic Dam多金屬礦區實際資料的計算, 證明了約束后的相關成像在橫向分辨率和縱向分辨率方面得到了顯著的提高, 驗證了該方法的正確性和有效性。

1 方法原理

1.1 重力數據相關成像

重力相關成像是通過計算地下某質量點在地面產生的重力異常與實測重力異常之間的歸一化互相關系數, 用以表示該質量點是異常源的概率大小。將地下待成像的空間剖分成均勻網格, 然后從地下淺層到深層計算所有網格節點的相關系數從而可以得到整個測區深度范圍內各個位置的概率分布。其結果能顯示出地下異常地質體的位置和視密度的相對大小[7,13]。

首先將地下成像空間剖分為均勻網格, 并將單個網格視為點質量, 則根據式(1)可以計算地下某一網格q(xq,yq,zq)在地面上任意點(x,y,z)處產生的重力異常:

(1)

其中:G為萬有引力常數;Δσ為剩余密度;v是網格的體積。

定義測區重力異常與第q點質量重力異常的歸一化互相關系數為式(2):

(2)

將式(1)代入式(2), 假設剩余密度為正, 可以約掉G、Δσ、v得到實測異常值Δg與質量基函數B的互相關系數計算式(3):

(3)

式(3)中:

(4)

相關系數Cq的取值范圍在-1和1之間, 其絕對值越大, 表明地表的異常是由該質量點產生的概率越大。Cq值越接近于1, 表明該點質量剩余的可能性越大; 反之Cq值越接近于-1, 表明該點質量虧損的可能性越大[7,11]。

1.2 深度加權約束

重力數據相關成像的結果在深部非常發散, 成像的位置與異常體實際的位置之間的誤差較大, 即深度方向上的分辨率較低, 故需要引入深度加權函數來提高重力數據相關成像的縱向分辨率[16,19]。在重力三維密度正則化反演中常用深度加權函數提升縱向反演分辨率, Li等[20]、Commer等[21]提出了不同形式的深度加權函數, Li等[20]提出的表達式為:

(5)

式中:z0是與異常體深度有關的常數;β是一個正數, 對于重力數據β=2;對于重力梯度數據β=3。 取z0=5,β=2時, 其函數圖像如圖1所示, 從圖中可以看出, 其權重的大小隨著深度的增加迅速減小, 在目標異常體的深度位置并沒有較高的權重; 而本文對于深度加權函數的需求是在目標深度范圍內具有較大權重, 其他深度范圍具有較小權重, 故它不適用于解決相關成像在縱向上存在的分辨率低以及深度發散的問題。

(z1=100m,z2=300m,zmax=500m)

Commer等[21]提出的空間梯度加權函數通過引入先驗信息來提高相關成像的分辨率, 其中z方向分量可以作為相關成像的深度加權函數使用[16,19],其表達式為:

(6)

式中:zmax為研究區域地下空間范圍內的最大深度;α為經驗值, 決定近地表處的權值, 為克服趨膚效應,一般取α=0.001;z1為模型頂部深度,z2為模型底部深度, 二者的取值應由研究區域的地質資料等先驗信息決定;r為縮放因子。該式與α、z1、z2、zmax、r這5個參數有關性且形式較為復雜。當先驗深度z1=100 m,z2=300 m,zmax=500 m,α分別取0.1、0.01、0.001,r分別取100、10、1、20時, 深度加權函數對于目標深度的加權效果差異很大(圖1)。

為了減少超參數對深度加權的影響, 本文提出了一種更簡潔的深度加權函數, 其表達式為:

(7)

該函數僅與z1,z2,k這3個參數有關。當先驗信息z1=100 m,z2=300 m,k分別取值為0.05、0.10、0.20、0.50時, 其函數圖像如圖2所示, 可見k的值越小, 函數值隨深度變化速率越慢,z1和z2之間的權重逐漸降低;k的值越大, 函數值隨深度變化速率越快,z1和z2之間的權重逐漸增大, 圖像越接近方波, 函數的光滑度降低。為了使深度加權函數更加光滑, 則k的值應越小; 為了保證z1和z2之間的權重更大, 則k的值應越大; 因此,k的取值應兼顧深度加權函數的光滑性和高權重性, 可通過分析其函數圖像, 選取適當的值。例如通過分析圖2中不同k值的函數圖像, 建議k的取值范圍為(0.1±0.05)時, 可以看出在先驗深度z1和z2之間具有更大權重, 而且此時函數圖像也較為光滑。

(z1=100m,z2=300m)

將上述深度加權函數應用于相關成像, 則深度加權后的相關系數可表示為:

CD=C·Wz,

(8)

式中:CD為深度加權后的相關系數;C為未加權的相關系數;Wz為深度加權函數。

1.3 邊緣特征約束

重力數據相關成像的結果在深度方向上通過引入深度加權函數可以提高成像的縱向分辨率, 然而在水平方向上存在多個異常體時, 重力數據相關成像無法準確地劃分出異常體邊界。本文通過引入邊緣特征信息來改善重力數據相關成像的橫向分辨率。邊緣特征信息一般通過垂向導數(VDR)或解析信號振幅(ASM)獲得。然而, 這類方法對于深部異常體的邊界識別結果較為模糊。因此, 通過反正切函數將邊緣特征的強、弱異常信號進行均衡放大, 實現深部邊界識別。本文選擇均衡垂向導數(VDR)和解析信號振幅(ASM)作為水平加權函數, 對比研究它們對重力數據相關成像橫向分辨率的改善效果。

反正切函數的邊界均衡識別方法, 其計算公式為:

Bf(x,y,z0)=arctan[R·f(x,y,z0)],

(9)

式中:R是均衡系數, 用來調節強、弱信號均衡放大的程度。如圖3所示,從R分別取1、5、50、100時的函數圖像可以看出,隨著R的增大函數變化的斜率增大, 對于0值附近的值的放大程度也越大, 即可實現對于強、弱信號的均衡放大。將垂向導數和解析信號振幅邊緣識別方法對相關成像進行水平加權, 則要求在識別的異常體位置區域具有較高的權重, 在非邊緣位置區域的權重迅速減小為0, 從而達到增強邊界信息的作用。

圖3 不同均衡系數的反正切均衡函數Fig.3 Arctangent balance functions with different balance coefficients

垂向導數(VDR)方法最初由Hood等[22]和Bhattacharyya[23]提出, 該方法利用零值位置來識別地質體的邊緣位置。將均衡垂向導數(取絕對值)歸一化后應用于重力數據相關成像水平加權, 其計算公式為式(10)。

解析信號振幅(ASM)又稱總梯度模量, 由Nabighian[24-25]提出, 該方法利用極大值位置來識別地質體的邊緣位置。將均衡解析信號振幅歸一化后應用于重力數據相關成像水平加權, 其計算公式為式(11)。

(10)

(11)

式中:x,y分別是觀測點坐標;z0是觀測面的高度。

本文對重力數據相關成像的加權處理可統一寫為:

CGT=C·Wz·Wh,

(12)

式中:CGT為總加權相關系數;C為未加權相關系數;Wz為深度加權函數;Wh為水平加權函數,Wh可以為nBVDR或nBASM。

2 理論模型試驗

2.1 含5%噪聲試驗

設置兩個剩余密度、大小相同的直立長方體模型, 模型參數如表1所示; 平面測網在x、y方向為0～1 000 m, 點距為10 m, 高度位于0 m; 成像范圍為0～500 m。 如圖4所示, 該模型的重力異常(圖4a)含5%高斯白噪聲, 計算出該模型的歸一化均衡垂向導數(圖4b)和歸一化均衡解析信號振幅(圖4c) , 其中均衡系數R=10.0。

表1 模型參數(含5%噪聲)Table 1 Model parameters (including 5% noise)

a—重力異常;b—歸一化均衡垂向導數(nBVDR);c—歸一化均衡解析信號振幅(nBASM)

使用Commer等[21]提出的空間梯度加權函數對相關成像進行深度加權?？臻g梯度加權函數的先驗深度信息z1=100 m,z2=300 m, 取zmax=500 m,α=0.001,r=50; 其相關成像深度加權的三維結果如圖5a所示,y=500 m處的垂向切片如圖5b所示,z=200 m處的水平切片如圖5c所示。 從圖中可以看出, Commer等[21]提出的空間梯度加權函數使成像位置集中在模型的先驗深度范圍內, 提高了相關成像的縱向分辨率; 但在水平方向上, 兩個異常體之間的邊界位置無法區分, 橫向分辨率很低。

a—深度加權約束的成像結果;b—y=500 m的垂向切片(深度加權約束);c—z=200 m的水平切片(深度加權約束);d—本文方法成像結果(Wz+VDR加權);e—y=500 m的垂向切片(Wz+VDR加權);f—z=200 m的水平切片(Wz+VDR加權);g—本文方法成像結果(Wz+ASM加權);h—y=500 m的垂向切片(Wz+ASM加權);i—z=200 m的水平切片(Wz+ASM加權)

使用本文提出的基于邊緣特征與深度加權約束的重力相關成像反演方法的三維結果如圖5d、g所示,y=500 m處的垂向切片如圖5e、h所示,z=200 m處的水平切片如圖5f、i所示。 其中縱向約束使用本文提出的深度加權函數, 先驗信息上底z1取100 m, 下底z2取300 m, 調節因子k=0.1。 橫向約束分別使用歸一化均衡垂向導數和歸一化均衡解析信號振幅。 從圖中可以看出, 在本文提出的深度加權函數的基礎上, 使用歸一化均衡垂向導數和歸一化均衡解析信號振幅的邊緣特征信息均使相關成像的橫向分辨率顯著提高; 其中, 使用歸一化均衡垂向導數的橫向分辨率更高。 基于文章篇幅的原因, 本文就不再展示該模型不含噪聲的試驗。通過對比可知含5%噪聲的結果與不含噪聲的結果基本相同, 表明相關成像的方法具有良好的抗噪性, 故后文不再展示其他含噪試驗結果。

2.2 復雜直立模型試驗

為了研究不同剩余密度成像的效果, 設置了不同剩余密度直立長方體的組合模型; 平面測網在x、y方向為0～1 000 m, 點距為10 m, 高度位于0 m, 成像范圍為0～500 m。模型參數如表2所示。該模型的重力異常為圖6a, 計算出該模型的歸一化均衡垂向導數(圖6b)和歸一化均衡解析信號振幅(圖6c), 其中均衡系數R=2.0。

表2 復雜直立模型參數Table 2 Complex upright model parameters

a—重力異常;b—歸一化均衡垂向導數(nBVDR);c—歸一化均衡解析信號振幅(nBASM)

使用Commer等[21]提出的空間梯度加權函數對相關成像進行深度加權。先驗深度信息z1=50 m,z2=300 m, 取zmax=500 m,α=0.001,r=50; 其相關成像深度加權的三維結果如圖7a所示, 在水平方向上, 剩余密度為一正一負的異常體之間的邊界位置可以區分, 但剩余密度都為正的異常體之間的邊界位置無法區分, 橫向分辨率很低。

a—深度加權約束的成像結果;b—y=300 m的垂向切片(深度加權約束);c—z=200 m的水平切片(深度加權約束);d—本文方法成像結果(Wz+VDR加權);e—y=300 m的垂向切片(Wz+VDR加權);f—z=200 m的水平切片(Wz+VDR加權);g—本文方法成像結果(Wz+ASM加權);h—y=300 m的垂向切片(Wz+ASM加權);i—z=200 m的水平切片(Wz+ASM加權)

使用本文方法的三維結果如圖7所示, 縱向約束所需先驗信息上底z1取50 m, 下底z2取300 m, 調節因子k=0.1。橫向約束分別使用歸一化均衡垂向導數和歸一化均衡解析信號振幅。經過加權后成像的縱向分辨率和橫向分辨率顯著提高; 其中, 使用歸一化均衡垂向導數比歸一化均衡解析信號振幅的縱向分辨率更高。

2.3 復雜傾斜模型試驗

為了驗證更為復雜的情況, 設置了不同剩余密度傾斜長方體和臺階的組合模型; 平面測網在x、y方向為0～1 000 m, 點距為10 m, 高度位于0 m, 成像范圍為0～500 m。模型參數如表3所示。該模型的重力異常為圖8a, 計算出該模型的歸一化均衡垂向導數(圖8b)和歸一化均衡解析信號振幅(圖8c), 其中均衡系數R=10.0。

表3 復雜傾斜模型參數Table 3 Complex tilt model parameters

a—重力異常;b—歸一化均衡垂向導數(nBVDR);c—歸一化均衡解析信號振幅(nBASM)

使用Commer等[21]提出的空間梯度加權函數對相關成像進行深度加權三維結果如圖9a所示。先驗深度信息z1=99 m,z2=430 m, 取zmax=500 m,α=0.001,r=50; 同樣在水平方向上, 剩余密度為一正一負的異常體之間的邊界位置可以區分, 但剩余密度都為正的異常體之間的邊界位置無法區分, 橫向分辨率也很低。使用本文方法的三維結果如圖9所示, 縱向約束所需先驗信息上底z1取99 m,下底z2取430 m, 調節因子k=0.1。橫向約束分別使用歸一化均衡垂向導數和歸一化均衡解析信號振幅。經過加權后成像的縱向分辨率和橫向分辨率顯著提高。

a—深度加權約束的成像結果;b—y=250 m的垂向切片(深度加權約束);c—z=250 m的水平切片(深度加權約束);d—本文方法成像結果(Wz+VDR加權);e—y=250 m的垂向切片(Wz+VDR加權);f—z=250 m的水平切片(Wz+VDR加權);g—本文方法成像結果(Wz+ASM加權);h—y=250 m的垂向切片(Wz+ASM加權);i—z=250 m的水平切片(Wz+ASM加權)

通過上述模型試驗可以發現, 使用本文提出的深度加權函數對相關成像在深度方向加權, 使成像的結果在深部更加收斂, 且使成像的位置集中在所給的先驗深度范圍內, 使相關成像的縱向分辨率顯著提高; 但深度加權函數受先驗深度信息z1、z2的影響較大, 故先驗深度信息的選取對成像結果至關重要。

由于均衡垂向導數和均衡解析信號振幅的邊緣特征方法可以增強對于異常體深部位置的識別, 故本文在水平方向上分別使用均衡垂向導數和均衡解析信號振幅對相關成像進行水平加權, 均使相關成像的橫向分辨率顯著提高; 其中, 使用均衡垂向導數加權的效果更好。

3 實際資料處理

澳大利亞的Olympic Dam (Cu-U-Au-Ag)金屬礦床是南澳大利亞太古宙元古代高勒克拉通邊緣的幾種氧化鐵—銅—金—鈾(IOCG-U)礦床中最大的礦床[26-27], 礦床位于高勒前寒武紀克拉通東部, 基底構造層為廣闊的隆起區, 區域巖石由區域變質巖、條帶狀磁鐵石英巖和花崗巖組成。礦床被深而狹窄的地塹沉積不整合覆蓋, 地塹由快速沉積的角礫巖、火山碎屑巖和長英質火山巖充填, 其上又被新元古代砂礫巖和寒武紀石灰巖層所覆蓋。角礫巖的主要成分是花崗巖碎塊和各種類型的赤鐵礦。礦床中主要的硫化礦物有黃銅礦、斑銅礦和輝銅礦。品位較高的礦帶由浸染狀輝銅礦和斑銅礦組成, 產于礦床較高部位。鈾品位較高的地區, 通常含有細脈和細粒浸染狀瀝青鈾礦。金和銀品位雖低,但與銅和鈾共生,因而具有經濟價值。在個別金礦化的地區, 稀土元素含量甚高[28-30]。

IOCG型礦床,即鐵氧化物—銅—金礦床, 指鐵氧化物含量大于20%的銅—金礦床, 其一般具有規模大、品位高、元素多、埋藏淺、易采選等特點[31-32]。由于IOCG型礦床富含鐵氧化物且蝕變范圍廣闊, 地球物理特征明顯, 尤其是重磁力特征[33], 因此,高精度的重磁勘探是尋找IOCG礦床重要的有效手段之一。由于Olympic Dam多金屬礦床受到含鐵氧化物的影響,產生大量的磁鐵礦、赤鐵礦或半生礦物,這些組合物中鐵氧化物具有不同的磁性,導致Olympic Dam多金屬礦床的磁性特征變化較大,該變化也伴隨較為顯著的密度變化和重力變化,伴隨鐵氧化物蝕變產生明顯的重力異常高[34]。

圖10顯示了高勒—克拉頓省的地質結構[35]。在高勒—克拉頓省的東段, 以黑色矩形標記的Olympic Dam多金屬礦區作為研究區域。圖11顯示了Olympic Dam多金屬礦床的地質示意圖[36]。圖12顯示了該區域的4口鉆井的巖心實測密度曲線。鉆井巖心實測密度曲線及地質剖面圖顯示高密度巖體主要分布集中在400～1 400 m范圍內。

圖10 高勒克拉頓省的基底地質圖[35]Fig.10 Basement geological map of the Gawler Craton Province [35]

圖11 超巨型Olympic Dam角礫巖群的簡化地質圖及剖面[36]Fig.11 Simplified geological map and cross-sectional subsurface structure of the supergiant Olympic Dam breccia complex [36]

圖12 鉆井巖心實測密度曲線Fig.12 Drill core measured density curve

為了獲得礦體引起的剩余重力異常, 使用二維小波多尺度分解方法[37]對該區域的布格重力異常數據進行異常的分離, 如圖13a所示為重力異常的觀測值, 圖13b為區域重力異常, 圖13c為剩余重力異常。對其剩余重力異常計算其歸一化均衡垂向導數和歸一化均衡解析信號振幅, 其中均衡系數R=2.0。

a—重力異常觀測值;b—區域重力異常;c—剩余重力異常

使用剩余重力異常數據進行相關成像, 由鉆井資料可知該高密度體的先驗深度信息大約在400～1 400 m。故在使用深度加權函數進行縱向約束時, 先驗信息z1=400 m,z2=1 400 m, 調節因子k=0.01。再分別使用歸一化均衡垂向導數和歸一化均衡解析信號振幅進行水平加權, 相關成像加權結果的三維圖、y=6 063.5 km處的垂向切片、z=750 m處的水平切片如圖14所示。成像結果的水平位置與5%Fe含量等值線基本一致, 深度位置集中在400～1 400 m范圍內; 其中均衡垂向導數比均衡解析信號振幅加權的分辨率更高; 本文方法成像的結果與前人反演的結果類似[38-39], 較好地反映了目標巖體的位置。

a—本文方法成像結果(Wz+VDR加權);b—y=6 063.5 km的垂向切片(Wz+VDR加權);c—z=750 m的水平切片(Wz+VDR加權);d—本文方法成像結果(Wz+ASM加權);e—y=6 063.5 km的垂向切片(Wz+ASM加權);f—z=750 m的水平切片(Wz+ASM加權)

4 結論

針對重力異常相關成像的結果存在深部發散、深度加權函數參數過多以及由于多個異常體的重力異常疊加導致相關成像橫向分辨率和縱向分辨率低的問題。本文提出了基于邊緣特征和深度加權約束的重力三維相關成像反演方法, 該方法引入了重力異常均衡垂向導數以及均衡解析信號振幅, 并且提出了一種更為簡潔的深度加權函數, 較好地改善了重力異常相關成像的橫向和縱向分辨率。通過模型試驗以及對澳大利亞Olympic Dam多金屬礦區實測重力資料的處理, 得到了如下結論:

1) 使用反正切函數對重力異常垂向導數和解析信號振幅識別的邊界進行均衡, 可以得到異常體深部的邊緣特征。將均衡垂向導數和均衡解析信號振幅作為重力相關成像的水平加權, 可以在橫向上明顯地區分出異常體的邊界位置, 提高了重力異常相關成像的橫向分辨率。

2) 本文提出的深度加權函數除了先驗深度信息z1和z2, 只與因子k有關, 需要給定的參數數量很少, 減少了超參數的選取對于深度加權的影響。通過給定先驗深度信息z1和z2有效地改善了相關成像深部發散的問題, 提高了重力異常相關成像的縱向分辨率。然而深度加權函數非常依賴先驗深度信息z1和z2, 可通過鉆井等其他地球物理資料得到較為準確的先驗深度信息, 從而獲得較好的深度加權相關成像結果。

3) 本文方法更適用于直立塊狀異常體的三維成像, 對于傾斜板狀異常體的成像效果有限。

4) 作為一種半定量的解釋方法, 后續可為三維密度反演提供初始模型。

致謝：感謝南澳大利亞資源信息網站提供的Olympic Dam多金屬礦區布格重力數據。