穩(wěn)定正則化參數(shù)估計方法及其在盆地基底重力異常反演中的應用

劉燕東，紀曉琳，孟小紅，王萬銀，王俊

1 中國地質(zhì)大學(北京)地球物理與信息技術學院，北京 100083 2 長安大學地質(zhì)工程與測繪學院, 西安 710054

0 引言

在位場處理和反演過程中，往往需要求解一個不適定問題.正則化方法實際上是不適定問題的正則化，也就是將不適定問題轉(zhuǎn)換為一個適定問題.正則化參數(shù)的選取是否適當，不僅影響著算法的收斂性以及收斂速度，而且影響正則解是否收斂于問題的真實解.因此，正則化參數(shù)的選取對方程組求解結(jié)果至關重要，如何選取最優(yōu)正則化參數(shù)一直以來都是正則化研究的難點和熱點.

自Tikhonov(1963)提出求解不適定問題的正則化理論和方法以來，地球物理工作者將此方法應用到地球物理的各個領域.欒文貴(1988)在國內(nèi)首次將正則化方法應用于地球物理反演.關于正則化參數(shù)的選取，國內(nèi)外的許多學者提出了不同的方法，這些方法都是基于統(tǒng)計與數(shù)學概念.

在數(shù)學方法中，使用最廣泛的是L曲線法(Hansen, 1992).L曲線法通過解的向量模間接引入穩(wěn)定性分析，應用很廣泛(Calvetti et al., 2000; 馬濤等, 2013).但解的范數(shù)與解的穩(wěn)定性之間的定量關系不總是有效的.此外L曲線方法有多種局限性.Vogel(1996)指出L曲線法再一類特定問題中的收斂不足.Hanke(1996)指出，當信噪比在數(shù)據(jù)中趨于0時正則化參數(shù)不收斂于真解；且當精確解非常平滑時，L曲線方法具有局限性.最后，Xu等(2016)指出L曲線除“全局”角點外還存在一個角點.

綜上，正則化參數(shù)的選擇方法主要存在兩個問題：(1)解的不穩(wěn)定性的難以定義；(2)數(shù)據(jù)誤差對正則化參數(shù)的靈敏度非常低.

Lima等(2019)提出的正則化參數(shù)的選擇方法中，通過“交互”來選擇正則化參數(shù)增加了很多工作量且存在很大的主觀性.然而，不穩(wěn)定性很大成分取決于研究人員的主觀意向，所以定量表示解的不穩(wěn)定性且能夠直接利用其來直接選取正則化參數(shù)是重要的任務.

本文基于Tikhonov正則化方法，針對如何定量度量解的不穩(wěn)定性以及如何產(chǎn)生這種定量度量，如何利用解的不穩(wěn)定性定量度量來直接估計使解穩(wěn)定的最優(yōu)正則化參數(shù)展開研究.通過利用不同的隨機噪聲，將該方法應用到二維沉積盆地基底模擬重力數(shù)據(jù)與實測重力數(shù)據(jù)正則化反演中，證明這種方法在重力數(shù)據(jù)應用中的正確性和實用性.

1 穩(wěn)定正則化參數(shù)估計方法

1.1 解的不穩(wěn)定性定量度量

解的不穩(wěn)定性定量度量ρ(μ)是在正則化參數(shù)取不同值時對于多組數(shù)據(jù)反演得到的多組解之間的向量差的無窮范數(shù)中的最大值.由于任意兩個獨立序列偶然產(chǎn)生兩個與彼此非常接近的解的概率很小，故如果一組數(shù)據(jù)產(chǎn)生的解與其他數(shù)據(jù)的解相差太遠，足以證明反演過程由于解的不穩(wěn)定性而不適定.另一方面，5組數(shù)據(jù)也足以使概率值小于0.1%(Lima et al., 2019).在理想情況下，多組數(shù)據(jù)來自多次野外測量.

1.1.1J組數(shù)據(jù)來自J次野外測量

計算給定正則化參數(shù)值μk產(chǎn)生的解的不穩(wěn)定性定量度量其求解步驟如下:

(1)J次野外測量得到的J組數(shù)據(jù)，表達式為

dj=d+rj,j=1,…,J,

(1)

