劉燕東，紀(jì)曉琳，孟小紅，王萬銀，王俊

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

0 引言

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

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

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

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

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

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

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

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

解的不穩(wěn)定性定量度量ρ(μ)是在正則化參數(shù)取不同值時對于多組數(shù)據(jù)反演得到的多組解之間的向量差的無窮范數(shù)中的最大值.由于任意兩個獨立序列偶然產(chǎn)生兩個與彼此非常接近的解的概率很小，故如果一組數(shù)據(jù)產(chǎn)生的解與其他數(shù)據(jù)的解相差太遠(yuǎn)，足以證明反演過程由于解的不穩(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ù)，表達(dá)式為

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

(1)

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

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

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

(2)

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

(3)

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

(4)

將其作為解離差的度量.

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

用于反演的J組數(shù)據(jù)本應(yīng)來自J次野外測量，但是，一般情況下的實際重力勘探只進(jìn)行一次野外測量.通過對一次野外測量所獲得的數(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)造隨機(jī)變量α(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ù)量，所以可以預(yù)測ρ(μ)是單調(diào)遞減的.也就是說，正則化參數(shù)μ越大，不穩(wěn)定性越小.那么，當(dāng)一個單位正則化數(shù)量時Δρk很小，可以認(rèn)為μk是使解穩(wěn)定的正則化參數(shù)：

(10)

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

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

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

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

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

M個二維并置矩形柱體在觀測點處產(chǎn)生的重力異?？梢越票硎境练e盆地基底起伏在觀測點處引起的重力異常，即：

(11)

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

(12)

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

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

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

(13)

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

3.1 目標(biāo)函數(shù)

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

(14)

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

(15)

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

3.2 高斯牛頓法求解

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

(16)

(17)

將目標(biāo)函數(shù)ψ(p)在pk附近進(jìn)行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)定的方法應(yīng)用到模擬二維沉積盆地重力數(shù)據(jù)的盆地基底正則化反演中，測試該方法在通過對一次野外測量所獲得的數(shù)據(jù)分別加入J組不同噪聲得到的J組數(shù)據(jù)與從J次野外測量中得到J組數(shù)據(jù)這兩種情況中的反演效果.

4.1 模型設(shè)置

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

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

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

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

4.2 一次野外測量

為了模擬實測數(shù)據(jù)的情況，對于已用單個單高斯偽隨機(jī)噪聲序列ε1(ri,t)(i=1,…,100，均值為零，標(biāo)準(zhǔn)差為0.1 mGal)擾動圖4a的重力理論觀測值得到含噪聲的模擬實測重力觀測值如圖4b所示.為了說明變量α(ri,t)的噪聲分布類型可能不同于ε1(ri,t)的分布，用25組不同的均勻隨機(jī)噪聲α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組高斯噪聲(均值為零，標(biāo)準(zhǔn)差為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.故本文的方法適用于一般情況下只進(jìn)行一次野外測量的實際重力勘探情況.

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

圖7 25組擾動重力解

5 實際資料測試

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

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

本次研究區(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ù)準(zhǔn)備

位于北加蓬次盆SW-NE走向的測線AB的布格重力異常如圖10a所示，其包含了鹽下、鹽層及鹽上3套沉積層序與基底產(chǎn)生的重力異常，故采用最小曲率法(紀(jì)曉琳等, 2015)對布格重力異常進(jìn)行位場分離后得到反映盆地基底的剩余布格重力異常(圖10b)作為基底起伏反演的基礎(chǔ)數(shù)據(jù).由于北加蓬次盆發(fā)育多套沉積層序且含油氣，且本次測試的主要目的是用實測數(shù)據(jù)證明穩(wěn)定正則化參數(shù)估計方法的可靠性與實用性，故將盆地的構(gòu)造簡化進(jìn)行反演.根據(jù)查得的密度資料(紀(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ù)紀(jì)曉琳等, 2019)

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

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

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

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

利用5.2節(jié)中為正則化反演所準(zhǔn)備的數(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é)果與測線對應(yīng)的地質(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ù)一致，且獲得的反演解非常接近，都能夠反演得到較為準(zhǔn)確的模型基底深度.故該方法適用于一般情況下只進(jìn)行一次野外測量的實際重力勘探.

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

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

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