高秀鶴,熊盛青,于長春,孫思源

(1.中國自然資源航空物探遙感中心, 北京 100083; 2.自然資源部 航空地球物理與遙感地質(zhì)重點實驗室, 北京 100083；3.中國地質(zhì)大學(北京) 地球物理與信息技術(shù)學院, 北京 100029)

0 引言

重力數(shù)據(jù)反演面臨嚴重的多解性問題，高斯定理指出其本質(zhì)原因：有無數(shù)種不同的密度分布可以引起相同的重力異常[1-3]；此外，異常體引起的重力異常是連續(xù)分布的，而我們采集的信息只是某一平面上的部分離散數(shù)據(jù)，這種信息量的不足進一步加重了反演結(jié)果的多解性[4-5]。這種多解性是導致重力數(shù)據(jù)反演結(jié)果趨膚效應嚴重、分辨率低的本質(zhì)原因。利用已知的先驗地質(zhì)信息對密度進行合理約束是改善重力數(shù)據(jù)反演效果的有效手段，常用的約束方法有最小質(zhì)量約束、最小體積約束、最大光滑約束、最小慣量約束等[6-10]。另一方面，增加觀測數(shù)據(jù)的數(shù)據(jù)量也能夠改善反演效果[11-12]，增加的觀測數(shù)據(jù)可以是同種地球物理數(shù)據(jù)加密或不同平面采集，也可以是不同種類地球物理數(shù)據(jù)[13-14]。但在實際工作中，由于經(jīng)費或者測量條件限制，無法獲得更多觀測數(shù)據(jù)時，通過延拓技術(shù)獲取多高度面觀測數(shù)據(jù)加入到反演中是可行的[15]。

垂直分辨率差是重力數(shù)據(jù)反演中最為突出的問題[16]，重力數(shù)據(jù)的核函數(shù)隨深度的增加逐漸衰減是密度反演結(jié)果分辨率差的直接原因。從方法技術(shù)角度出發(fā)，引入深度加權(quán)抵消核函數(shù)衰減速率，能夠有效克服反演結(jié)果的趨膚效應，提高垂直分辨率[3,17-19]；從實測數(shù)據(jù)角度出發(fā)，引入鄰近異常源的井中觀測重力數(shù)據(jù)、井中巖石密度數(shù)據(jù)，可提高反演結(jié)果的垂直分辨率[20]。

不同高度觀測面的重力數(shù)據(jù)包含的波場信息不同[21]，因此，單獨反演效果存在差異，本文首先設(shè)計模型試驗分析了不同高度面觀測重力數(shù)據(jù)的反演效果；然后，融合反演不同高度面觀測的重力數(shù)據(jù)，與單一高度面觀測重力數(shù)據(jù)反演結(jié)果對比，分析融合反演的效果，同時，多觀測面數(shù)據(jù)不可得時，利用延拓技術(shù)獲得不同高度重力數(shù)據(jù)用于反演；最后，分析井中重力數(shù)據(jù)、井中密度數(shù)據(jù)約束下，不同高度重力數(shù)據(jù)反演效果。

1 融合反演原理

離散重力數(shù)據(jù)與離散密度塊體的關(guān)系如下：

g=Aρ，

(1)

式中：ρ是密度向量，g是某一觀測面的觀測重力數(shù)據(jù)向量，A是正演算子矩陣。反演過程中，已知觀測重力數(shù)據(jù)向量g和正演算子矩陣A，求解密度向量ρ。由于觀測數(shù)據(jù)中不可避免含有噪聲干擾，使得反演問題病態(tài)，導致解不穩(wěn)定，吉洪諾夫正則化方法作為一種常規(guī)的處理手段，可有效解決這一問題，求式(1)中ρ的反演問題可轉(zhuǎn)化為求解目標函數(shù)最優(yōu)解的問題，目標函數(shù)如下：

φ=φg+αφρ，

(2)

式中：φg是數(shù)據(jù)擬合函數(shù)，φρ是模型擬合函數(shù)，α代表正則化參數(shù)。數(shù)據(jù)擬合函數(shù)為

(3)

模型擬合函數(shù)為最小模型約束：

(4)

