基于預(yù)條件共軛梯度法的寬帶約束波阻抗反演

尹德強(qiáng)

中國(guó)電力工程顧問(wèn)集團(tuán)東北電力設(shè)計(jì)院有限公司，長(zhǎng)春 130021

0 引言

寬帶約束反演是一種地震巖性反演方法，能夠充分利用原始地震資料提供地層的空間分布和斷層等有效信息，利用測(cè)井資料補(bǔ)償?shù)卣鹳Y料缺失的低頻和高頻信息，拓寬了地震資料的頻帶，有效解決地震數(shù)據(jù)的帶限問(wèn)題，恢復(fù)出寬頻帶波阻抗的空間分布[1]。地震數(shù)據(jù)和波阻抗參數(shù)之間呈現(xiàn)非線性數(shù)值關(guān)系，實(shí)質(zhì)要解決非線性反演問(wèn)題[2--3]。而非線性反演理論體系并不完善，反演過(guò)程中的收斂速度慢、解的穩(wěn)定性差及其反演精度低等存在的問(wèn)題常限制其應(yīng)用[4]。

針對(duì)求解雅克比系數(shù)矩陣病態(tài)問(wèn)題，研究者對(duì)共軛梯度算法改進(jìn)策略不同。采用正則化加入阻尼因子降低矩陣條件數(shù)；對(duì)系數(shù)矩陣采用預(yù)處理方法降低矩陣條件數(shù)；對(duì)搜索方向的改進(jìn)措施，采用P--R--P方法和DY方法等，但多數(shù)用于稀疏約束反褶積間接反演波阻抗。波阻抗反演存在帶限問(wèn)題，主要采取測(cè)井資料進(jìn)行約束，對(duì)于橫向阻抗急劇變化反演效果較差。筆者以地震資料解釋的標(biāo)志層作為控制點(diǎn)，依據(jù)測(cè)井資料提供的井點(diǎn)信息，以井點(diǎn)處進(jìn)行外推內(nèi)插，建立合適的寬頻帶初始波阻抗模型。以地震數(shù)據(jù)與合成地震記錄的殘差為基礎(chǔ)建立反演目標(biāo)函數(shù)，利用參數(shù)置換法對(duì)原始波阻抗進(jìn)行參數(shù)變換，反演解加入硬約束條件，對(duì)轉(zhuǎn)化的擬線性矩陣方程采用改進(jìn)的預(yù)條件共軛梯度算法進(jìn)行數(shù)值求解以獲得高分辨率波阻抗模型。

1 構(gòu)建反演目標(biāo)函數(shù)

采用阻尼最小二乘方差擬合函數(shù)建立地震反演目標(biāo)函數(shù)：地震數(shù)據(jù)與合成地震記錄的殘差為基礎(chǔ)，采用模型參數(shù)最小長(zhǎng)度解為約束：

J=s-Wr2+μrTr

(1)

式中：s為地震數(shù)據(jù)；r為反射系數(shù)序列；W為子波矩陣；μ為阻尼因子。

(GTG+μI)·ΔZ=GTΔS

(2)

針對(duì)式(2)求解波阻抗矩陣方程，可采用線性化算法進(jìn)行求解。但在實(shí)際迭代過(guò)程中，系數(shù)矩陣G每一次迭代會(huì)因初始模型的改變而發(fā)生變化，而系數(shù)矩陣G往往矩陣的條件數(shù)很大，呈現(xiàn)病態(tài)特征。同時(shí)搜索方向不一定是極值方向，會(huì)影響收斂速度和反演解的精度。依據(jù)參數(shù)置換法，設(shè)置波阻抗對(duì)數(shù)為新的參數(shù)，L(i)=ln[(Zi)]，反射系數(shù)與波阻抗映射關(guān)系為：

r=DL

(3)

式中：r=(r1,r2,…,rn)T,L=(L1,L2,…,Ln)T，為波阻抗對(duì)數(shù)矩陣；D為常數(shù)矩陣，形成的矩陣如式(4)所示：

(4)

原始問(wèn)題目標(biāo)函數(shù)式(1)可轉(zhuǎn)化為：

J=‖s-WDL‖2+μDLTDL

(5)

目標(biāo)函數(shù)式(5)在新參數(shù)波阻抗對(duì)數(shù)進(jìn)行泰勒級(jí)數(shù)展開，可得矩陣方程：

