劉富強

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

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

劉富強

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

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

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

1 引言

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

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

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

2 研究方法

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

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

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

計算第i個觀測點()xi,yi,zi的非線性函數(shù)

圖1 解釋模型

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

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

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

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

2.1 基底深度反演

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

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

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

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

2.2 估算Δρ0和α

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

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

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

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

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

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

3.1 估算Δρ0和α

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

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

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

3.2 估算三維基底起伏

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

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

圖2 組合測試

4 結(jié)論

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

[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