基于混合快速共軛梯度法的有限差分對比源反演

王豆豆 王守東 鄒少峰 高艷霞 劉 晗

(①中國石化石油物探技術研究院，江蘇南京 211103； ②中國石油大學(北京)油氣資源與探測國家重點實驗室，北京 102249)

0 引言

全波形反演是一種精確的地下彈性參數重建方法，它充分利用地震記錄中包含的動力學和運動學信息，以波動方程為數學模型建立目標泛函，通過不斷迭代優化目標泛函使正演波場與實際波場差值最小，得到目標泛函的最優解，重構地下物性參數，服務油氣勘探開發。全波形反演既可以在時間域進行也可以在頻率域進行[1-3]。全波形反演主要包括正演和反演兩個步驟。正演過程中常用的方法有迭代法和直接求接法；反演過程利用迭代優化方法逐步對初始模型進行修改直到滿足誤差要求。迭代優化是全波形反演中最重要的部分。全波形反演的優化算法主要包括梯度法和牛頓法。梯度法主要包括最速下降法[4-6]和非線性共軛梯度法[7-8]；牛頓法主要包括高斯牛頓法[9]、擬牛頓法[10]、截斷牛頓法[11]和全牛頓法[9,12]。牛頓法收斂快但計算量和內存消耗較大，在實際反演過程中主要采用非線性共軛梯度法和高斯牛頓法。在每次迭代優化中都需要多次正演計算，這會大大增加計算量，因此研究一種高效的反演算法尤為重要。

為解決擾動較大的逆散射問題，同時提高解逆散射問題的效率，van den Berg等[13]從逆散射理論出發，基于源型積分方程[14]提出積分型對比源反演(IE-CSI)方法。該方法是一種解逆散射問題的高效方法，通過不斷優化目標函數更新對比源(散射源)和對比度(擾動)，在迭代優化的過程中只需進行一次全正演運算，從而提高了反演效率。van Dongen等[15]將該方法應用到三維聲波成像，驗證了該方法的高效性。Wang等[16]采用WKBJ近似方法求解Green函數的思路研究了對比源反演方法。由于積分型對比源反演是基于積分方程提出的，在實際反演中需要求解Green函數，在均勻或者水平漸變的背景模型條件下可以得到Green函數，但當模型比較復雜時就很難求解Green函數。為了進一步拓展對比源反演方法的適用范圍，Abubakar等[17-18]研究了有限差分對比源反演(FDCSI)方法，并且將該方法成功推廣到地震全波形反演。FDCSI適用于較復雜的背景模型，其構造的正演算子只與背景模型和頻率有關。由于迭代求解過程中頻率和背景模型保持不變，因此不需要重新構造正演算子，只需進行一次LU分解，提高了計算效率。基于以上優點，FDCSI具有處理大尺度數據的優勢[19]。Han等[20]和He等[21]分別對地震彈性波和聲波數據進行FDCSI反演，實現了多參數反演。段曉亮等[22]利用對比源反演研究了地層衰減對地震波速度逆散射反演的影響。

FDCSI計算效率較高，但是對目標函數的優化采用的是常規共軛梯度法，因而收斂較慢。為了解決收斂速度的問題，進一步提高FDCSI效率，本文在反演中采用混合快速共軛梯度(HFCG)法優化目標函數。HFCG法是在快速迭代收縮閾值算法(FISTA)[23-25]和快速共軛梯度法(FCG)[26]的基礎上改進得到的一種適用于FDCSI的快速收斂優化方法。FISTA是一種一階優化近似分裂算法，該方法能夠加速目標泛函收斂[25]。HFCG比常規共軛梯度法收斂更快。該方法通過修改迭代公式加速目標泛函收斂，因此不會增加總的計算量，能夠很好地平衡計算量與收斂速度。模型測試結果表明，該方法能夠明顯加速收斂，并提高反演結果精度。

1 FDCSI基本理論

1.1 問題描述

對比源反演是為了解決逆散射問題[12](圖1)。當目標域中的物性參數(如速度、密度等)不同于總區域背景模型參數時，對比度χ≠0(存在擾動)。地震波傳播到目標域時會產生散射波，并被布置在S域上的檢波器接收到，利用接收到的散射波反演對比度χ、重建目標域D的物性參數即是解逆散射問題。

