一種重力多界面反演方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明涉及地球物理勘探技術(shù)領(lǐng)域，特別是涉及一種重力多界面反演方法，其為 重力勘探中的密度界面反演方法，可以得到地下多層界面的起伏分布，且具有較高的運算 效率。
【背景技術(shù)】
[0002] 地質(zhì)勘探主要應(yīng)用的地球物理勘探方法包括重力勘探、磁法勘探、電法勘探和地 震勘探等。重力勘探通過利用重力儀觀測地下物質(zhì)密度差異引起的重力異常以查明地下的 地質(zhì)構(gòu)造和巖性異常體。相較于地震等勘探方法而言，重力勘探具有橫向分辨率好、成本低 廉等優(yōu)點。重力勘探分為三個大的環(huán)節(jié):野外重磁異常采集、重磁異常處理和重磁異常解 釋。
[0003] 從地質(zhì)角度，重力反演問題的目標(biāo)在于尋找、研究或推斷金屬或非金屬礦體和研 究地質(zhì)構(gòu)造。從地球物理角度，重力反演問題的目標(biāo)在于確定地質(zhì)體的幾何和物性參數(shù)，確 定物性分界面的深度和起伏，以及密度的分布等等。
[0004] 總體來看，當(dāng)前重力反演方法研究主要可分為兩大類。第一類是密度分布研究，著 重于物性（即密度)在地下半空間中的分布情況的研究;第二類是界面反演研究，關(guān)注點在 于明顯的密度分界面，一般用于處理密度已知(或密度變化規(guī)律已知）的情形下密度分界面 的起伏變化情況。
[0005] 重力界面反演方法的研究起步較早，自20世紀(jì)60年代起，Bott(1960)、Cordell et al. (1968)等即對重力異常在空間域的正反演問題進行了初步研究。這個時期的計算方法 都是基于空間域的，通過長方體模型進行擬合，逐次逼近消去剩余值，無法利用先驗信息。 隨著計算機技術(shù)的發(fā)展，快速傅里葉變換技術(shù)(FFT)投入使用，Parker(1973)借助FFT，給出 了頻率域計算地下單層物質(zhì)重力異常的正演公式，較以往的離散模型方法速度提高了至少 一個數(shù)量級(馮銳等，1986)，極大地提高了計算效率。隨后由01denburg(1974)加以發(fā)展，建 立了反演問題迭代求解方法，后人稱為Parker-Oldenburg反演法，目前已經(jīng)成為了地球物 理領(lǐng)域最常用的重力反演方法之一，具有理論嚴(yán)謹(jǐn)、適應(yīng)性強和快速有效的特點。除此以 外，也有其他科研人員提出了一些反演方法，如Spector et al. (1970)的功率譜分析方法、 Chavez et al. (1985)的Parker頻譜展開線性反演法和Pederson(1977)的廣義矩陣反演法 等。這些方法里有一些是對Parker-Oldenburg方法的改進或引申，另外一些方法則常常用 于某些特定情況，比如功率譜方法適用于對物質(zhì)界面的大體預(yù)估。不過總體而言，以上方法 都存在同樣的局限性:假定地下物質(zhì)為單一界面模型，如此通過重力異常的反演只能得到 單一界面的深度變化情況。針對這一問題，王萬銀等（1993)給出了雙界面反演方法，實質(zhì)采 用了單層地殼雙界面模型，方法上與Parker-Oldenburg方法非常類似，需要額外考慮的是 "調(diào)整深度"的合理選取以加快收斂速率。該方法理論較為簡易，也可以引入控制點進行約 束(黃松，2010)，但模型僅針對單層物質(zhì)的頂?shù)捉缑?，?yīng)用范圍有限，對諸如俯沖帶地區(qū)"雙 莫霍"的反演問題無能為力，且只能反演兩個界面，而且密度、頂界面、底界面三者的反演進 行中，對其中一項的反演總是相當(dāng)于假定其他兩個參數(shù)為已知確定值，使用上造成了一定 的局限性。為避開這種將未知界面"假定為已知確定值"的情況，一些基于雙界面模型發(fā)展 的方法(馮旭亮等，2014)干脆以地形作為雙界面模型的頂界面，實質(zhì)又退化為了單界面反 演問題，且因模型本身受限，關(guān)注點轉(zhuǎn)變?yōu)槌练e盆地基底研究，方法無法用于深層界面反 演。其他的一些非線性全局最優(yōu)化方法(朱自強等，1995; Snopek，2005;于鵬等，2007;李倩 等，2010;明圓圓等，2012)采用非線性全局最優(yōu)化方法，包括模擬退火算法、遺傳算法和人 工魚群算法等，這類算法為避免線性/非線性方程組大規(guī)模矩陣求解問題而采用了零階算 法，其優(yōu)點是對初始模型的依賴程度低，但參數(shù)取值范圍和選取參數(shù)存在隨機性，采用啟發(fā) 式算法導(dǎo)致每次的反演結(jié)果變動很大，且求解效率不高。
[0006] Martins et al. (2011)和馮旭亮等(2014)針對沉積盆地基底反演問題，采用過全 變差泛函與控制點信息相結(jié)合的方式進行研究。但需要指出，Martins et al. (2011)和馮 旭亮等（2014)的方法針對的僅是單一界面反演問題，而且反演對象為淺層物質(zhì)，無法對深 層物質(zhì)進行反演。
[0007] 綜上所述，現(xiàn)有的地球物理勘探方法難以處理更為復(fù)雜的多界面反演情況，并且 運算效率不高。

【發(fā)明內(nèi)容】