其中d是一組無噪聲的理論觀測值，rj是偽隨機時變噪聲序列.

(3)計算J組估計解之間的差的絕對值，表達式為

i=j+1,…,J,l=1,…,M.

(2)

(4)計算標量ρji作為Δpji的無窮范數(shù)，即：

(3)

(5)當前正則化參數(shù)μk產(chǎn)生的解不穩(wěn)定性定量度量，表達式為

(4)

將其作為解離差的度量.

1.1.2 從一次野外測量得到J組數(shù)據(jù)

用于反演的J組數(shù)據(jù)本應來自J次野外測量，但是，一般情況下的實際重力勘探只進行一次野外測量.通過對一次野外測量所獲得的數(shù)據(jù)分別加入J組不同噪聲也可得到的J組數(shù)據(jù)，下面從方差角度證明通過人為加入噪聲獲得的J組數(shù)據(jù)與J次測量獲得的J組數(shù)據(jù)對于正則化反演是等價的.

一次野外測量中得到的一組數(shù)據(jù)是位置與時間的函數(shù)，由理論值do(ri),i=1,…,N與噪聲ε1(ri,t)組成：

d1(ri,tj)=do(ri)+ε1(ri,tj),j=1，

(5)

構(gòu)造隨機變量α(ri,t)與ε(t)互相獨立且方差相同或相近，那么數(shù)據(jù)加入J組噪聲后的方差為

V(d1(ri,tj)+αJ(ri,tj))=V(do(ri))

+V(ε1(ri,tj))+V(αJ(ri,tj)),

(6)

由于do(ri)與ε1(ri,tj)是實數(shù)，方差為零，且：

V(αJ(ri,tj))≈V(ε(t)),

(7)

故得到：

V(d1(ri,tj)+αJ(ri,tj))≈V(ε(t))，

(8)

因為不穩(wěn)定性的分析與正則化參數(shù)的估計主要取決于數(shù)據(jù)中的噪聲方差，故式(8)表明，多次測量獲得的J組數(shù)據(jù)，與用J組噪聲分別擾動一次測量數(shù)據(jù)得到的J組數(shù)據(jù)是等價的.

1.2 確定使解穩(wěn)定的正則化參數(shù)的最優(yōu)值

根據(jù)解的不穩(wěn)定性定量度量ρ(μ)的定義，可以利用ρ(μ)來確定正則化參數(shù)的最優(yōu)值.定義解的不穩(wěn)定性定量度量的變化量為

Δρk=ρk-1(μk-1)-ρk(μk),k=2,3,….

(9)

從理論上說，因為μ控制正則化的數(shù)量，所以可以預測ρ(μ)是單調(diào)遞減的.也就是說，正則化參數(shù)μ越大，不穩(wěn)定性越小.那么，當一個單位正則化數(shù)量時Δρk很小，可以認為μk是使解穩(wěn)定的正則化參數(shù)：

(10)

在Δρk滿足式(10)時，認為Δρk很小.其中，β為穩(wěn)定系數(shù)，β根據(jù)求解精度來確定.若β取值過大，則解不穩(wěn)定；若β取值過小，則解過穩(wěn)定.

在實際情況中，可能會出現(xiàn)解的不穩(wěn)定性定量度量偏離單調(diào)行為的情況，通常是由于計算不精密或違背重要的前提造成的.此時需要停止計算，檢查并修改問題所在，重新開始.

2 二維沉積盆地重力異常正演

沉積盆地由上界面、沉積層與基底組成.假設沉積層的上界面是一個平面，在直角坐標系中，z軸向下為正.將沉積層剖分為M個二維矩形棱柱，模擬沉積盆地基底起伏.二維矩形柱體的寬度已知且為常數(shù)，厚度為參數(shù)pj，二維矩形柱體的頂面與沉積盆地的上界面重合，底面與沉積盆地的基底重合(如圖1所示).

圖1 二維沉積盆地模型示意圖

M個二維并置矩形柱體在觀測點處產(chǎn)生的重力異常可以近似表示沉積盆地基底起伏在觀測點處引起的重力異常，即：

(11)

其中fi(p)表示第i個測點O(xi,zi)處的重力異常，第j個柱體的厚度為pj，Δg(xi,zi,pj)表示第j個柱體在第i個測點處引起的重力異常，其表達式為(Rao et al., 1994)：