(GTG+μI)·L=GTS

(6)

式中：G為子波矩陣W與常數(shù)矩陣D形成的系數(shù)矩陣；S為地震數(shù)據(jù)矩陣；阻尼因子μ與地震數(shù)據(jù)和模型的方差矩陣有關(guān)。

由式(6)可知，在應(yīng)用共軛梯度法進(jìn)行迭代求解時(shí)，系數(shù)矩陣G為常數(shù)矩陣，合理的選擇阻尼因子μ，能夠降低系數(shù)矩陣的條件數(shù)，一定程度上減少在反演迭代過(guò)程中對(duì)收斂速度和搜索方向的影響。應(yīng)用傳統(tǒng)Fletcher--Reeves方法迭代求解式(6)矩陣方程時(shí)，算法后期迭代過(guò)程中，容易產(chǎn)生連續(xù)小步長(zhǎng)，收斂速度很慢。

2 改進(jìn)的預(yù)條件共軛梯度法

共軛梯度法是利用初始點(diǎn)處的梯度方向構(gòu)造一組共軛方向，沿著共軛方向進(jìn)行搜索目標(biāo)函數(shù)的極值，具有收斂速度快和二次終止性。傳統(tǒng)Fletcher--Reeves方法算法依賴于初始模型的選擇，易陷入局部極值。針對(duì)系數(shù)矩陣的病態(tài)性，常用的方法有正則化加入阻尼因子以降低矩陣條件數(shù)；對(duì)系數(shù)矩陣進(jìn)行不完全Cholesky分解、不完全的LU分解等降低系數(shù)矩陣條件數(shù)；采用重開始策略，令第n次迭代結(jié)果為新的初始點(diǎn)重開始，以使其最終達(dá)到收斂[6--7]。

采用Cauchy稀疏分布構(gòu)造反射系數(shù)約束項(xiàng)，建立誤差的最小二乘擬合函數(shù)[8--9]，矩陣方程為：

(7)

對(duì)式(7)進(jìn)行極小化并將矩陣方程進(jìn)行分解，預(yù)條件矩陣R以乘積的形式作用于系數(shù)矩陣為：

(8)

筆者采用一種改進(jìn)的近似預(yù)條件共軛梯度算法，令：G=WTW，d=WTs，改進(jìn)的措施是對(duì)負(fù)梯度d-Gr進(jìn)行預(yù)條件處理，具體的迭代步驟如下所示：

采用改進(jìn)的預(yù)條件共軛梯度法可改善系數(shù)矩陣的病態(tài)特征，提高收斂速度和反演精度，一定程度上可提高程序運(yùn)行效率。

3 寬帶約束波阻抗反演

寬帶約束波阻抗反演以原始地震數(shù)據(jù)資料為基礎(chǔ)，在先驗(yàn)知識(shí)的約束下，以地質(zhì)信息、測(cè)井資料為約束條件，建立寬頻帶的初始波阻抗模型；采用合適的數(shù)值方法求解，最終獲取具有寬頻帶的高分辨率波阻抗模型[10--12]。

采用改進(jìn)的預(yù)條件共軛梯度算法，使得搜索方向靠近負(fù)梯度方向，改進(jìn)的搜索方向?yàn)椋?/p>

p(k+1)=g(k+1)+βkp(k)

(9)

式中：p(k+1)為新的搜索方向；g(k+1)為負(fù)梯度方向；參數(shù)βk更改為：

(10)

從系數(shù)矩陣G本身性質(zhì)出發(fā)，其中?Ji/?Lj表示數(shù)據(jù)si在解分量Lj方向上的變化程度，則合成地震記錄J在解分量Lj總的變化率可以表示為：

(11)

式中：T為衰減因子，T=a·(k-1)；a為常數(shù)；k為迭代次數(shù)。

本文所采用的寬帶約束波阻抗反演主要步驟：

①由于地震資料屬于帶限信號(hào)，寬帶約束反演的關(guān)鍵是建立寬頻帶的初始波阻抗模型，補(bǔ)償缺失的低頻信息；以地震資料解釋的標(biāo)志層作為控制點(diǎn)，依據(jù)測(cè)井資料提供的井點(diǎn)信息，以井點(diǎn)處進(jìn)行外推內(nèi)插，建立合適的寬頻帶初始波阻抗模型?？稍诘卣鸬乐羞x取井資料進(jìn)行約束。

