郭一豪 陳 曉* 楊海燕 張志勇 曾志文 周 勇

（①東華理工大學核技術(shù)應用教育部工程研究中心，江西南昌330013；②東華理工大學地球物理與測控技術(shù)學院，江西南昌330013）

0 引言

Tikhonov正則化反演理論［1］適用于不適定的地球物理反問題。正則化因子的選取直接影響反演結(jié)果，故此，如何正確選取正則化因子一直是地球物理正則化反演領(lǐng)域的研究熱點和難點。

利用一定的準則確定近似最佳正則化因子，且在整個反演過程中取值保持不變的方法，稱為經(jīng)驗定值法。傳統(tǒng)的定值方法包括：無偏風險估計法［2］，廣義交叉驗 證法［3-4］和L曲 線 法［5-6］。Lima等［7］考慮了噪聲的隨機性，提出解的不穩(wěn)定性評價方式，結(jié)合解的穩(wěn)定性確定正則化因子，并分別利用地震、重力和大地電磁測深（MT）數(shù)據(jù)反演試驗驗證了方法的有效性。由于該類方法操作復雜，需要開展大量的計算以確定近似最佳的正則化因子，故此，自適應正則化算法是正則化反演的發(fā)展趨勢。

所謂自適應正則化算法，是指在反演過程中根據(jù)一定的準則自適應地調(diào)整正則化因子的取值的方法。Zhdanov［8］提出了一種根據(jù)初始數(shù)據(jù)擬合函數(shù)與模型約束函數(shù)的比值確定正則化因子的初值，在數(shù)據(jù)擬合效率下降時對正則化因子進行衰減的自適應方法。陳小斌等［9］以MT反演為例，提出基于最平緩約束的自適應算法，該算法強調(diào)數(shù)據(jù)擬合項與模型穩(wěn)定項的平衡，反演結(jié)果穩(wěn)定性高。向陽等［10］以Zhdanov自適應算法為基本框架，提出了改進的自適應算法，該算法對Zhdanov算法所確定的正則化因子的初值做放大處理，并提出了衰減因子與數(shù)據(jù)擬合效率之間的關(guān)聯(lián)方案，基于共軛梯度算法的MT數(shù)據(jù)光滑反演試驗表明，該算法在一定程度上可以降低初始模型對反演結(jié)果的影響，提高解的穩(wěn)定性。黃賢陽等［11］將向陽等［10］提出的自適應算法應用于基于光滑和聚焦約束的MT極值邊界正則化反演，驗證了該算法具有適用于其他穩(wěn)定器的潛質(zhì)。

需要指出是，以往關(guān)于正則化因子的文獻大多是基于線性優(yōu)化算法的。與線性優(yōu)化算法不同，概率思想是非線性優(yōu)化算法的核心，如模擬退火算法［12-13］按照概率接受“高能量”解；遺傳算法［14］和差分進化算法［15］按照概率進行交叉和變異操作；粒子群算法［16］按照概率更新粒子的速度和位置。目前，鮮見文獻在確定正則化因子時，對非線性優(yōu)化算法中的隨機性進行針對性的考量。

綜上所述，從提高反演的穩(wěn)定性方面考量而確定正則化因子取值是本文自適應正則化算法的思路。此外，以往的自適應算法對非線性算法中的隨機性考慮不足，基于此，本文提出一種階段式自適應正則化算法，并分別在MT線性和非線性反演試驗中檢驗其適用性。

1 MT正則化反演

1.1 正則化目標函數(shù)

MT數(shù)據(jù)的正則化反演目標函數(shù)為

式中：α為正則化因子；m為待解模型；φ為數(shù)據(jù)擬合函數(shù)；S為模型約束函數(shù)；A為正演算子；dobs為觀測數(shù)據(jù)；W為模型加權(quán)矩陣；We為約束函數(shù)矩陣；mapr為先驗信息。本文采用有限元法實現(xiàn)MT數(shù)據(jù)正演。

1.2 階段式自適應算法

1.2.1 Zhdanov自適應方案

Zhdanov自適應方案是一種無需復雜計算的自適應正則化算法，其核心步驟包括計算正則化因子初值和自適應衰減兩步［8，17-18］。

正則化因子初值為

式中m0為初始模型。

若反演過程中數(shù)據(jù)擬合函數(shù)上升，或效率低于給定閾值，按照下式對正則化因子進行衰減

式中：n為迭代次數(shù)；q為衰減系數(shù)，取值參考區(qū)間為［0.5，0.9］。

1.2.2 階段式自適應方案

階段式自適應正則化算法的核心思想是，將反演的迭代過程分為若干“階段”，利用“階段”的平均擬合效率代替某次迭代的擬合效率，判斷正則化因子的衰減，旨在降低正則化因子的衰減速度，降低非線性優(yōu)化算法中隨機性的影響，進而提高反演的穩(wěn)定性。該算法基本流程如下。

