劉 潔 張建中*② 江 麗 萬 麗 胡加山

(①中國海洋大學(xué)海洋地球科學(xué)學(xué)院海底科學(xué)與探測技術(shù)教育部重點實驗室,山東青島 266100；②青島海洋科學(xué)與技術(shù)國家實驗室海洋礦產(chǎn)資源評價與探測技術(shù)功能實驗室，山東青島 266061；③中國石化勝利油田分公司勘探開發(fā)研究院，山東東營 257015)

0 引言

重力勘探是通過探測地下介質(zhì)的密度分布研究地層結(jié)構(gòu)或目標(biāo)體的重要地球物理手段之一。目前，重力數(shù)據(jù)反演主要包括兩方面的內(nèi)容。一是界面(幾何參數(shù))反演[1-3]，通常假設(shè)地層密度已知，通過反演觀測重力數(shù)據(jù)確定某個密度體的邊界或地層界面深度，廣泛應(yīng)用于求取沉積層的基底深度[4]或深部不連續(xù)界面等；重力反演的另外一個內(nèi)容就是地層或目標(biāo)體的密度(物性參數(shù))[5-6]，通常將地下介質(zhì)網(wǎng)格剖分成許多矩形單元，每個矩形單元的密度值設(shè)為一個常數(shù)，利用觀測的重力異常值，通過反演獲得密度分布，據(jù)此識別斷裂帶、深部巖漿體或者描繪淺部巖脈和礦體展布等。

現(xiàn)有的密度反演方法通常是在水平方向和垂向上剖分網(wǎng)格單元，假設(shè)每個單元為常密度，通過反演每個單元的密度值獲得地下空間的密度分布。較為典型的方法是L2模約束反演[7]，該方法通過數(shù)據(jù)項匹配和平滑約束獲得具有地質(zhì)意義的解，再通過密度全正或全負(fù)約束壓縮解空間，但反演模型結(jié)果較發(fā)散，尤其是當(dāng)正、負(fù)密度異常同時存在時，此現(xiàn)象更嚴(yán)重。為了更好地研究異常體邊界，把稀疏約束、最小支撐約束等加入反演，發(fā)展了聚焦反演法[8-9]，通過數(shù)學(xué)手段實現(xiàn)模型體積最小化。然而，上述方法均基于橫向和垂向上同時進(jìn)行網(wǎng)格剖分。若所劃分網(wǎng)格單元數(shù)目過少，則反演結(jié)果必然不夠精細(xì)；若垂向網(wǎng)格單元劃分?jǐn)?shù)目過多，反演未知量的數(shù)目和占用內(nèi)存也會成倍增加，多解性問題也更嚴(yán)重。為此，姚長利等[10-11]研究了重磁異常快速計算和有效存儲技術(shù)。

近年來，變密度方法逐漸受到學(xué)者們的關(guān)注。該方法假設(shè)密度隨深度變化，地層密度的垂向變化可以表示為深度的函數(shù)，簡單密度函數(shù)產(chǎn)生的重力異常可直接推導(dǎo)出相應(yīng)的解析表達(dá)式[12-14]，無需進(jìn)行網(wǎng)格剖分即可計算變密度體的重力異常。鑒于變密度體重力異常解析表達(dá)精度高、耗時少等優(yōu)點，低階的密度函數(shù)形式(如二次函數(shù)、指數(shù)函數(shù)、拋物線函數(shù)、雙曲線函數(shù)等)已廣泛應(yīng)用于密度界面的反演[15-16]。

但是，變密度遲遲未能在反演中發(fā)揮作用，原因是低階的密度函數(shù)變化趨勢相近且單調(diào)，無法描述地下復(fù)雜的密度變化，而且對應(yīng)的重力異常解析表達(dá)式各不相同，這為重力反演帶來諸多不便。針對上述問題，Zhang等[17]利用任意階多項式表示地下密度異常分布，且推導(dǎo)了計算二維變密度體重力異常的解析表達(dá)式。近些年，二維和三維變密度體重力異常正演計算得到進(jìn)一步完善和發(fā)展[18-21],為任意不規(guī)則變密度體的重力異常正演和垂向無網(wǎng)格剖分的變密度反演奠定了堅實基礎(chǔ)。