②在反演過(guò)程中，要對(duì)反演解進(jìn)行約束，不然會(huì)造成反演解偏離初始猜測(cè)的模型處太遠(yuǎn)，獲得的相對(duì)波阻抗期望特征相同；但絕對(duì)波阻抗值會(huì)出現(xiàn)嚴(yán)重誤差，不是根據(jù)初始猜測(cè)波阻抗模型推出數(shù)據(jù)趨勢(shì)得到。反演過(guò)程中對(duì)反演解加入硬約束：

Lmx(i)≤L(i)≤Lmz(i)

(12)

式中：Lmx(i)為約束反演解最小值；Lmz(i)為約束反演解最大值，取值的范圍可約束在偏離初始猜測(cè)值±15%。

③建立反演目標(biāo)函數(shù)，采用參數(shù)置換法對(duì)原始波阻抗進(jìn)行參數(shù)變換，得到式(6)近似線性系統(tǒng)；采用更改的預(yù)處理共軛梯度法求解，對(duì)初始波阻抗模型進(jìn)行迭代修改，直到求解的模型波阻抗合成的地震記錄與原始地震數(shù)據(jù)最佳匹配為止，進(jìn)而獲得最優(yōu)的高分辨率波阻抗模型，寬帶約束波阻抗反演流程見圖1。

圖1 寬帶約束波阻抗反演流程圖Fig.1 Inversion flow chart of broadband constrained wave impedance

4 模型試算及結(jié)果分析

設(shè)計(jì)一個(gè)12層單道理論速度模型及一個(gè)6層初始速度模型，相關(guān)系數(shù)僅為0.785(圖2)，采樣時(shí)間為2 ms，密度常數(shù)1.0，采樣點(diǎn)數(shù)235個(gè)。選取地震子波為零相位阻尼余弦子波，采樣率2 ms，子波的長(zhǎng)度為40 ms，主頻為60 Hz，地震記錄是由地震子波與反射系數(shù)相褶積形成，加入5%的隨機(jī)噪聲。

利用fortran語(yǔ)言編制計(jì)算機(jī)程序進(jìn)行波阻抗反演，其反演效果如圖3所示：

圖3 地震記錄中加入5%隨機(jī)噪聲反演結(jié)果對(duì)比Fig.3 Inversion results comparison with 5% random noise in seismic record

由圖3可知，在原始地震記錄中加入5%的隨機(jī)噪聲，反演結(jié)果與理論模型相關(guān)系數(shù)可達(dá)0.997，剖面與絕對(duì)波阻抗值均與理論值大體相同；對(duì)比采用稀疏約束反褶積法[13--14]，由迭代求出的反射系數(shù)存在誤差，小的反射系數(shù)誤差經(jīng)過(guò)合并后，遞推公式產(chǎn)生的波阻抗會(huì)產(chǎn)生很大累積誤差，造成深部波阻抗值偏離真實(shí)值很大。而對(duì)于矩陣求逆法、不完全Cholesky分解及不完全的LU分解，存在迭代過(guò)度，收斂不穩(wěn)定等問(wèn)題。筆者直接對(duì)波阻抗模型進(jìn)行反演，改進(jìn)搜索方向收斂速度更快，反演精度高，抗燥性強(qiáng)。

以一個(gè)61道砂泥巖互層理論速度模型模擬實(shí)際油氣儲(chǔ)層。該模型大致可分為5層：地表為速度較低的第四系覆蓋層；第二層為泥巖層；第三層為楔狀砂巖體儲(chǔ)層，局部含有油氣資源；而第四層泥巖層局部存在逆斷層，使得下部第五層含水砂巖向上侵入，各層的巖性及速度值見表1。選取的地震子波為阻尼余弦子波，地震道的長(zhǎng)度為0.4 s，采樣點(diǎn)數(shù)為200個(gè)，采樣率為2 ms，理論速度模型如圖4所示。依據(jù)第35道速度值建立虛擬井?dāng)?shù)據(jù)，構(gòu)建與理論模型相關(guān)性較差的初始速度模型(圖5)。

表1 理論模型速度參數(shù)Table 1 Speed parameters of theoretical model

