大地電磁測(cè)深階段式自適應(yīng)正則化反演

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

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

0 引言

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

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

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

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

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

1 MT正則化反演

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

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

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

1.2 階段式自適應(yīng)算法

1.2.1 Zhdanov自適應(yīng)方案

Zhdanov自適應(yīng)方案是一種無(wú)需復(fù)雜計(jì)算的自適應(yīng)正則化算法，其核心步驟包括計(jì)算正則化因子初值和自適應(yīng)衰減兩步［8，17-18］。

正則化因子初值為

式中m0為初始模型。

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

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

1.2.2 階段式自適應(yīng)方案

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

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

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

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

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

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

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

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

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

1.3.1 共軛梯度反演

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

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

1.3.2 差分進(jìn)化算法

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

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

（1）變異操作

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

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

（2）交叉操作

利用下式開展交叉操作

（3）選擇操作

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

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

2 模型試驗(yàn)

2.1 模型一

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

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

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

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

2.2 模型二

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

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

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

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

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

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

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

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

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

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

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

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

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

表1 解釋的高阻隆起界面與真實(shí)界面的絕對(duì)平均偏差

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

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

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

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

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

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

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

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

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

3 實(shí)測(cè)資料處理

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

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

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

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

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

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

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

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

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

圖8 地震時(shí)間剖面

4 結(jié)論與展望

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

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

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

（3）通過(guò)實(shí)測(cè)資料處理，檢驗(yàn)了算法的實(shí)用性。

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