（1）將迭代反演過程設(shè)置為若干“階段”。

（2）將初始“階段”模型m的平均數(shù)據(jù)擬合函數(shù)與平均模型約束函數(shù)之比值作為正則化因子的初值

式中上劃線表示取均值。對于差分進化（DE）算法，則利用初始種群數(shù)據(jù)擬合函數(shù)的平均值和模型約束函數(shù)的平均值之比值作為正則化因子的初值。

（3）利用當前“階段”的平均數(shù)據(jù)擬合效率自適應調(diào)整下一階段的正則化因子，即計算反演過程中相鄰“階段”的平均數(shù)據(jù)擬合下降效率，并按照下式確定下一“階段”的正則化因子

式中：k＝1，2，3，…，Kmax為“階段”數(shù)；表示“階段”內(nèi)的平均數(shù)據(jù)擬合函數(shù)的擬合效率；ε表示給定的閾值。

（4）判斷是否達到終止條件，若是，則結(jié)束迭代；若否，則返回步驟（3）。

1.3 目標函數(shù)優(yōu)化

為了驗證本文提出的階段式自適應算法在線性和非線性優(yōu)化算法中的適用性，將其應用于二維MT共軛梯度反演和一維MT差分進化反演。

1.3.1 共軛梯度反演

本文采用成熟通用的共軛梯度算法［10］實現(xiàn)二維MT反演，基本流程如下：

式中：Rn為第n次數(shù)據(jù)擬合殘差；mn為待解模型；為梯度方向；為共軛梯度方向；J為靈敏度矩陣；為系數(shù)；為迭代步長。

1.3.2 差分進化算法

DE（差分進化）算法是Storn等［19］基于遺傳算法等進化思想提出的一種非線性全局優(yōu)化算法，其本質(zhì)是多目標函數(shù)優(yōu)化算法，用于求解多維空間中的整體最優(yōu)解。差分進化算法與遺傳算法具有一定的相似性，都是利用隨機的方式生成種群，主要都包括變異、交叉和選擇三個步驟。由于DE算法具有受控參數(shù)少、原理簡單、魯棒性強的特點，受到越來越多的關(guān)注。但與成熟的非線性算法（如模擬退火、遺傳算法）相比，該算法在地球物理反演領(lǐng)域的應用偏少。

本文采用全局尋優(yōu)的DE算法［19-21］實現(xiàn)一維MT反演，該算法的核心環(huán)節(jié)如下。

（1）變異操作

在當前種群中隨機選擇兩個不同的個體，將它們的差向量縮放后與另外的個體進行向量運算，生成經(jīng)過變異操作的中間種群為

式中：i＝1，2，3，…，NP，NP為種群規(guī)模數(shù)；r1、r2和r3均為區(qū)間［1，NP］內(nèi)的隨機整數(shù)，且i≠r1≠r2≠r3；it表示進化代數(shù)；為中間種群的第i個個體；F表示變異因子；表示當前種群的第i個個體。

（2）交叉操作

利用下式開展交叉操作

（3）選擇操作

采用貪婪算法，根據(jù)目標函數(shù)的大小選擇進入新種群的個體

式中：f為適應度函數(shù)，即式（1）中的目標函數(shù)；為第it＋1代的第i個交叉?zhèn)€體；為第it代新種群的第i個個體。

2 模型試驗

2.1 模型一

為了驗證階段式自適應算法在非線性算法中的可行性，選擇文獻［9］中的一維電阻率模型（圖1）開展試驗。在1000～0.0001Hz選取40個對數(shù)等間距頻點，采用隨機生成的模型為初始模型。

圖1 一維電阻率模型［9］

設(shè)計了兩套對比試驗方案：Zhdanov自適應正則化反演和階段式自適應正則化反演。兩種方案僅正則化因子確定方式不同，其余參數(shù)均一致，反演最大迭代次數(shù)為1000，衰減因子q＝0.5，采用最小模型約束進行反演。反演結(jié)果見圖2。

需要指出的是，為了直觀地對比不同正則化方案反演結(jié)果隨迭代次數(shù)的變化，本次試驗特意選擇了較大的迭代次數(shù)。從圖2可以看出，即使Zhdanov自適應方案的目標函收斂效率較高，但其反演結(jié)果隨迭代次數(shù)的增加出現(xiàn)了較大的波動，甚至逐漸變差；而階段式自適應方案的目標函數(shù)曲線平穩(wěn)下降，反演結(jié)果穩(wěn)定性高，與真實模型偏離程度始終在可接受范圍內(nèi)。對比兩個方案的正則化因子迭代曲線，可以看出階段式自適應方案的正則化因子衰減速度更慢，可以更充分地發(fā)揮模型約束函數(shù)的作用，反演過程更加穩(wěn)定。