本文在變密度體重力異常正演的基礎(chǔ)上，提出一種基于高階多項式密度函數(shù)的重力異常反演方法，該方法無需在垂向上進(jìn)行網(wǎng)格剖分，通過反演多項式密度函數(shù)的系數(shù)就可以確定地下物質(zhì)的密度分布。理論模型測試和實際資料應(yīng)用結(jié)果證明了方法的有效性。

1 方法原理

1.1 變密度直立矩形二度體重力異常正演解析式

與常規(guī)密度反演的網(wǎng)格剖分方法不同，將地下空間劃分為R個橫向排列的長條矩形單元(圖1)，將第i個矩形單元的密度差ρi表示為隨深度z變化的N階多項式函數(shù)

(1)

式中ai,j是第i個矩形單元多項式的系數(shù)。

如圖1所示，假設(shè)第i個矩形單元的四個頂點坐標(biāo)分別為A(xi，1,zi，1),B(xi，1,zi，2),C(xi，2,zi，2)和D(xi，2,zi，1)，該矩形單元的橫截面用Si表示，則該單元在測點P(xp,zp)處產(chǎn)生的重力異常為

圖1 地下介質(zhì)網(wǎng)格劃分示意圖

(2)

式中G是萬有引力常數(shù)。將式(1)代入式(2)，得

Δgi(xp,zp)

(3)

根據(jù)斯托克斯定理和文獻(xiàn)[17-18,22]，可將式(3)中的面積分轉(zhuǎn)換為矩形單元四條邊沿逆時針方向的閉環(huán)積分，整理可得

(4)

Ei(j,m,k)=

(5)

其中

(6)

(7)

式(4)為一個矩形單元在測點P處產(chǎn)生的重力異常，該測點測得的重力異常為所有直立矩形單元產(chǎn)生重力異常的疊加，即

(8)

1.2 多項式密度函數(shù)系數(shù)的約束反演

1.2.1 目標(biāo)函數(shù)

假設(shè)Δg=[Δgx1,z1,Δgx2,z2,…,ΔgxM,zM]T為在M個測點位置處實測的重力異常值，利用式(8)可以得到下列方程組

FM×(N+1)Ra(N+1)R×1=ΔgM×1

(9)

式中

(10)

為重力核函數(shù)矩陣；

a(N+1)R×1=[a1,0,…,a1,N,…,aR,0,…,aR,N]T

(11)

為由R個矩形單元密度函數(shù)多項式系數(shù)組成的矩陣。

與其他位場反演一樣，重力反演具有很強(qiáng)的多解性，降低其多解性的有效方法之一就是施加各種約束。為此，結(jié)合密度差上下限約束、平滑約束、深度加權(quán)和聚焦約束等建立目標(biāo)函數(shù)，以期在盡量少的先驗信息和人為干預(yù)下產(chǎn)生符合地質(zhì)規(guī)律的解。目標(biāo)函數(shù)為

(12)

式中：λ1、λ2、λ3分別為各項的權(quán)重系數(shù)；Wd=diag{1/σ1,1/σ2,…,1/σM}為數(shù)據(jù)權(quán)重矩陣，其中σi為第i個測點觀測數(shù)據(jù)的標(biāo)準(zhǔn)偏差；第一項為數(shù)據(jù)殘差項；第二項為先驗?zāi)Ｐ图s束項，在先驗密度差模型信息未知時，該項用零模型約束；第三項為橫向平滑約束項；第四項為垂向平滑約束項。式中其他參數(shù)及約束施加具體描述如下。

1.2.2 上、下限約束

一般來說，密度差分布在一定區(qū)間范圍內(nèi)。密度差的上下限約束有利于減少反演多解性。因此構(gòu)建對角分塊矩陣P=diag{Pi}(i=1,2,…,R)，其中

(13)

這里的z1、z2、…、zK為K個不同的深度點。那么，式(12)中的Pa即可表示每個矩形單元不同深度點的密度差。令ρmin≤Pa≤ρmax，其中ρmin和ρmax分別為密度差的上、下限， 從而實現(xiàn)對密度差絕對幅值的約束。

1.2.3 平滑約束

平滑約束，即幅值相對約束，有利于產(chǎn)生最小結(jié)構(gòu)解。與平滑相關(guān)的差分矩陣廣泛應(yīng)用于反演問題中。式(12)中，H為密度差橫向二階差分矩陣，具體形式為

(14)

式(12)中V為密度差縱向二階差分矩陣，具體形式為

V=diag{Vi}i=1,2,…，R