[0008] 本發(fā)明主要解決的技術(shù)問題是提供一種重力多界面反演方法，可以得到地下多層 界面的起伏分布，且具有較高的運算效率和反演精度。
[0009] 為解決上述技術(shù)問題，本發(fā)明采用的一個技術(shù)方案是：提供一種重力多界面反演 方法，該方法包括以下步驟：
[0010] (1)獲取觀測點的信息、重力觀測值、初始模型范圍、地下介質(zhì)信息(地下介質(zhì)的相 關(guān)數(shù)據(jù)）、待反演層位信息、與模型泛函各部分的正則化權(quán)重；
[0011] (2)對地下介質(zhì)進行網(wǎng)格剖分；
[0012] (3)計算初始模型網(wǎng)格剖分后的各柱狀體對觀測點的重力效應(yīng)并進行累加，得到 重力正演值；
[0013] (4)根據(jù)重力觀測值與重力正演值，計算失配泛函？:^)，其中，X為待求的各柱狀 體底界面深度；
[0014] (5)計算全變差函數(shù)F2(x);
[0015] (6)計算模型界面反演偏差函數(shù)F3(x);
[0016] (7)根據(jù)步驟(5)~(6)的計算結(jié)果構(gòu)成模型泛函，并根據(jù)模型泛函各部分的正則 化權(quán)重λ2~λ 3與失配泛函Fl(X)-起構(gòu)成目標(biāo)函數(shù)：
[0017] F(x) =Fi(x)+A2F2(x)+hF3(x);
[0018] (8)采用迭代方法使目標(biāo)函數(shù)最小，期間根據(jù)迭代的結(jié)果與給定的不等約束，不斷 修改模型以使反演結(jié)果趨向真實值。
[0019] 其中，步驟(3)進一步包括利用長方體重力正演公式計算初始模型網(wǎng)格剖分后的 各柱狀體對觀測點的重力效應(yīng)并進行累加，得到重力正演值。
[0020] 其中，步驟(1)進一步包括：
[0021] 獲取用于約束反演模型的信息，其包括控制點信息、模型各層上下界信息、失配泛 函協(xié)方差矩陣以及模型界面反演偏差矩陣，其中，控制點信息包括控制點位置信息；
[0022] 步驟(4)和步驟(5)之間還包括:步驟(45)存在控制點信息時，計算控制點信息函 數(shù)F 4(x);
[0023] 步驟(7)包括:根據(jù)步驟(45)、（5)以及(6)計算結(jié)果構(gòu)成模型泛函，并根據(jù)模型泛 函各部分的正則化權(quán)重λ2~λ 4與失配泛函FKx)-起構(gòu)成目標(biāo)函數(shù)：
[0024] F(X) =Fl(X)+A2F2(X)+X3F3(X)+X4F4(X)。
[0025] 其中，步驟(8)具體包括:不斷重復(fù)步驟(3)~(8)，直至達到迭代退出條件為止，最 終求得各層界面深度反演結(jié)果，其中，迭代退出條件包括迭代次數(shù)達到預(yù)設(shè)的次數(shù)以及迭 代的結(jié)果誤差小于預(yù)設(shè)的誤差。
[0026] 其中，觀測點的信息包括觀測點數(shù)目和各觀測點位置(x，y，z)，其中，用于表示X值 的X軸和用于表示y值的Y軸位于同一平面，用于表示z值的Z軸向下為正，X軸、Y軸以及Z軸三 者構(gòu)成右手直角坐標(biāo)系；
[0027]重力觀測值包括重力異常觀測值和/或重力梯度觀測值，當(dāng)進行模型實驗研究時， 由標(biāo)準(zhǔn)模型計算出重力異常正演值和/或重力梯度正演值，若使用重力梯度觀測值，則重力 梯度觀測值至少包括梯度張量的其中一個分量；
[0028] 地下介質(zhì)的相關(guān)數(shù)據(jù)包括地下介質(zhì)層數(shù)、各層介質(zhì)的密度和底界面初始深度；
[0029] 待反演層位信息包括要反演的層位編號，層位編號為由淺部到深部排列或自定義 的編號。
[0030] 其中，重力梯度張量的形式為：
[0032]其中G為萬有引力常數(shù)，P為密度，v為體積微元，尸和f分別為觀測點位置和場源 位置，算子▽作用于觀測點，在χ-γ-ζ坐標(biāo)系下的重力梯度為
[0036]其中，（114243) = (1，7,2)，.沒為偏微分符號，&為在坐標(biāo)1上的偏微分。
[0037] 其中，步驟(4)中，失配泛函形式為
[0039] 其中，fi(x)的取值如下：
[0040] 當(dāng)無重力梯度觀測值時：;1；'1(1) = 8瓜￥。1)；^)-8瓜￥?！?1(叉），
[0041 ]當(dāng)有重力梯度觀測值時：
7
[0042 ] 其中，gra Vobs (X)為重力觀測值，gra Veal (X)為重力正演值，gra vgraobs (X)為重力梯 度觀測值，gravgra^U)為重力梯度正演值，Cddg為失配泛函協(xié)方差矩陣，Μ為觀測點數(shù)目。 [0043]其中，步驟(45)中，控制點信息函數(shù)為
[0045] 其中
[0046] f4(x)=k_Wx
[0047] 其中，W是一個BXN的稀疏矩陣，B為控制點數(shù)目，N為待求柱狀體個數(shù)，稀疏矩陣?yán)?面元素非0即1，保證在相應(yīng)控制點i上的柱狀體底界面深度^趨近于k里的對應(yīng)元素值，而n 4 為歸一化