式中：H和Hb分別為與模型和背景有關的稀疏阻抗矩陣；Sj為震源項。

圖1 逆散射示意圖

定義對比度χ(模型擾動)為[18]

(3)

式中：r為空間坐標；k和c分別為模型的波數和速度；kb和cb分別為背景模型的波數和速度。

根據總波場、背景波場和散射波場的關系，可得[16]

(4)

(5)

在S域上接收到的散射波為[17]

(6)

式中MS為將T域中的散射波映射到S域的算子。式(6)即為數據方程。

根據式(4)可以得到狀態方程

(7)

MD為將T域中的散射波映射到D域的算子。

引入對比源(散射源)[18]

(8)

式(6)和式(7)可以改寫為[18]

(9)

(10)

在地震勘探中，由于目標空間等同于地下空間，因此實際情況中T域與D域是相等的。

1.2 有限差分對比源反演

在頻率域FDCSI過程中，需要同時滿足數據方程式(9)和狀態方程式(10)，不斷優化求解對比源Wj和對比度χ，目標泛函為[17-18]

C(χ,Wj)=CS(Wj)+CD(Wj,χ)

(11)

定義S域和D域的L2范數為[17]

(12)

從式(11)可以看出，反演過程中需要不斷更新對比源Wj和對比度χ。采取先更新對比源Wj、再更新對比度χ的方式完成每一次迭代。求解目標泛函關于對比源Wj的Fréchet導數，得到第n次迭代目標泛函關于對比源Wj的梯度[17]

(13)

式中

(14)

采用常規PR(Polak-Ribiere)共軛梯度法[30]更新對比源，其迭代公式為

(15)

(16)

(17)

式中Re(·)表示求實部。

初始模型會影響反演結果，因此采取波場反傳方法得到初始對比源Wj,0和初始對比度χ0[14]

(18)

(19)

(20)

FDCSI算法能夠高效求解逆散射問題，但是常規共軛梯度法的優化方法收斂速度慢，降低了反演效率。因此本文將快速迭代收縮閾值算法(FISTA)[23]引入FDCSI，提高反演效率。

2 混合快速共軛梯度法

HFCG法與常規共軛梯度法相比，主要是修改了對比源更新迭代公式，這僅增加少量點乘運算，但大大提高了反演的收斂速度。常規共軛梯度法的迭代公式為式(15)，式中的vj,n即是通過PR共軛梯度法得到的搜索方向。HFCG法參考快速迭代收縮閾值算法，通過修改式(15)中的Wj,n-1達到快速收斂的效果，用Qj,n-1替代式(15)中的Wj,n-1。Qj,n-1的表達式為[23]

(21)

(22)

式中t1=1，修改后的迭代形式如下[22]

(23)

當目標泛函殘差趨于收斂時，式(23)的迭代收斂曲線存在抖動，因此在趨于收斂時采用常規共軛梯度法進行優化迭代，進一步減小目標泛函殘差。則HFCG的迭代公式為

(24)

式中Nunstable為迭代收斂曲線出現不穩定現象時的迭代次數。

3 模型算例

為了驗證方法的有效性，分別針對水平層狀模型和重采樣后的Marmousi模型進行數值模擬實驗。實驗中把上一個頻率反演的結果作為下一個頻率反演的背景模型。

3.1 水平層狀模型

水平層狀模型如圖2a所示，該模型共有七層，速度范圍是1800～2600m/s，模型深度為1225m。模型被離散為100(x)×49(z)個網格。采用地表放炮，共13炮，炮間距為200m；101個檢波器均勻分布在地表，道間距為25m，每炮激發共101道接收。根據頻率選取原則[31]選擇6個頻率進行反演，分別為3.0、4.2、6.1、8.7、12.5、17.8Hz，每個頻率迭代30次，把低頻點反演結果作為高頻點反演的背景模型。背景模型是真實模型經二維平滑處理得到的，如圖2b所示。

圖2 水平層狀速度模型

分別采用HFCG法和常規共軛梯度法得到的反演結果見圖3。對比圖3a與圖3b可以看出，采用HFCG法得到的結果分辨率更高、收斂效果更好。兩種方法反演結果與真實速度模型的L2范數殘差分別為4292.3和5376.5，可見HFCG法反演結果更接近真實速度。圖4為x=1000m和x=1600m處縱向速度對比，可以看出HFCG法對界面具有更好的識別和刻畫能力。

