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

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

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

0 引言

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

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

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

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

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

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

1 穩定正則化參數估計方法

1.1 解的不穩定性定量度量

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

1.1.1J組數據來自J次野外測量

計算給定正則化參數值μk產生的解的不穩定性定量度量其求解步驟如下:

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

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

(1)

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

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

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

(2)

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

(3)

(5)當前正則化參數μk產生的解不穩定性定量度量，表達式為

(4)

將其作為解離差的度量.

1.1.2 從一次野外測量得到J組數據

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

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

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

(5)

構造隨機變量α(ri,t)與ε(t)互相獨立且方差相同或相近，那么數據加入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)是實數，方差為零，且：

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

(7)

故得到：

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

(8)

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

1.2 確定使解穩定的正則化參數的最優值

根據解的不穩定性定量度量ρ(μ)的定義，可以利用ρ(μ)來確定正則化參數的最優值.定義解的不穩定性定量度量的變化量為

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

(9)

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

(10)

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

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

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

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

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

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

(11)

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

(12)

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

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

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

(13)

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

3.1 目標函數

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

(14)

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

(15)

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

3.2 高斯牛頓法求解

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

(16)

(17)

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

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

(18)

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

(19)

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

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

(20)

4 模型測試

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

4.1 模型設置

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

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

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

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

4.2 一次野外測量

為了模擬實測數據的情況，對于已用單個單高斯偽隨機噪聲序列ε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)，模擬從一組勘探數據獲得J組數據的情況.

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

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

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

4.3 J次野外測量

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

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

圖6 模擬實測數據的25組擾動重力解

圖7 25組擾動重力解

5 實際資料測試

5.1 區域地質特征

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

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

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

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

5.2 反演數據準備

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

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

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

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

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

5.3 實際資料反演結果

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

圖11 (a)解不穩定定量度量隨正則化參數單調遞減；(b)Δρk/Δμ曲線

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

圖12 北加蓬次盆AB剖面基底反演結果與地質剖面對比圖

6 結論

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

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

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

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