(12)

其中，G為牛頓萬有引力常數(shù)；σ為沉積層與基底之間的密度差，為常數(shù)；(ξ,ζ)為二維矩形柱體內(nèi)微元的坐標，如圖2所示第j個柱體的坐標范圍為ξj1～ξj2，ζj1～ζj2；柱體厚度參數(shù)pj=ζj2-ζj1；柱體水平寬度dx=ξj2-ξj1.

圖2 二維矩形柱體微元示意圖

在模型測試中，假設重力異常的觀測面位于z=0的平面上，且與沉積層上界面重合，則式(12)可以化簡為

(13)

3 二維沉積盆地基底正則化反演

3.1 目標函數(shù)

基于Tikhonov正則化方法來構(gòu)造正則化反演的目標函數(shù)，計算使得目標函數(shù)極小的模型參數(shù)p，目標函數(shù)為

(14)

φ(p)=(p-p0)TLTL(p-p0)，

(15)

其中p0為先驗信息，可以用已知的深度信息或前一次反演的結(jié)果；L為光滑度矩陣.

3.2 高斯牛頓法求解

式(14)構(gòu)造的目標函數(shù)ψ(p)的一階導數(shù)(即梯度)用g(p)表示；目標函數(shù)的二階導數(shù)(即Hessian矩陣)用H(p)表示.假設當前的模型向量為pk，那么梯度矩陣與Hessian矩陣的表達式為

(16)

(17)

將目標函數(shù)ψ(p)在pk附近進行Taylor級數(shù)展開，并且忽略二次以上的高階項，得到：

ψ(p)≈ψ(pk)+[g(pk)]T(p-pk)

(18)

H(pk)·δpk=-g(pk).

(19)

用奇異值分解法求解線性方程(19)得到模型修正量δpk，模型參數(shù)的迭代公式可表示為

pk+1=pk+δpk，k=0,1,…，

(20)

4 模型測試

通過將這種穩(wěn)定的方法應用到模擬二維沉積盆地重力數(shù)據(jù)的盆地基底正則化反演中，測試該方法在通過對一次野外測量所獲得的數(shù)據(jù)分別加入J組不同噪聲得到的J組數(shù)據(jù)與從J次野外測量中得到J組數(shù)據(jù)這兩種情況中的反演效果.

4.1 模型設置

本次模型測試設計的二維模擬沉積盆地的基底起伏如圖3所示，其中沉積物和基底之間的密度差為-0.3 g·cm-3.

圖3 二維模擬沉積盆地的基底起伏

正演二維模擬沉積盆地的基底起伏產(chǎn)生的重力異常時，進行100個點的模擬觀測，觀測點間隔為0.5 km.將二維沉積盆地剖分為100個寬度為0.5 km的二維矩形棱柱，且棱柱中心的x坐標與每一個觀測值的x坐標一致.正演得到理論重力異常值(圖4a)與含噪聲的重力異常(圖4b).給定3個點的水平坐標與已知深度，作為二維沉積盆地基底深度正則化反演的先驗信息(表1).

表1 反演二維模擬沉積盆地基底深度時的先驗信息

4.2 一次野外測量

為了模擬實測數(shù)據(jù)的情況，對于已用單個單高斯偽隨機噪聲序列ε1(ri,t)(i=1,…,100，均值為零，標準差為0.1 mGal)擾動圖4a的重力理論觀測值得到含噪聲的模擬實測重力觀測值如圖4b所示.為了說明變量α(ri,t)的噪聲分布類型可能不同于ε1(ri,t)的分布，用25組不同的均勻隨機噪聲α25(ri,tj)(對于i=1,…,100，平均值為零并且均方差為0.1 mGal)擾動一組模擬實測觀測值(圖4b)，模擬從一組勘探數(shù)據(jù)獲得J組數(shù)據(jù)的情況.

圖4 沉積盆地基底起伏正演重力異常圖

