劉富強(qiáng)

（新疆維吾爾自治區(qū)有色地質(zhì)勘查局物探隊(duì)烏魯木齊830011）

利用重力數(shù)據(jù)和深度控制同時(shí)估算三維基底深度和密度差

劉富強(qiáng)

（新疆維吾爾自治區(qū)有色地質(zhì)勘查局物探隊(duì)烏魯木齊830011）

本文提出了一種重力反演方法，它結(jié)合一些已知的基底深度，同時(shí)估算了沉積盆地的三維基底起伏和沉積包裹體的密度差隨深度按拋物線衰減的參數(shù)。沉積物被并列在水平方向的三維垂向棱柱體網(wǎng)格所近似。棱柱體的厚度代表基底的深度，它是從重力數(shù)據(jù)中被估算的參數(shù)。為了估算密度差隨深度拋物線衰減的參數(shù)，并得到穩(wěn)定的基底深度的估算，我們對(duì)基底深度強(qiáng)行平滑并近似鉆孔的估計(jì)深度和已知深度。將該方法應(yīng)用于復(fù)雜的兩個(gè)沉積截面的三維基底起伏的組合數(shù)據(jù)中，能清楚地用拋物線定律描述密度差隨深度的變化。它的結(jié)果良好的估算密度差隨深度拋物線衰減和基底起伏的真實(shí)參數(shù)。

重力數(shù)據(jù)深度控制三維基底深度密度差

1 引言

大多數(shù)的重力反演方法估算沉積盆地的基底起伏可以劃分為兩類：第一類假定在沉積層和基底之間的一個(gè)密度差是常量（Bott,1960;Oldenburg, 1974;Le?oetal.,1996;Barbosaetal.,1997）；第二類是假定由壓實(shí)引起的沉積物的密度差隨深度變化。一些基底深度的重力數(shù)據(jù)反演的新發(fā)展假定密度差隨深度是單調(diào)衰減的，如指數(shù)定律（Cordell,1973），二次方程定律（Rao,1990;Gallardo-Delgado,2003），拋物線定律（_Chakravarthi和Sundararajan,2004），三次多項(xiàng)式定律（García-Abdeslem,2005）和雙曲線定律（Silvaetal.,2006）。

用拋物線定律方法（如Chakravarthi和Sundararajan,2007）估計(jì)區(qū)域重力背景和沉積盆地的三維基底起伏深度，需假定密度差隨深度衰減遵循拋物線定律。滿足兩個(gè)假設(shè)，一是假設(shè)各個(gè)重力觀測(cè)站的重力異常是由重力觀測(cè)站以下的密度差變化符合拋物線定律的水平無限平板產(chǎn)生的，這個(gè)初始假設(shè)半定量地遵循幾何學(xué)基底起伏。二是假設(shè)：基底起伏很簡(jiǎn)單，基底的深度很淺，密度差隨深度的衰減率很低。第一個(gè)假設(shè)在地質(zhì)學(xué)角度上是很重要的；當(dāng)結(jié)合第二個(gè)假設(shè)時(shí)，它會(huì)得到一個(gè)穩(wěn)定的解。

通過假設(shè)密度差隨深度衰減符合拋物線定律來提出一個(gè)三維基底深度重力反演。為了確定定義的密度差隨深度拋物線衰減的參數(shù)，通過平滑基底深度和已知深度，估算一個(gè)穩(wěn)定的基底起伏。重復(fù)這些步驟是為了得到這些參數(shù)的不同值，并選擇那些與鉆孔資料的基底起伏相關(guān)性較好的參數(shù)值。通過正則化參數(shù)的穩(wěn)定函數(shù)最小化，使目標(biāo)函數(shù)成為全局凸函數(shù)來獲得基底深度的估算。結(jié)果表明，該方法不依賴于基底深度的初始值，并且可以用來估算那些地質(zhì)構(gòu)造復(fù)雜的基底起伏。

2 研究方法

假設(shè)一個(gè)由均勻基底巖石和不均勻沉積物所組成的沉積盆地，基底和沉積物之間的密度差隨深度z減小符合拋物線定律（Raoetal.,1994）：

式中，Δρo是地表面的密度差；α是通過深度控制密度差梯度的因子。