在模型擬合項中加入深度加權(quán)函數(shù)W，可有效克服重力反演結(jié)果中的趨膚問題。反演結(jié)果中異常體的埋深很大程度上受深度加權(quán)函數(shù)影響，弱的深度加權(quán)函數(shù)容易獲得較淺的異常體埋深，強的深度加權(quán)則容易獲得較深的異常體埋深。

引入深度加權(quán)函數(shù)后，模型擬合函數(shù)可以寫成

(5)

不同平面(高度或深度為z)重力數(shù)據(jù)融合反演時，數(shù)據(jù)擬合函數(shù)不同。觀測重力數(shù)據(jù)gz表示觀測面高度或深度為z的重力數(shù)據(jù)。將gz作為觀測數(shù)據(jù)加入到目標函數(shù)的數(shù)據(jù)擬合函數(shù)中，則數(shù)據(jù)擬合函數(shù)為

(6)

式中：Az是與gz對應的正演算子矩陣。此外，還可以在數(shù)據(jù)擬合函數(shù)中加入井中重力數(shù)據(jù)gb、井中巖石密度數(shù)據(jù)ρd，則有

(7)

其中：gb代表井中重力數(shù)據(jù)，Ab是與gb對應的正演算子，ρd代表井中巖石密度數(shù)據(jù)，Ad可以看成是與ρd對應的正演算子，它是一個稀疏矩陣，矩陣每行僅有一個非0元素，與取樣巖石位置有關(guān)，其值為1。式(7)所示的數(shù)據(jù)擬合函數(shù)加入到目標函數(shù)中，可實現(xiàn)不同高度觀測面、井中數(shù)據(jù)的融合反演，通過調(diào)整系數(shù)λ1、λ2、λb和λd的大小可調(diào)整不同高度面觀測重力數(shù)據(jù)、井中重力數(shù)據(jù)、井中巖石密度數(shù)據(jù)對目標函數(shù)的作用大?。合禂?shù)越小，相應數(shù)據(jù)對目標函數(shù)影響越小；系數(shù)越大，相應數(shù)據(jù)對目標函數(shù)影響越大；當系數(shù)為0時，相應數(shù)據(jù)對目標函數(shù)不起作用，即相應數(shù)據(jù)不參與反演過程。

2 模型試驗

2.1 模型一：單立方體模型

為了驗證不同平面觀測重力數(shù)據(jù)三維反演效果，建立一個立方體模型，如圖1a所示，立方體模型的大小為400 m×400 m×400 m，頂面埋深為300 m，中心位置為(2.5 km, 2.5 m, 0.5 km)，密度差為1 g·cm-3，背景密度為0。將地下三維空間剖分為50×50×15個大小為100 m×100 m×100 m的立方體單元，50×50個觀測點位于每個單元的中心。正演計算得到的不同高度觀測面的重力數(shù)據(jù)如圖1b、c、d、e、f所示。

利用正則化方法反演不同高度面觀測重力數(shù)據(jù)。為了對比反演效果，對于不同高度面的觀測重力數(shù)據(jù)反演，引入相同的密度上下限 (0, 1 g·cm-3)和相同的模型加權(quán)函數(shù)W=(ATA)1/2，其中，A表示地面重力數(shù)據(jù)的正演算子。利用共軛梯度算法迭代求解目標函數(shù)最優(yōu)解，得到的反演結(jié)果切片如圖2所示。觀測面高度越高，相應的反演結(jié)果越發(fā)散；觀測面高度越低，相應的反演結(jié)果越收斂，且反演結(jié)果與理論模型的分布范圍、密度值大小越接近。

不同高度面的觀測重力數(shù)據(jù)融合反演結(jié)果如圖3所示。圖3a與圖2c相比，由于地面以下300 m觀測數(shù)據(jù)的加入，密度分布更收斂，密度值更接近理論密度值；圖3a與圖2f相比，由于地面觀測數(shù)據(jù)的加入，壓制了異常體兩側(cè)的假異常。這初步驗證了不同高度面觀測重力數(shù)據(jù)融合反演的良好效果。相比于圖2c，圖3a顯然與圖2f更接近，說明近源平面數(shù)據(jù)對反演結(jié)果影響更大。圖3b與圖2b、d相比，密度分布更收斂，密度值更接近理論密度值，與圖2f相比，壓制了異常體兩側(cè)的假異常。這再次驗證了不同高度面觀測重力數(shù)據(jù)融合反演的良好效果。圖3c與圖3b相比，增加了兩個面的觀測數(shù)據(jù)，反演效果并未見明顯提高。這說明為保證計算效率，需選擇適當、適量平面的觀測重力數(shù)據(jù)進行融合反演，過多的數(shù)據(jù)將導致計算量增大，但反演效果并沒有明顯的提高。

