基于高階多項(xiàng)式密度函數(shù)的重力反演

劉 潔 張建中*② 江 麗 萬(wàn) 麗 胡加山

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

0 引言

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

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

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

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

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

1 方法原理

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

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

(1)

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

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

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

(2)

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

Δgi(xp,zp)

(3)

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

(4)

Ei(j,m,k)=

(5)

其中

(6)

(7)

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

(8)

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

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

假設(shè)Δg=[Δgx1,z1,Δgx2,z2,…,ΔgxM,zM]T為在M個(gè)測(cè)點(diǎn)位置處實(shí)測(cè)的重力異常值，利用式(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個(gè)矩形單元密度函數(shù)多項(xiàng)式系數(shù)組成的矩陣。

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

(12)

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

1.2.2 上、下限約束

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

(13)

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

1.2.3 平滑約束

平滑約束，即幅值相對(duì)約束，有利于產(chǎn)生最小結(jié)構(gòu)解。與平滑相關(guān)的差分矩陣廣泛應(yīng)用于反演問(wèn)題中。式(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為聚焦約束矩陣，聚焦可通過(guò)多次反演改變模型空間內(nèi)部的權(quán)重實(shí)現(xiàn)

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

(18)

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

(19)

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

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

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

2 理論模型試驗(yàn)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖7 實(shí)際資料反演

4 結(jié)論

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

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

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

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

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

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

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