(15)

式中

1.2.4 深度加權(quán)

深度加權(quán)可以減弱重力核函數(shù)隨深度增加而衰減的影響。式(12)中Cd為深度加權(quán)矩陣

Cd=diag{cd(i,k)}i=1,2,…,R，k=1,2,…,K

(16)

其元素值與深度有關(guān)，表示為

(17)

式中：β一般取1.5～2；z0值大小取決于目標(biāo)體的中心埋深與重力異常幅值之間的關(guān)系，重力異常幅值越大，目標(biāo)體埋深越淺，則z0取值越大。

1.2.5 聚焦約束

平滑約束產(chǎn)生的結(jié)果往往較為發(fā)散，聚焦約束則可以壓縮模型體積，突出異常體的邊界特征。式(12)中Cm為聚焦約束矩陣，聚焦可通過多次反演改變模型空間內(nèi)部的權(quán)重實現(xiàn)

Cm=diag{cm(i,k)}i=1,2,…,R，k=1,2,…,K

(18)

利用最小支撐泛函,Cm中的元素值表示為

(19)

式中：ρi,k表示某次反演獲得的第i個長條矩形在zk深度處的密度值；γ為一個防止分母無意義的很小值。

關(guān)于多項式最大階數(shù)的選擇，一般地，多項式最高階數(shù)越大，能夠模擬的模型復(fù)雜程度越高，Jiang等[19]指出8階多項式便可以較好地表示密度突變的情況。為兼顧計算效率，取N在8～11范圍內(nèi)。在利用多項式函數(shù)模擬密度變化時，多項式曲線在階次較高時容易出現(xiàn)振蕩現(xiàn)象，通過式(12)中第二項施加的先驗?zāi)Ｐ图s束和上下限約束可較好地解決這一問題。

由于通過反演多項式系數(shù)可間接獲得地下密度分布，約束的施加需要取一系列深度點并通過多項式函數(shù)將其轉(zhuǎn)換為密度值。深度取樣間隔不可過大，過大則缺乏約束能力，正如同網(wǎng)格剖分法中單元劃分過于稀疏，從而降低反演的分辨率。與網(wǎng)格剖分法不同的是，本文方法可適當(dāng)增加約束點數(shù)而不增加反演未知參量的數(shù)量。

2 理論模型試驗

為了驗證本文提出的密度多項式系數(shù)反演方法，設(shè)計了3個模型進(jìn)行試驗。研究區(qū)范圍均設(shè)置為8km×3km，理論重力異常數(shù)據(jù)利用式(8)計算，并添加均值為0、標(biāo)準(zhǔn)差為0.01mGal的高斯白噪聲。使用本文方法將地下等間距劃分為60個長條矩形，每個矩形的縱向長度均為3km；將密度隨深度變化用9階多項式函數(shù)表示，即N=9； 設(shè)置密度差值約束的上、下限分別為500kg/m3和-500kg/m3；設(shè)置權(quán)重參數(shù)λ2=λ3=0.01，初次反演時取λ1=0.01，后續(xù)迭代中取λ1=1； 迭代次數(shù)一般取5～10； 初次反演時Cm取為單位矩陣，即無聚焦約束，相當(dāng)于僅L2模約束反演；后續(xù)迭代中Cm由式(19)不斷更新，進(jìn)而得到最終的反演結(jié)果。

2.1 模型Ⅰ：中心埋深相同、水平間距不同的雙密度體模型

模型Ⅰ(圖2左)截面由兩個邊長為1 km的正方體組成，中心埋深均為1.5 km，兩異常體密度差均為300kg/m3。為了測試算法在水平方向上分辨兩異常體的能力，設(shè)置兩密度體的水平間距d分別為3、2、1、0km。模擬測線沿水平方向，在0～8km等間距布設(shè)120個測點，理論重力異常如圖3所示。

圖2右顯示了不同水平間距下模型的迭代反演結(jié)果。由圖可見，在水平間距較大時，本文反演方法能較好地分辨兩異常體(圖2a)；隨著水平間距減小，兩密度體產(chǎn)生的重力異常疊加效應(yīng)增強(qiáng)，此時，兩異常體仍然能被分辨出來(圖2b和圖2c)；直到兩異常體合并在一起(圖2d)，反演結(jié)果能正確反映異常體的規(guī)模和邊界。圖3中灰色點線是不同水平間距下利用反演模型計算的重力異常值，可見與理論值吻合很好。