用高斯牛頓法反演這25組數(shù)據(jù).圖5a中實線為25組反演解的不穩(wěn)定性度量ρ隨正則化參數(shù)μ的變化曲線，以左軸為單位；虛線為數(shù)據(jù)均方誤差RMS隨μ的變化曲線，以右軸為單位；圖5b為單位正則化數(shù)量產(chǎn)生的解不穩(wěn)定性定量度量的變化量Δρk/Δμ隨正則化參數(shù)μ的變化曲線.根據(jù)Δρk/Δμ(穩(wěn)定系數(shù)β取0.005)，直接選取到的正則化參數(shù)的最優(yōu)值為μ*=8.而在圖5a數(shù)據(jù)誤差均方根曲線中，正則化參數(shù)μ從0.01變化到15，而均方誤差RMS僅從0.1497變化到0.1626，由此可見數(shù)據(jù)誤差對正則化參數(shù)的靈敏度非常低.圖6為正則化參數(shù)μ=0.01、0.5、2、8時的25組反演解，μ=0.01時反演解極不穩(wěn)定，μ=0.5、2時反演解較為不穩(wěn)定，μ=8時反演解較穩(wěn)定.且μ=8時，均方根誤差為0.1592 mGal，在可接受誤差范圍內(nèi).故正則化參數(shù)的最優(yōu)估計μ*=8是同時產(chǎn)生穩(wěn)定的解和可接受的誤差的最小的正則化參數(shù).

圖5 (a)不穩(wěn)定性曲線和數(shù)據(jù)誤差均方根曲線；(b)Δρk/Δμ曲線

4.3 J次野外測量

用25組高斯噪聲(均值為零，標準差為0.1 mGal)擾動理論重力異常(圖4a)，模擬從J次野外測量獲得J組數(shù)據(jù)，以證明上述模擬實測數(shù)據(jù)情況的合理性.反演這25組數(shù)據(jù)，根據(jù)單位正則化數(shù)量產(chǎn)生的解不穩(wěn)定性定量度量的變化量Δρk/Δμ(穩(wěn)定系數(shù)β取0.005)，直接選取到的正則化參數(shù)的最優(yōu)值為μ*=8.

反演結(jié)果圖7與圖6對比，可說明模擬從一組勘探數(shù)據(jù)獲得J組數(shù)據(jù)所得到的解與模擬的J次野外測量所得的解非常接近.且在該種情況下，利用穩(wěn)定正則化參數(shù)估計方法確定的最優(yōu)正則化參數(shù)同樣是μ*=8.故本文的方法適用于一般情況下只進行一次野外測量的實際重力勘探情況.

圖6 模擬實測數(shù)據(jù)的25組擾動重力解

圖7 25組擾動重力解

5 實際資料測試

5.1 區(qū)域地質(zhì)特征

為了檢驗穩(wěn)定正則化參數(shù)估計方法在二維沉積盆地實測重力數(shù)據(jù)正則化反演中的應用效果，選擇非洲西海岸的加蓬盆地的部分衛(wèi)星測高重力異常數(shù)據(jù)進行了實際資料測試.

本次研究區(qū)位于非洲西海岸加蓬盆地，加蓬盆地是西非海岸盆地群中的一個富油氣盆地(劉深艷等, 2011).加蓬盆地由3個二級構(gòu)造單元組成(圖8)，以恩科米轉(zhuǎn)換斷層、蘭巴雷內(nèi)高地為界分為南加蓬次盆、北加蓬次盆和內(nèi)次盆(趙紅巖等, 2017).收集到的SW-NE走向的剖面CD的剖面位置如圖8所示，地質(zhì)構(gòu)造資料如圖9所示.本次研究北加蓬次盆中的剖面AB.

圖8 加蓬盆地構(gòu)造劃分與剖面位置示意圖(據(jù)趙紅巖等，2017)

北加蓬次盆(North Gabon sub-Basin)沉積組合具有明顯的三分性，發(fā)育了鹽下、鹽層及鹽上三套沉積層序，含油氣層序以鹽上層系為主(郭念發(fā), 2015).故北加蓬次盆具有十分重要的研究意義.北加蓬次盆從西向東構(gòu)造單元劃分為外裂谷盆地、中央隆起帶與內(nèi)裂谷盆地，共三個構(gòu)造單元(圖9)(劉亞雷 等, 2019).

5.2 反演數(shù)據(jù)準備