僅有一個觀測面的重力數(shù)據(jù)時，可通過延拓技術(shù)得到多個不同高度面的延拓重力數(shù)據(jù)。使用延拓數(shù)據(jù)反演需注意幾個問題:一是為提高反演效果，使用向下延拓重力數(shù)據(jù)進行反演，原因是：單個平面重力數(shù)據(jù)反演驗證，觀測面高度越低，反演效果越好；二是為保證反演可靠性，下延平面應高于異常體頂面，原因是：穿過或低于異常體的延拓數(shù)據(jù)不準確，不適合用于反演；三是為了擬合原始觀測重力數(shù)據(jù)，使用延拓數(shù)據(jù)與原始觀測數(shù)據(jù)融合反演，原因是：向下延拓得到的重力數(shù)據(jù)存在誤差，單獨反演延拓數(shù)據(jù)將導致反演結(jié)果可能無法擬合原始觀測數(shù)據(jù)。

圖1 立方體模型及正演重力數(shù)據(jù)Fig.1 Cube model and its forward gravity data

圖2 不同高度面觀測重力數(shù)據(jù)單獨反演結(jié)果的豎直切片(y=2.5 km) Fig.2 Vertical slices of inversion results of gravity data at different heights

圖3 不同高度面觀測重力數(shù)據(jù)融合反演結(jié)果的豎直切片(y=2.5 km)Fig.3 Vertical slices of inversion results of joint gravity data at different heights(y=2.5 km)

觀測重力數(shù)據(jù)及其延拓重力數(shù)據(jù)融合反演結(jié)果如圖4所示，圖4a與圖3a相比，密度分布較發(fā)散，密度值較小，這是由于向下延拓數(shù)據(jù)中存在誤差，說明延拓數(shù)據(jù)不能取代觀測數(shù)據(jù)；與圖2c相比，密度分布更收斂，密度值更接近真實密度值，反演效果明顯提高，這說明，在多平面數(shù)據(jù)不可得時，延拓數(shù)據(jù)的加入可提高反演效果。圖4b與圖3b相比，密度分布較發(fā)散；但與圖2b相比，反演效果明顯提高，再次驗證了加入延拓數(shù)據(jù)的意義。圖4c與圖4b對比，反演效果提高并不明顯。因此，為保證計算效率，選擇適當、適量的延拓數(shù)據(jù)加入反演是至關(guān)重要的，過多的延拓數(shù)據(jù)導致計算量增大，但反演效果并沒有明顯提高。

圖4 觀測重力數(shù)據(jù)及其延拓重力數(shù)據(jù)融合反演結(jié)果的豎直切片(y=2.5 km) Fig.4 Inverted density models using the joint observed gravity and continuation gravity data(y=2.5 km)

利用單立方體模型驗證加入井數(shù)據(jù)約束對反演效果的影響。模擬設(shè)計3口豎直井，3口豎直井分別位于立方體中心、立方體外100 m、立方體外300 m，其坐標分別為(2.5 km,2.5 km)、(2.8 km,2.5 km)、(3.0 km,2.5 km)。井中觀測重力數(shù)據(jù)、巖石密度數(shù)據(jù)如圖5所示。在地面、地面以下300 m重力數(shù)據(jù)融合反演中，分別加入3口井數(shù)據(jù)做約束，得到的反演結(jié)果如圖6所示。與圖3a相比，反演效果明顯提高，這說明加入位于立方體中心的井數(shù)據(jù)做約束可有效改善反演效果；加入位于立方體外100 m的井數(shù)據(jù)做約束對鄰近井位的密度分布有一定改善作用，對整體反演結(jié)果影響較?。患尤胛挥诹⒎襟w外300 m的井數(shù)據(jù)做約束對反演結(jié)果的影響很小。這說明，井口穿過異常體，能有效改善反演效果；井口距離異常體越遠，對反演效果的有益影響越小。

圖5 不同位置井的井中重力和巖石密度Fig.5 Well gravity and well density at different locations

