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

(①中國石化石油物探技術(shù)研究院，江蘇南京 211103； ②中國石油大學(xué)(北京)油氣資源與探測國家重點(diǎn)實(shí)驗(yàn)室，北京 102249)

0 引言

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

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

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

1 FDCSI基本理論

1.1 問題描述

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

式中：H和Hb分別為與模型和背景有關(guān)的稀疏阻抗矩陣；Sj為震源項(xiàng)。

圖1 逆散射示意圖

定義對(duì)比度χ(模型擾動(dòng))為[18]

(3)

式中：r為空間坐標(biāo)；k和c分別為模型的波數(shù)和速度；kb和cb分別為背景模型的波數(shù)和速度。

根據(jù)總波場、背景波場和散射波場的關(guān)系，可得[16]

(4)

(5)

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

(6)

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

根據(jù)式(4)可以得到狀態(tài)方程

(7)

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

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

(8)

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

(9)

(10)

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

1.2 有限差分對(duì)比源反演

在頻率域FDCSI過程中，需要同時(shí)滿足數(shù)據(jù)方程式(9)和狀態(tài)方程式(10)，不斷優(yōu)化求解對(duì)比源Wj和對(duì)比度χ，目標(biāo)泛函為[17-18]

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

(11)

定義S域和D域的L2范數(shù)為[17]

(12)

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

(13)

式中

(14)

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

(15)

(16)

(17)

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

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

(18)

(19)

(20)

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

2 混合快速共軛梯度法

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

(21)

(22)

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

(23)

當(dāng)目標(biāo)泛函殘差趨于收斂時(shí)，式(23)的迭代收斂曲線存在抖動(dòng)，因此在趨于收斂時(shí)采用常規(guī)共軛梯度法進(jìn)行優(yōu)化迭代，進(jìn)一步減小目標(biāo)泛函殘差。則HFCG的迭代公式為

(24)

式中Nunstable為迭代收斂曲線出現(xiàn)不穩(wěn)定現(xiàn)象時(shí)的迭代次數(shù)。

3 模型算例

為了驗(yàn)證方法的有效性，分別針對(duì)水平層狀模型和重采樣后的Marmousi模型進(jìn)行數(shù)值模擬實(shí)驗(yàn)。實(shí)驗(yàn)中把上一個(gè)頻率反演的結(jié)果作為下一個(gè)頻率反演的背景模型。

3.1 水平層狀模型

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

圖2 水平層狀速度模型

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

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

圖5為兩種方法6個(gè)頻率的目標(biāo)泛函收斂曲線對(duì)比。從圖中可以看出，HFCG法具有非常明顯的快速收斂優(yōu)勢，該方法在很少的迭代次數(shù)時(shí)就可實(shí)現(xiàn)目標(biāo)泛函的收斂。

3.2 Marmousi模型

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

分別采用HFCG法和常規(guī)共軛梯度法對(duì)圖6所示模型進(jìn)行反演，結(jié)果見圖7?？梢钥闯?，采用HFCG法的結(jié)果優(yōu)于常規(guī)共軛梯度法反演結(jié)果，其反演結(jié)果分辨率更高。兩種方法反演結(jié)果與真實(shí)模型的L2范數(shù)殘差分別為49115和55991。圖8為模型縱向速度曲線對(duì)比，可以看出HFCG法反演結(jié)果更接近真實(shí)模型。

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

圖5 不同頻率時(shí)HFCG法(紅線)與常規(guī)共軛梯度法(藍(lán)線)目標(biāo)泛函收斂曲線對(duì)比

圖6 抽稀的Marmousi模型

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

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

圖9為HFCG法和常規(guī)共軛梯度法在不同頻率的目標(biāo)泛函收斂曲線對(duì)比。可以看出，每個(gè)頻率迭代100次后，HFCG法的目標(biāo)泛函已經(jīng)收斂，而常規(guī)共軛梯度法的目標(biāo)泛函還未收斂，可見前者具有非常明顯的快速收斂優(yōu)勢，只需很少的迭代次數(shù)就可實(shí)現(xiàn)目標(biāo)泛函的收斂。

3.3 反演效率測試

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

圖9 HFCG法(紅線)與常規(guī)共軛梯度法(藍(lán)線)在不同頻率的目標(biāo)泛函收斂曲線對(duì)比

(25)

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

采用HFCG和常規(guī)共軛梯度法分別對(duì)圖6所示Marmousi模型進(jìn)行反演，滿足式(25)時(shí)的反演結(jié)果見圖11。從圖中可以看出，當(dāng)?shù)`差終止條件相同時(shí)，兩種方法反演得到Marmousi模型基本相同。經(jīng)統(tǒng)計(jì)，HFCG法和常規(guī)共軛梯度法反演的迭代次數(shù)分別為617和1183，前者較后者的計(jì)算效率提高了近49%。

圖10 水平層狀模型反演結(jié)果

圖11 Marmousi模型反演結(jié)果

4 結(jié)論

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

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