假設(shè)一個(gè)均勻基底中分離出來的任意界面，為了估計(jì)該界面的起伏，在x-y空間選擇一個(gè)包含盆地的整個(gè)上表面的有限區(qū)域，將它沿x和y方向離散成一個(gè)mx×my的三維垂向并列分布的棱柱體的網(wǎng)格（圖1）。每個(gè)棱柱體的頂部與地表面一致，并且所有棱柱體沿x方向的大小為dx、沿y方向的大小為dy。M個(gè)棱柱體的厚度(M=mx·my)便是待估算的參數(shù)。棱柱體的厚度pj,j=1,…,M，代表在M個(gè)離散點(diǎn)沉積基底的深度，它與第i個(gè)觀測(cè)點(diǎn)(x=xi,y=yi,z=zi)重力異常gi非線性相關(guān)

計(jì)算第i個(gè)觀測(cè)點(diǎn)()xi,yi,zi的非線性函數(shù)

圖1 解釋模型

重力異常值（灰色輪廓線）是由下伏于沉積包裹體（圖上未給出）的均勻基底構(gòu)建的規(guī)則網(wǎng)格（黑點(diǎn)）插值得到的。包含沉積包裹體的次級(jí)界面是離散成厚度參數(shù)被估計(jì)的三維垂向棱柱體的網(wǎng)格。右邊的插圖展示了第j個(gè)三維棱柱體和第i個(gè)重力異常的垂向分量的坐標(biāo)。

式中，γ是牛頓萬有引力常數(shù);xoj和yoj是第j個(gè)棱柱體中心的x,y坐標(biāo)。

公式（3）即是Chakravarthietal.（2002）介紹的積分方程的封閉形式。

假設(shè)重力數(shù)據(jù)是規(guī)則網(wǎng)格（圖1）插值得到的，它的每個(gè)觀測(cè)點(diǎn)的x,y坐標(biāo)與每個(gè)棱柱體中心的水平坐標(biāo)相一致。定義gi為矢量g→≡() g1,…,gMT的第i個(gè)分量（T代表轉(zhuǎn)置），它包含了由M個(gè)棱柱體模擬的包裹體密度差隨深度按方程（1）衰減所引起的理論重力異常。

2.1 基底深度反演

在方程（6）中Δρ0和α?xí)谙乱还?jié)通過一些步驟的敘述而獲得。一階Tikhonov正則化函數(shù)式（4）TikhonovandArsenin,1977）利用平滑來約束基底起伏。該函數(shù)式中，‖‖·是歐幾里得范數(shù)，R→是一個(gè)矩陣代表一階離散的微分算子（Twomey,1963;Vogel, 2002）。

表達(dá)式（5）中的函數(shù)式介紹了水平坐標(biāo)上使已知深度和估計(jì)深度接近的基底深度的鉆孔信息。該函數(shù)式中，B是貫穿基底起伏面深度zBk,k=1,…,B的鉆孔個(gè)數(shù)，wBkj是B×M階矩陣WB的第kj個(gè)元素，該矩陣的行只包含一個(gè)非零元素，相當(dāng)于一個(gè)單位元素。這個(gè)非零元素位于WB的第k行，與矢量p→的元素有關(guān)系，其相應(yīng)的水平坐標(biāo)是最接近第k個(gè)鉆孔的水平坐標(biāo)。在擬合程度較差的函數(shù)式（6）中，是N維矢量的第i個(gè)元素包含了第i個(gè)觀測(cè)點(diǎn)（2）上計(jì)算的重力異常，δ2是重力數(shù)據(jù)中的干擾識(shí)別的均方差期望。

帶約束的反問題由函數(shù)式（4）和（5）的極小值給出，方程（6）給出了滿足的約束條件，可以改寫成非約束函數(shù)式的極小值的最優(yōu)化問題:

其中，μ(δ)是控制擬合程度的函數(shù)式（8）和先驗(yàn)條件（4)和(5）之間權(quán)重衡量的非負(fù)標(biāo)量。函數(shù)式（8）的最小值通過Marquardt（1963）的方法獲得，在每次迭代中加入海森矩陣的高斯-牛頓近似值（Silvaetal., 2001）。最后，非負(fù)性約束（方程7）由一個(gè)同胚變換給出（例如Barbosaetal.,1999）。然而，這些約束也能由Haskell和Hanson（1981）提出的通過非負(fù)性約束最小二乘算法給出（SilvaDiasetal.2007）。

2.2 估算Δρ0和α

為了估算基底起伏（即為了獲得由方程8給出的非約束問題的最小值），我們第一步需要獲得表面的密度差Δρ*0和密度差隨深度拋物線衰減的因子α*（方程1）的估計(jì)值。在下文中會(huì)得到一對(duì)的值。給定這對(duì)Δρ0和α的值，并且獲得一個(gè)估計(jì)矢量參數(shù)，使?jié)M足重力觀測(cè)值在測(cè)量誤差（6）和非負(fù)性約束（7）以內(nèi)的一階Tikhonov穩(wěn)定函數(shù)式（4）最小化。估算矢量∧p(Δρ0,α)之后，來求方程的值:

接下來，把Θ()Δρ0，α畫到平面Δρ0×α上。重復(fù)這個(gè)過程，該過程是由解釋者設(shè)定的對(duì)不同點(diǎn)的一對(duì)() Δρ0，α產(chǎn)生一對(duì)增量為Δρ0和Δα的離散圖形Θ()Δρ0，α的最小() Δρ0，α。該離散圖形允許一個(gè)Θ()0,α*的視覺估計(jì)。然后，檢查下面的不等式是否成立：Δρ*0,α*是觀測(cè)值滿足已知L的極限的最優(yōu)值；否則，離散區(qū)域Δρ0×α可能不能得到真實(shí)的一對(duì)()

如果滿足上式，那么估算值p∧() Δρ*0,α*，既然這樣，那么離散化的網(wǎng)格是被精細(xì)化的且上述過程是被反復(fù)的。Δρ*

3 組合模型的應(yīng)用

圖2a給出了由模擬沉積盆地的復(fù)雜基底起伏（圖2b）產(chǎn)生的噪聲-干擾布格異常（藍(lán)色實(shí)線）。在有26×78個(gè)節(jié)點(diǎn)的x和y方向（分別代表北南和東西向）網(wǎng)格間距為1km的網(wǎng)格上，重力數(shù)據(jù)是在其z= 0km平面上計(jì)算出的。用標(biāo)準(zhǔn)偏差為0.1mGal的均值為零的高斯偽噪聲干擾理論異常。

圖2b給出了等高線圖和基底起伏的實(shí)際深度的透視圖。假設(shè)一個(gè)非均勻沉積包裹體覆蓋在均勻的基底上，該基底有一個(gè)復(fù)雜的結(jié)構(gòu)骨架被由北西走向的正斷層將基底起伏面分割成復(fù)雜的鑲嵌構(gòu)造凹陷和凸起的兩個(gè)區(qū)域強(qiáng)有力地控制著（如圖2a和b，Ⅰ和Ⅱ）。假設(shè)密度差隨深度衰減符合拋物線定律（方程1）。我們定義區(qū)域Ⅰ和Ⅱ密度變化隨深度的不同拋物線由每個(gè)區(qū)域的不同對(duì)點(diǎn)()Δρ0，α賦值（區(qū)域Ⅰ為-0.6g/cm3和0.10g/cm3/km，區(qū)域Ⅱ?yàn)?0.4g/cm3和0.05g/cm3/km）。區(qū)域Ⅱ是西北狹長(zhǎng)的次級(jí)盆地包含兩個(gè)獨(dú)立的、最大深度達(dá)到7.2km的斷層邊界構(gòu)造凹陷。區(qū)域Ⅰ中，有一系列鄰近的交替出現(xiàn)的構(gòu)造凹陷和構(gòu)造凸起被西北走向的斷層控制。這一區(qū)域的深度范圍為3.5到7.2km。圖2a最突出的特點(diǎn)是，無論是區(qū)域Ⅱ存在的兩個(gè)構(gòu)造凹陷還是區(qū)域Ⅰ存在的四個(gè)構(gòu)造凹陷，都不能很容易的從重力異常檢查區(qū)中推斷出。

3.1 估算Δρ0和α

為了估算三維基底起伏，需要通過圖2中每個(gè)區(qū)域Ⅰ和Ⅱ離散的Θ() Δρ0，α的圖形（9）知道最合適的一對(duì)() Δρ0，α的估計(jì)值。我們建立一個(gè)初始的解釋模型，它由三維垂向分布的棱柱體在x、y方向各自有相同的1.0km的網(wǎng)格間距的26×78的網(wǎng)格組成。然后，我們定義調(diào)查間隔和在區(qū)域Ⅰ和Ⅱ中Δρ0與α的增量從而產(chǎn)生函數(shù)式Θ()Δρ0，α（9）的離散圖形，見圖2a；對(duì)于區(qū)域Ⅰ，我們用到了三個(gè)鉆孔（圖2a中的黑色星號(hào)）的基底面深度的信息。