位于北加蓬次盆SW-NE走向的測線AB的布格重力異常如圖10a所示，其包含了鹽下、鹽層及鹽上3套沉積層序與基底產(chǎn)生的重力異常，故采用最小曲率法(紀曉琳等, 2015)對布格重力異常進行位場分離后得到反映盆地基底的剩余布格重力異常(圖10b)作為基底起伏反演的基礎數(shù)據(jù).由于北加蓬次盆發(fā)育多套沉積層序且含油氣，且本次測試的主要目的是用實測數(shù)據(jù)證明穩(wěn)定正則化參數(shù)估計方法的可靠性與實用性，故將盆地的構(gòu)造簡化進行反演.根據(jù)查得的密度資料(紀曉琳等, 2019)(如表2)，取沉積層密度為3套沉積層的平均密度2.35 g·cm-3，盆地基底密度為2.8 g·cm-3，所以本次實際資料處理所取的沉積層與基底的剩余密度為-0.45 g·cm-3.根據(jù)北加蓬次盆的剖面CD的地質(zhì)構(gòu)造信息(圖9)，給出用于正則化反演的先驗信息(位置及深度)如表3所示.

表2 密度資料(單位：g·cm-3)(據(jù)紀曉琳等, 2019)

圖9 北加蓬次盆CD剖面構(gòu)造劃分(據(jù)劉亞雷等，2019)

表3 正則化反演的深度先驗信息

圖10 (a)測線布格重力異常；(b)剩余布格重力異常

5.3 實際資料反演結(jié)果

利用5.2節(jié)中為正則化反演所準備的數(shù)據(jù)，用穩(wěn)定正則化參數(shù)估計方法在正則化反演中選擇最優(yōu)正則化參數(shù)得到最優(yōu)的反演結(jié)果.解的不穩(wěn)性定量度量隨正則化參數(shù)的變化曲線如圖11a所示，根據(jù)產(chǎn)生的Δρk/Δμ(穩(wěn)定系數(shù)β取0.005)(圖11b)，直接選取的最優(yōu)正則化參數(shù)為μ*=13.最優(yōu)反演結(jié)果如圖12b所示，反演的15組解穩(wěn)定性較好，解不穩(wěn)定性定量度量為0.1786.

圖11 (a)解不穩(wěn)定定量度量隨正則化參數(shù)單調(diào)遞減；(b)Δρk/Δμ曲線

盆地基底反演結(jié)果的深度范圍在6.0～7.4 km，反演結(jié)果與測線對應的地質(zhì)構(gòu)造基底深度基本吻合，基底整體起伏特征基本與地質(zhì)剖面中基底特征一致，且反演的15組結(jié)果穩(wěn)定性較好(如圖12所示).說明穩(wěn)定正則化參數(shù)估計方法在實際數(shù)據(jù)中確定的最優(yōu)正則化參數(shù)可以得到穩(wěn)定最優(yōu)的反演解.

圖12 北加蓬次盆AB剖面基底反演結(jié)果與地質(zhì)剖面對比圖

6 結(jié)論

(1)通過對一次野外測量數(shù)據(jù)加入不同噪聲得到的J組數(shù)據(jù)與從J次野外測量中得到的J組數(shù)據(jù)這兩種情況的模型測試結(jié)果表明，在這兩種情況下利用穩(wěn)定正則化參數(shù)估計方法直接確定的最優(yōu)正則化參數(shù)一致，且獲得的反演解非常接近，都能夠反演得到較為準確的模型基底深度.故該方法適用于一般情況下只進行一次野外測量的實際重力勘探.

(2)將穩(wěn)定正則化參數(shù)估計方法應用到北加蓬次盆剩余布格重力異常中進行盆地基底反演，最優(yōu)反演結(jié)果是穩(wěn)定的且符合當?shù)氐刭|(zhì)構(gòu)造認識，基底整體深度與起伏特征基本與地質(zhì)剖面中基底特征一致.故在實測重力數(shù)據(jù)中，穩(wěn)定正則化參數(shù)估計方法也可以直接確定最優(yōu)正則化參數(shù)以得到穩(wěn)定最優(yōu)反演解.

(3)應用穩(wěn)定正則化參數(shù)估計方法所需的前提條件非常少，且與隨機噪聲分布無關，故該方法具有廣泛的適用性.

致謝感謝評審專家提出的修改意見.