圖6 地面、地面以下300 m重力數(shù)據(jù)和井中重力、巖石密度數(shù)據(jù)融合反演結(jié)果(白色實線代表鉆井位置)Fig.6 Recovered density models in the constraint of well data(the solid white line represents the drilling location)

2.2 模型二：雙長方體模型

為進一步驗證不同高度重力數(shù)據(jù)和井數(shù)據(jù)融合反演的效果，設(shè)置不同埋深的雙長方體模型，左側(cè)長方體頂面埋深300 m，中心位置坐標為(550 m,1 050 m,500 m)，大小為300 m×300 m×400 m，密度為1 g·cm-3，右側(cè)長方體頂面埋深為500 m，中心位置坐標為(1 350 m,1 050 m,700 m)，其大小、密度與左側(cè)長方體相同。地下空間剖分為38×42×15個大小為50 m×50 m×100 m的立方體單元，38×42個地面觀測數(shù)據(jù)采樣間距為50 m，井中觀測數(shù)據(jù)采樣間距為50 m。加入3%的高斯噪聲的正演理論地面重力數(shù)據(jù)如圖7a所示，將地面重力數(shù)據(jù)向下延拓300 m 得到的延拓重力數(shù)據(jù)如圖7b所示。模擬兩口豎直井：井A、井B，其水平位置分別為(550 m,1 050 m)、(1 350 m,1 050 m)，即位于兩個長方體中心，填加3%高斯噪聲的井中重力、巖石密度數(shù)據(jù)如圖8所示。

地面重力數(shù)據(jù)單獨反演結(jié)果如圖9b所示，左側(cè)埋深較淺的異常體中心位置識別準確，右側(cè)埋深較深的異常體未被有效識別，且無法區(qū)分兩個異常體，反演結(jié)果整體較為發(fā)散，分辨率較低。加入向下延拓300 m重力數(shù)據(jù)與地面重力數(shù)據(jù)融合反演結(jié)果如圖9c，與圖9b相比，左側(cè)異常體中心位置定位準確，且異常體邊界清晰，密度值的最大值達到真實值，右側(cè)異常體反演效果仍不理想，且由于下延數(shù)據(jù)中噪聲很大，導致反演結(jié)果中也出現(xiàn)了假異常。再加入兩口井的井中重力、巖石密度數(shù)據(jù)約束，反演結(jié)果如圖9d，準確定位了兩個不同埋深的異常體的中心位置，邊界清晰，分辨率高，兩個異常體的最大異常值都達到了理論值，與圖9a所示的理論模型剖面匹配度較高。

圖7 雙長方體模擬重力數(shù)據(jù)Fig.7 Gravity data based on double prisms

圖8 含3%高斯噪聲的模擬井中觀測數(shù)據(jù)Fig.8 The simulated data observed in the wells

圖9 雙立方體模型反演結(jié)果 (白色實線代表鉆井位置)Fig.9 Recovered double-prism models(the solid white line represents the drilling location)

3 結(jié)論

1) 單一觀測面重力數(shù)據(jù)三維反演效果與測量面高度有關(guān)，觀測面高度越低(觀測面距離異常源越近)，反演效果越好。

2) 與單一觀測面重力數(shù)據(jù)反演相比，多個高度面觀測重力數(shù)據(jù)融合反演的效果可得到改善，其中，近源高度面數(shù)據(jù)比遠源高度面數(shù)據(jù)對反演結(jié)果影響大。

3) 過多個觀測面的重力數(shù)據(jù)對反演效果改善不明顯，因此，為保證計算效率，應選擇適當、適量的觀測面的重力數(shù)據(jù)進行融合反演。

4) 當多個高度面的觀測重力數(shù)據(jù)不可得時，通過向下延拓技術(shù)可得到多個平面的延拓重力數(shù)據(jù)(延拓平面應高于異常體頂面)，觀測重力數(shù)據(jù)與向下延拓重力數(shù)據(jù)融合反演，也可改善反演效果。

5) 井中重力、巖石密度數(shù)據(jù)約束，可一定程度改善反演效果，改善程度受鉆井位置影響明顯，穿過異常體的井數(shù)據(jù)可明顯改善反演效果，鉆井位置距離異常體越遠，對反演效果影響越小。