劉富強

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

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

劉富強

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

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

重力數據深度控制三維基底深度密度差

1 引言

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

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

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

2 研究方法

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

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

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

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

圖1 解釋模型

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

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

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

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

2.1 基底深度反演

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

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

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

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

2.2 估算Δρ0和α

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

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

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

3 組合模型的應用

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

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

3.1 估算Δρ0和α

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

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

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

3.2 估算三維基底起伏

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

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

圖2 組合測試

4 結論

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

[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