2.2 模型Ⅱ：中心埋深相同、形狀不同的正負(fù)密度體模型

模型Ⅱ由一個邊長為1km的正方體和一傾斜的截面為平行四邊形密度異常體組成(圖4a)，兩個異常體的中心埋深均為1.5km，兩異常體密度差分別為300kg/m3和-300kg/m3。

首先，利用本文方法在無聚焦約束下進(jìn)行初次反演(圖4c)，可見該結(jié)果呈現(xiàn)出一正一負(fù)的兩個密度異常區(qū)域，但異常范圍較實際范圍稍大，且異常幅值比實際值偏小，解較發(fā)散，傾斜異常體的傾向無法識別。另外，由于正負(fù)異常同時存在，反演結(jié)果表現(xiàn)為正負(fù)極性之間平緩過渡。將模型劃分為60×30個常密度單元，采用L2模約束進(jìn)行反演，其結(jié)果(本文未展示)與圖4c基本一致。多項式系數(shù)反演方法僅需要反演60×10=600個系數(shù)，而網(wǎng)格剖分法需要反演1800個密度值。采用本文方法迭代反演的結(jié)果如圖4d所示，可見反演結(jié)果的異常體位置和形狀都與模型較吻合，且異常體邊界更加清晰，傾向也更易識別。

圖2 模型Ⅰ兩個異常體在不同水平間距d時的真實模型(左)及反演結(jié)果(右)

圖3 模型Ⅰ兩個異常體在不同水平間距

2.3 模型Ⅲ：形狀、中心埋深均不同的正、負(fù)密度體模型

模型Ⅲ由規(guī)模、埋深均不相同的兩個正、兩個負(fù)密度異常體組成(圖5a)，異常體參數(shù)見表1。該模型產(chǎn)生的重力異常如圖5b所示，由圖可見，淺部異常體產(chǎn)生的重力異常更為高陡，而埋深較深的密度體產(chǎn)生的重力異常較為寬緩。圖5c是本文方法僅L2模約束下的初次反演結(jié)果。由圖可見，與模型二模擬結(jié)果類似，密度模型較為彌散，且重力異常幅值較真實值偏小，淺部異常體產(chǎn)生的異常較深部異常體誤差稍小一些；深度越大，異常彌散范圍越大，僅能大體判斷異常體的中心位置；圖5d為本文方法反演結(jié)果，可以看出淺部與深部異常體的重力異常在位置和形狀上均與真實模型吻合較好，異常體規(guī)模和邊界更加突出。圖5a中利用反演結(jié)果(圖5d)正演計算的重力異常值與理論重力異常值吻合較好，兩者誤差小于0.5%，表明了本文反演方法在復(fù)雜多源情況下是可行的。

表1 模型Ⅲ密度異常體位置及密度差參數(shù)

圖4 模型Ⅱ正、反演結(jié)果

圖5 模型Ⅲ正、反演結(jié)果

3 實際資料應(yīng)用

研究區(qū)位于渤海灣盆地濟(jì)陽坳陷，研究目標(biāo)是該區(qū)的古潛山。渤海灣盆地歷經(jīng)印支期褶皺隆升、燕山期斷塊抬升和喜山期拉張斷塊改造及覆蓋掩埋等三個重要的演化階段[23]。研究區(qū)地勢平坦，由地形引起的重力異?？珊雎圆挥嫛|營運(yùn)動之后，盆地漸趨穩(wěn)定，地形逐漸夷平，淺部地層厚度基本一致，且構(gòu)造運(yùn)動不明顯，內(nèi)部密度分布均勻。因此，第四系底界面及新近系與古近系之間的界面對重力場形態(tài)的影響較小[24]。古潛山頂面是該區(qū)最主要的物性界面，也是本文主要的研究對象。該區(qū)潛山主要由下古生界碳酸鹽巖和太古界變質(zhì)巖組成，表現(xiàn)為高密度、高阻抗的特征。區(qū)內(nèi)局部凹陷主要由中、新生代(部分地區(qū)含石炭、二疊系)碎屑巖組成，相對于潛山呈現(xiàn)低密度特征。潛山與上覆低密度體的密度差最大可達(dá)180kg/m3，這為利用重力資料研究潛山構(gòu)造提供了物性基礎(chǔ)。