2.2 模型二

選擇文獻［22］的二維電阻率模型（圖3）開展模型試驗。在320～0.00055Hz范圍內(nèi)選取40個對數(shù)等間距頻點，測點數(shù)共101個，用均勻半空間作為初始模型，利用共軛梯度算法進行光滑反演，最大迭代次數(shù)設(shè)為40。針對TE模式進行反演。

圖2 模型一Zhdanov自適應方案和階段式自適應方案反演結(jié)果

圖3 二維電阻率模型［22］

本文對Zhdanov自適應正則化反演與階段式自適應正則化反演進行對比。這兩種方案僅正則化因子衰減方式不同，其余參數(shù)均一致。初始模型為均勻半空間（電阻率為100Ω·m），采用光滑約束進行反演，以相鄰兩次迭代為一個階段。選取衰減系數(shù)q為0.5，反演結(jié)果見圖4。

分析圖4可以看出：q取0.5時，Zhdanov自適應方案反演結(jié)果中高阻隆起位置偏高，且兩側(cè)也存在虛假的低阻條帶；而階段式自適應方案反演結(jié)果中高阻隆起位置與模型更加接近，且反演模型更加光滑。另外，相對于Zhdanov自適應方案，階段式自適應方案的目標函數(shù)曲線下降更平緩均勻，具有更強的反演穩(wěn)定性。

為了研究衰減系數(shù)q對反演結(jié)果的影響，進一步選取q＝0.6和q＝0.7，分別利用Zhdanov自適應方案與階段式自適應方案進行反演，結(jié)果見圖5和圖6。

對比圖4和圖5，可以看出，Zhdanov自適應方案反演結(jié)果中高阻隆起位置依舊偏高，兩側(cè)出現(xiàn)了虛假的低阻條帶，反演模型光滑性較差。而階段式自適應方案反演結(jié)果兩側(cè)沒有虛假的低阻條帶，且高阻隆起位置與設(shè)計模型更加接近。相對于Zhdanov自適應方案，階段式自適應方案的目標函數(shù)曲線下降更加穩(wěn)定平緩。

綜合分析圖4、圖5和圖6，不難發(fā)現(xiàn)，Zhdanov自適應反演的目標函數(shù)收斂效率更高，但階段式自適應反演結(jié)果更加穩(wěn)定，具體表現(xiàn)在：不同衰減因子取值條件下，階段式自適應反演不存在Zhdanov自適應反演中一直存在的條帶狀假異常，而且其反演結(jié)果中高阻隆起位置與真實模型更加接近。這說明階段式自適應算法在一定程度上為解決衰減因子不易確定的難題提供了新思路，也提示地球物理反演不能僅考慮目標函數(shù)的收斂速度。

為了定量評價反演結(jié)果與真實模型的偏離程度，根據(jù)反演結(jié)果解釋了高阻隆起界面（圖4、圖5和圖6a和圖6b中最右圖中的黑線），利用MATLAB進行線性插值獲得界面起伏函數(shù)，并利用

圖4 模型二q＝0.5時Zhdanov自適應方案與階段式自適應方案的反演結(jié)果

圖5 模型二q＝0.6時Zhdanov自適應方案與階段式自適應方案的反演結(jié)果

圖6 模型二q＝0.7時Zhdanov自適應方案與階段式自適應方案的反演結(jié)果

定量地計算了兩種自適應算法在不同衰減因子條件下的反演結(jié)果與真實模型的偏離度。式中：l表示反演區(qū)域橫向長度；Zreal表示真實高阻隆起界面函數(shù)；Zinv表示反演高阻隆起界面函數(shù)。計算結(jié)果見表1。

表1 解釋的高阻隆起界面與真實界面的絕對平均偏差

通過表1可以看出，在相同衰減因子條件下，階段式自適應方案反演結(jié)果的高阻隆起位置與模型偏差更小，結(jié)果更準確。

此次試驗模型中，高阻隆起的還原很困難。為了定量評價反演過程中的穩(wěn)定性，記錄了每一次迭代反演結(jié)果，記錄了高阻隆起出現(xiàn)雛形時的數(shù)據(jù)擬合函數(shù)和模型約束函數(shù)，并按照式（13）計算了數(shù)據(jù)擬合函數(shù)絕對誤差和模型約束函數(shù)絕對誤差，利用式（14）計算了數(shù)據(jù)擬合函數(shù)相對誤差和模型約束函數(shù)相對誤差。

式中：“Δ”表示絕對誤差；“δ”表示相對誤差；k為迭代次數(shù)；下標“前”、“后”和“均”分別代表相關(guān)參數(shù)的當次迭代值與前一次迭代值、當次迭代值與后一次迭代值、以及前一次迭代值與后一次迭代值的絕對誤差或相對誤差。兩種自適應方案反演結(jié)果中高阻隆起出現(xiàn)時的數(shù)據(jù)擬合函數(shù)和模型約束函數(shù)的變化情況參見表2。