我們注意到一個(gè)寬的最小值的區(qū)域（圖2a中的黑灰色區(qū)域）包含實(shí)際的一對(duì)() Δρ0，α的值（圖2a中的白色十字）。所有位于最小值區(qū)域里成對(duì)的() Δρ0，α的值與實(shí)際基底起伏好的估算值有關(guān)系，就像插圖中四個(gè)基底起伏（圖2b中的灰色虛線）的估計(jì)，分別是由在1～4點(diǎn)（圖2a中的白點(diǎn)和白色十字）上用成對(duì)的() Δρ0，α的值得到的。這些成對(duì)的值服從可接受的異常擬合（平均誤差≈0.08mGal）。

這些結(jié)果表明函數(shù)式Θ()Δρ0，α（9）的最小值區(qū)域可能不包括實(shí)際的()Δρ0，α的值。它依賴于鉆孔的數(shù)量和分布以及對(duì)重力的響應(yīng)。

3.2 估算三維基底起伏

圖2c給出了估算基底起伏深度的等高線圖和透視圖，它用Δρ0和α的估計(jì)值估計(jì)基底起伏的深度，Δρ0和α是由區(qū)域Ⅰ和Ⅱ履行的一個(gè)系統(tǒng)搜索函數(shù)式Θ() Δρ0，α得到的，并檢查是否滿足不等式(10)，它定義了一個(gè)可能的不確定區(qū)域。在Θ() Δρ0，α的最小區(qū)域中運(yùn)用成對(duì)（Δρ*0=-0.6g/cm3,α*=0.10g/cm3/km）值，它符合數(shù)據(jù)約0.08mGal的平均誤差。將基底深度估計(jì)和實(shí)際深度對(duì)比，我們證明我們的方法在恢復(fù)復(fù)雜基底起伏時(shí)有很好的效果。

該方法可以改進(jìn)呈現(xiàn)不連續(xù)起伏的基底的沉積盆地的基底估計(jì)，例如Barbosaetal（1999）提出的方法。該方法也可以在比沉積盆地更小規(guī)模的類似盆地特征中應(yīng)用，例如（Silvaetal.,2009）廢棄的填埋區(qū)。

圖2 組合測(cè)試

4 結(jié)論

眾所周知，僅僅用重力數(shù)據(jù)確定密度和場(chǎng)源值是不可能的，因此，利用很少點(diǎn)上的基底深度信息去同時(shí)估計(jì)三維基底起伏（場(chǎng)源值）和沉積盆地的密度差隨深度拋物線衰減的參數(shù)。此外，在一些點(diǎn)上用場(chǎng)源體的信息去克服涉及場(chǎng)源體物理性質(zhì)的不明確性。組合模型實(shí)例表明，即使在復(fù)雜的地質(zhì)條件下，該方法也得到了良好的沉積基底起伏估計(jì)。

[1]CristianoM.Martins1,ValeriaC.F.Barbosa,andJo o B.C.Silva,Simultaneous3Ddepth-to-basementanddensitycontrastestimatesusinggravitydataanddepthcontrolatfew points.GEOPHYSICS,VOL.75,NO.3,2010;P121-128.

[2]Barbosa,V.C.F.,P.T.L.Menezes,andJ.B.C.Silva, 2007,Gravitydataasatoolfordetectingfaults:In-depthenhancementofsubtleAlmada’sbasementfaults,Brazil:Geophysics,72,no.3,B59-B68.

[3]Barbosa,V.C.F.,J.B.C.Silva,andW.E.Medeiros, 1997,Gravityinversionofbasementreliefusingapproximate equalityconstraintsondepths:Geophysics,62,1745-1757.

[4]Beltrao,J.F.,J.B.C.Silva,andJ.C.Costa,1991,Robustpolynomialfittingforregionalgravityestimation:Geophysics, 56,80-89.

[5]Chakravarthi,V.,H.M.Raghuram,andS.B.Singh,2002, 3Dforwardgravitymodelingofdensityinterfacesabovewhichthe densitycontrastvariescontinuouslywithdepth:Computersand Geosciences,28,53-57.

[6]Gallardo-Delgado,L.A.,M.A.Perez-Flores,andE.Gomez-Trevino,2003,Aversatilealgorithmforjoint3Dinversionof gravityandmagneticdata:Geophysics,68,949-959.

[7]Garcia-Abdeslem,J.,2005,Thegravitationalattraction ofarightrectangularprismwithdensityvaryingwithdepthfollowingacubicpolynomial:Geophysics,70,no.6,J39-J42.

[8]Leao,J.W.D.,P.T.L.Menezes,J.F.Beltrao,andJ.B.C. Silva,1996,Gravityinversionofbasementreliefconstrainedby theknowledgeofdepthatisolatedpoints:Geophysics,61,1702-1714.

收稿：2014-12-24

10.16206/j.cnki.65-1136/tg.2015.04.007