圖4 砂泥巖互層巖性油氣藏模型Fig.4 Lithologic reservoir model in sand-shale interbed

圖5 構(gòu)建速度初始模型Fig.5 Initial model of established velocity

設(shè)各個(gè)巖層的密度均為常數(shù)，采用褶積模型合成地震記錄：阻尼余弦子波與理論速度模型相褶積合成人工地震記錄剖面，并在原始地震記錄中加入能量比為5%的隨機(jī)噪聲，則形成的擬實(shí)際地震剖面如圖6所示。在無(wú)噪音及5%隨機(jī)噪聲下，采用改進(jìn)的預(yù)條件共軛梯度法進(jìn)行寬帶約束波阻抗反演，反演重構(gòu)地下地質(zhì)模型，恢復(fù)含油儲(chǔ)層和斷層等地質(zhì)構(gòu)造信息，并反演出較為真實(shí)的速度模型。

圖6 加入5%隨機(jī)噪聲的合成地震記錄Fig.6 Synthetic seismic record with 5% random noise

如圖7所示，無(wú)噪聲情況下，反演得到的速度模型較好反映出楔狀含油儲(chǔ)層和下部地層存在的逆斷層位置，準(zhǔn)確描述出各個(gè)反射層位，相關(guān)系數(shù)可達(dá)0.998，各道反演速度點(diǎn)的絕對(duì)誤差最大僅為真實(shí)值的±4%。

圖7 無(wú)噪音反演得到的速度模型Fig.7 Velocity model obtained from inversion without noise

如圖8所示，存在5%的隨機(jī)噪聲干擾時(shí)，反演得到速度模型也能夠描述出儲(chǔ)層和逆斷層的位置，由于噪音強(qiáng)度分布不同，造成各道反演得到的絕對(duì)速度值存在偏差，但最大偏差值僅為真實(shí)值的±5%?；趯拵Ъs束波阻抗反演研究方法對(duì)比中，筆者依據(jù)系數(shù)矩陣的特征巧妙構(gòu)造預(yù)處理矩陣，既不改變?cè)汲跏紗?wèn)題，又能利用預(yù)處理技巧；在迭代過(guò)程中，改進(jìn)搜索方向，對(duì)反演解進(jìn)行加入硬約束條件，充分利用測(cè)井信息，使得反演解是由初始猜測(cè)的波阻抗模型推出的數(shù)據(jù)趨勢(shì)而得到。

圖8 5%隨機(jī)噪音反演得到的速度模型Fig.8 Velocity model obtained from inversion with 5% random noise

由圖9可知，基于預(yù)條件的寬帶約束波阻抗反演較大的偏差均集中在15～25道之間，是由于給出初始模型相關(guān)性較差的原因，以單道進(jìn)行反演，沒(méi)有考慮相鄰道相互約束，造成相鄰道數(shù)據(jù)存在不同誤差。通過(guò)擬實(shí)際數(shù)據(jù)驗(yàn)證，采用的改進(jìn)的預(yù)條件共軛梯度算法應(yīng)用于寬帶約束反演方法，驗(yàn)證反演算法的正確性和優(yōu)越性。

圖9 反演結(jié)果的平均誤差曲線Fig.9 Average absolute error curves of inversion results

5 結(jié)論

(1)寬帶約束反演建立優(yōu)化的初始波阻抗模型，可以減少多解性問(wèn)題，依賴于地震資料的品質(zhì)和高信噪比，適用于井資料豐富的地質(zhì)區(qū)域。

(2)用改進(jìn)的預(yù)處理共軛梯度法，對(duì)反演解加入硬約束條件，實(shí)現(xiàn)基于模型反演的寬帶約束波阻抗反演。在5%隨機(jī)噪聲干擾時(shí)，反演最大偏差值僅為5%，且真實(shí)反映儲(chǔ)層和逆斷層位置。筆者采用的反演算法收斂速度快、反演精度高、數(shù)值穩(wěn)定性好，具有一定的抗噪性。

(3)線性化的數(shù)值算法進(jìn)行迭代求解，反演的結(jié)果常與初始模型選擇有關(guān)，通過(guò)復(fù)雜理論模型試算，偏差值較大集中在15～25道，收斂到局部極值，未考慮相鄰道約束。因此，發(fā)展和完善完全非線性反演理論是地震反演方法發(fā)展必然趨勢(shì)。