圖6是過研究區(qū)的一條南北向地震剖面，從圖中可以比較容易地識別地層分界面和正斷層等。古近紀(jì)該區(qū)經(jīng)歷了斷陷、斷拗和拗陷等階段，在地震剖面上可觀察到左側(cè)凹陷至凸起(潛山)的過渡帶上廣泛發(fā)育超覆尖滅和錯斷。潛山頂界面為速度和密度的突變面，由于地層的高速屏蔽效應(yīng)，地震波難對潛山內(nèi)部進(jìn)行成像。右側(cè)潛山的頂界面成像較為清晰，而左側(cè)潛山右翼頂界面成像模糊。僅憑現(xiàn)有地震資料，兩潛山之間是否存在低密度溝槽尚存在爭議。

由圖7a所示的經(jīng)區(qū)域背景場分離后獲得的剩余重力異常(實線)可見，剩余重力異常曲線可以突出潛山頂面起伏及沉積蓋層的厚度變化。重力低對應(yīng)基底斷陷或坳陷，蓋層增厚；重力高對應(yīng)基底隆起，蓋層減薄。重力異常曲線明顯出現(xiàn)兩個重力高值區(qū)，在約9.5km和19.0km處分別存在一重力極大值點，該重力高解釋為地下潛山的響應(yīng)。采用本文方法進(jìn)行反演，將地下等間距劃分為70個長條矩形，矩形單元的最大深度為4km，將密度隨深度變化用10階多項式函數(shù)表示，反演各長條矩形多項式密度函數(shù)的系數(shù)以獲得地下密度分布。

圖7b為利用本文方法反演得到的密度差剖面，結(jié)果顯示研究區(qū)存在明顯的“兩正兩負(fù)”的密度異常區(qū)，正密度區(qū)域指示了潛山的分布位置，而負(fù)密度區(qū)域則指示了凹陷或溝槽的存在。對比發(fā)現(xiàn)反演的密度差剖面(圖7b)與地震剖面(圖6)有良好的對應(yīng)關(guān)系： 正密度體規(guī)模和邊界與地震潛山成像及潛山頂界面吻合；兩潛山之間的溝槽在反演結(jié)果也有體現(xiàn)，呈“兩凸夾一凹”的特征。對比利用反演密度模型計算的重力異常(圖7a)與實測的剩余重力異常，兩者相對誤差小于1%。基于本文方法的重力反演結(jié)果彌補(bǔ)了該區(qū)現(xiàn)有地震資料的不足，證實了溝槽的存在。

圖6 研究區(qū)重力測線(圖7)對應(yīng)的疊后地震剖面

圖7 實際資料反演

4 結(jié)論

本文提出了一種基于高階多項式密度函數(shù)的重力異常反演方法，該方法具有以下特征：

(1)地下復(fù)雜的密度變化可利用多項式函數(shù)表示，變密度體產(chǎn)生的重力異常采用解析式表達(dá)，計算效率較高;

(2)將地下劃分為豎直并排的矩形單元，只需反演各個矩形單元的多項式系數(shù)，無需在垂向上剖分網(wǎng)格，求解未知量相對較少，一定程度上降低了反演的多解性;

(3)構(gòu)建的目標(biāo)函數(shù)靈活，可以方便地施加約束條件和其他地質(zhì)、地球物理先驗信息;

(4)將隨深度連續(xù)變化的密度多項式函數(shù)與多種約束相結(jié)合，既能研究異常體的空間位置和規(guī)模，又能突出異常體邊界信息，得到符合地質(zhì)規(guī)律的模型解。

該方法成功應(yīng)用于渤海灣盆地濟(jì)陽坳陷區(qū)古潛山的識別，重力反演結(jié)果證實了該區(qū)潛山與溝槽呈高低相間的展布特征。本文僅就基于高階多項式密度函數(shù)的二維重力反演的原理方法做了詳細(xì)闡述，并進(jìn)行初步應(yīng)用，但所述反演思路和約束策略同樣適用于三維情況。

需要指出的是，由于多項式函數(shù)是連續(xù)平滑的，本文密度表示方法可能更適合描述密度連續(xù)變化的地層。在反演局部密度體時，由于多解性的存在，施加合理的約束是必要的。本文只對高階多項式密度函數(shù)展開研究，對低階多項式討論較少，而低階多項式可以模擬密度的整體趨勢，不同階次信息的結(jié)合和充分利用是下一步需要研究和完善的方向。