表2 Zhdanov自適應方案與階段自適應方案在高阻隆起出現(xiàn)時的數(shù)據(jù)擬合函數(shù)和模型約束函數(shù)的變化情況對比

此外，本文還嘗試定義了模型改變與數(shù)據(jù)改變的匹配率，即模型穩(wěn)定函數(shù)的改變率與數(shù)據(jù)擬合函數(shù)改變率之比值，來衡量反演過程中的不穩(wěn)定性。為了規(guī)避量綱不一致的問題，采用模型穩(wěn)定函數(shù)與數(shù)據(jù)擬合函數(shù)的相對改變率，采用下式計算匹配率

若匹配率較低，說明在同樣的數(shù)據(jù)擬合函數(shù)相對改變率的條件下，模型的改變更小，即不容易出現(xiàn)反演不穩(wěn)定的情況；反之亦然，若相對匹配率較高，說明在同樣的數(shù)據(jù)擬合函數(shù)相對改變率的條件下，模型的改變更大，即更容易出現(xiàn)反演不穩(wěn)定的情況。

根據(jù)表2數(shù)據(jù)利用式（15）計算兩種反演方案的匹配率見表3。

表3 反演結(jié)果出現(xiàn)高阻隆起時兩種方案的匹配率

分析表3，可以看出在不同的衰減因子條件下，本文提出的階段式自適應算法的匹配率均小于Zhdanov自適應算法，說明階段式自適應算法可以降低反演過程的不穩(wěn)定性。

3 實測資料處理

為了驗證階段式自適應算法的實用性，選取下?lián)P子地區(qū)的某條測線進行處理。

下?lián)P子地區(qū)海相地層具有較好的油氣資源成藏條件，但該地區(qū)地質(zhì)條件復雜、構(gòu)造活動性強、埋深大，遭受多期強烈構(gòu)造活動疊加，勘探難度較大。選取的測線長度約30km，MT測點14個，頻點數(shù)為38，頻率范圍為320～0.001Hz。測線附近有一口鉆井，鉆孔深度為4060m，井底地層為下二疊統(tǒng)灰?guī)r，自上而下地層包括第四系、新近系、古近系、白堊系和下二疊統(tǒng)［23］。

根據(jù)文獻［24］，整理出了研究區(qū)電阻率統(tǒng)計信息，如表4所示。

根據(jù)該區(qū)域電阻率統(tǒng)計以及鉆孔信息，可以看出，該區(qū)域的電性整體呈現(xiàn)“高—低—高”分布特征，即第四系—新近系（Q～N）為高阻層，古近系—白堊系泰州組（E-K2t）為低阻層，白堊系浦口組（K2p）及其下伏地層整體呈現(xiàn)高阻。

表4 測區(qū)電阻率信息統(tǒng)計表

使用均勻半空間（10Ω·m）作為電阻率反演的初始模型，利用共軛梯度算法進行最小模型反演，最大迭代次數(shù)為80，衰減系數(shù)設(shè)置為0.8，以相鄰兩次迭代為一個“階段”。選擇TE模式進行反演，反演結(jié)果見圖7，圖中黑線是根據(jù)電阻率反演結(jié)果、電阻率統(tǒng)計信息以及鉆孔信息初步解釋的電性界面。

為了評價反演結(jié)果，搜集到一條與該MT測線平行的地震時間剖面［25］（圖8）。

綜合圖7和圖8可以發(fā)現(xiàn)，電阻率反演結(jié)果所揭示的界面信息與地震剖面基本吻合，驗證了本文電阻率反演結(jié)果的可靠性，也說明本文階段式自適應算法具有實用性。

圖7 階段式自適應反演結(jié)果

圖8 地震時間剖面

4 結(jié)論與展望

與傳統(tǒng)的自適應正則化算法按照“迭代次數(shù)”自適應調(diào)整正則化因子取值不同，階段式自適應正則式算法按照“階段”自適調(diào)整正則化因子。模型試驗和實測資料處理表明：

（1）階段式自適應正則式算法，可更充分地發(fā)揮模型穩(wěn)定函數(shù)的作用，提高解的穩(wěn)定性，同時也在一定程度上為解決衰減因子不易確定的難題提供了新思路；

（2）階段式自適應正則式算法適用于線性和非線性優(yōu)化算法，具有較高的推廣前景；

（3）通過實測資料處理，檢驗了算法的實用性。

需要指出的是，本文僅對一維MT差分進化算法和二維MT共軛梯度算法驗證了階段式自適應算法的可行性，如何將該算法推廣至高維反演，以及合理地設(shè)計“階段”是今后的研究方向。