圖3 采用HFCG(a)和常規共軛梯度法(b)反演的速度剖面

圖5為兩種方法6個頻率的目標泛函收斂曲線對比。從圖中可以看出，HFCG法具有非常明顯的快速收斂優勢，該方法在很少的迭代次數時就可實現目標泛函的收斂。

3.2 Marmousi模型

對Marmousi模型進行抽稀，抽稀后的速度模型如圖6a所示。模型被離散為368×160個網格。采用地表激發和地表接收的觀測系統，共47炮，炮間距為200m，369道檢波器均勻分布在地表，道間距為25m，每炮激發369道接收。選取9個頻率進行反演，分別為2、3、4、5、6、8、10、12、16Hz，每個頻率均迭代100次，把低頻點反演結果作為高頻點反演的背景模型。與水平層狀模型一樣，背景模型也是真實模型經二維平滑處理得到的，如圖6b所示。

分別采用HFCG法和常規共軛梯度法對圖6所示模型進行反演，結果見圖7。可以看出，采用HFCG法的結果優于常規共軛梯度法反演結果，其反演結果分辨率更高。兩種方法反演結果與真實模型的L2范數殘差分別為49115和55991。圖8為模型縱向速度曲線對比，可以看出HFCG法反演結果更接近真實模型。

圖4 x=1000m(a)和x=1500m(b)處不同方法縱向速度曲線對比

圖5 不同頻率時HFCG法(紅線)與常規共軛梯度法(藍線)目標泛函收斂曲線對比

圖6 抽稀的Marmousi模型

圖7 采用HFCG法(a)和常規共軛梯度法(b)的Marmousi模型速度反演剖面

圖8 x=4000m(a)和x=5200m(b)處HFCG法和常規共軛梯度法反演縱向速度曲線對比

圖9為HFCG法和常規共軛梯度法在不同頻率的目標泛函收斂曲線對比。可以看出，每個頻率迭代100次后，HFCG法的目標泛函已經收斂，而常規共軛梯度法的目標泛函還未收斂，可見前者具有非常明顯的快速收斂優勢，只需很少的迭代次數就可實現目標泛函的收斂。

3.3 反演效率測試

為了驗證HFCG法的有效性，分別對水平層狀模型(圖2)和Marmousi模型(圖6)運用HFCG法和常規共軛梯度法進行反演效率定量測試，比較、分析在相同的迭代終止條件下兩種方法的迭代次數。本文采用的迭代終止條件是數據殘差

圖9 HFCG法(紅線)與常規共軛梯度法(藍線)在不同頻率的目標泛函收斂曲線對比

(25)

采用HFCG和常規共軛梯度法分別對圖2所示水平層狀模型進行反演，滿足式(25)時的反演結果見圖10。從圖中可以看出，當迭代誤差終止條件相同時，兩種方法反演得到水平層狀模型基本相同。HFCG法和常規共軛梯度法的反演迭代次數分別為122和271，前者較后者的計算效率提高了近55%。

采用HFCG和常規共軛梯度法分別對圖6所示Marmousi模型進行反演，滿足式(25)時的反演結果見圖11。從圖中可以看出，當迭代誤差終止條件相同時，兩種方法反演得到Marmousi模型基本相同。經統計，HFCG法和常規共軛梯度法反演的迭代次數分別為617和1183，前者較后者的計算效率提高了近49%。

圖10 水平層狀模型反演結果

圖11 Marmousi模型反演結果

4 結論

有限差分對比源反演利用有限差分算子構造反演算子，可求解逆散射問題，適用于復雜介質的對比源反演。運用該方法對單個頻率數據進行反演，正、反演算子只與背景模型有關，無需改變背景模型。因此只需構造一次正反演算子，進行一次LU分解，減少了正反演計算量，提高了反演效率。本文在常規共軛梯度法對比源反演的基礎上引入適用于對比源反演的混合快速共軛梯度法，在不增加計算量的基礎上加速目標泛函收斂，進一步提高了反演效率。基于HFCG法的FDCSI方法具有明顯加速目標泛函收斂的優勢，適用于大數據量反演。

本文只研究了介質密度為常數的頻率域聲波方程FDCSI方法。對于復雜地下介質，后續可針對密度參數，開展彈性波方程FDCSI的研究。