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

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

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

0 引言

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

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

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

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

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

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

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

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

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

1.1.1J組數(shù)據(jù)來(lái)自J次野外測(cè)量

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

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

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

(1)

其中d是一組無(wú)噪聲的理論觀(guān)測(cè)值，rj是偽隨機(jī)時(shí)變?cè)肼曅蛄?

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

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

(2)

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

(3)

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

(4)

將其作為解離差的度量.

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

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

一次野外測(cè)量中得到的一組數(shù)據(jù)是位置與時(shí)間的函數(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)互相獨(dú)立且方差相同或相近，那么數(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í)數(shù)，方差為零，且：

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

(7)

故得到：

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

(8)

因?yàn)椴环€(wěn)定性的分析與正則化參數(shù)的估計(jì)主要取決于數(shù)據(jù)中的噪聲方差，故式(8)表明，多次測(cè)量獲得的J組數(shù)據(jù)，與用J組噪聲分別擾動(dòng)一次測(cè)量數(shù)據(jù)得到的J組數(shù)據(jù)是等價(jià)的.

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

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

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

(9)

從理論上說(shuō)，因?yàn)棣炭刂普齽t化的數(shù)量，所以可以預(yù)測(cè)ρ(μ)是單調(diào)遞減的.也就是說(shuō)，正則化參數(shù)μ越大，不穩(wěn)定性越小.那么，當(dāng)一個(gè)單位正則化數(shù)量時(shí)Δρk很小，可以認(rèn)為μk是使解穩(wěn)定的正則化參數(shù)：

(10)

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

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

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

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

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

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

(11)

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

(12)

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

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

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

(13)

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

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

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

(14)

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

(15)

其中p0為先驗(yàn)信息，可以用已知的深度信息或前一次反演的結(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級(jí)數(shù)展開(kāi)，并且忽略二次以上的高階項(xiàng)，得到：

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

(18)

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

(19)

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

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

(20)

4 模型測(cè)試

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

4.1 模型設(shè)置

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

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

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

表1 反演二維模擬沉積盆地基底深度時(shí)的先驗(yàn)信息

4.2 一次野外測(cè)量

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

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

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

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

4.3 J次野外測(cè)量

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

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

圖6 模擬實(shí)測(cè)數(shù)據(jù)的25組擾動(dòng)重力解

圖7 25組擾動(dòng)重力解

5 實(shí)際資料測(cè)試

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

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

本次研究區(qū)位于非洲西海岸加蓬盆地，加蓬盆地是西非海岸盆地群中的一個(gè)富油氣盆地(劉深艷等, 2011).加蓬盆地由3個(gè)二級(jí)構(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)造單元?jiǎng)澐譃橥饬压扰璧?、中央隆起帶與內(nèi)裂谷盆地，共三個(gè)構(gòu)造單元(圖9)(劉亞雷 等, 2019).

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

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

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

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

表3 正則化反演的深度先驗(yàn)信息

圖10 (a)測(cè)線(xiàn)布格重力異常；(b)剩余布格重力異常

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

利用5.2節(jié)中為正則化反演所準(zhǔn)備的數(shù)據(jù)，用穩(wěn)定正則化參數(shù)估計(jì)方法在正則化反演中選擇最優(yōu)正則化參數(shù)得到最優(yōu)的反演結(jié)果.解的不穩(wěn)性定量度量隨正則化參數(shù)的變化曲線(xiàn)如圖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/Δμ曲線(xiàn)

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

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

6 結(jié)論

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

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

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

致謝感謝評(píng)審專(zhuān)家提出的修改意見(jiàn).