尹德強(qiáng)

中國電力工程顧問集團(tuán)東北電力設(shè)計院有限公司，長春 130021

0 引言

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

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

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

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

J=s-Wr2+μrTr

(1)

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

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

(2)

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

r=DL

(3)

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

(4)

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

J=‖s-WDL‖2+μDLTDL

(5)

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

(GTG+μI)·L=GTS

(6)

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

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

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

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

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

(7)

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

(8)

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

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

3 寬帶約束波阻抗反演

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

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

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ù)。

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

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

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

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

(12)

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

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

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

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

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

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

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

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

以一個61道砂泥巖互層理論速度模型模擬實際油氣儲層。該模型大致可分為5層：地表為速度較低的第四系覆蓋層；第二層為泥巖層；第三層為楔狀砂巖體儲層，局部含有油氣資源；而第四層泥巖層局部存在逆斷層，使得下部第五層含水砂巖向上侵入，各層的巖性及速度值見表1。選取的地震子波為阻尼余弦子波，地震道的長度為0.4 s，采樣點(diǎn)數(shù)為200個，采樣率為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è)各個巖層的密度均為常數(shù)，采用褶積模型合成地震記錄：阻尼余弦子波與理論速度模型相褶積合成人工地震記錄剖面，并在原始地震記錄中加入能量比為5%的隨機(jī)噪聲，則形成的擬實際地震剖面如圖6所示。在無噪音及5%隨機(jī)噪聲下，采用改進(jìn)的預(yù)條件共軛梯度法進(jìn)行寬帶約束波阻抗反演，反演重構(gòu)地下地質(zhì)模型，恢復(fù)含油儲層和斷層等地質(zhì)構(gòu)造信息，并反演出較為真實的速度模型。

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

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

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

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

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

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

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

5 結(jié)論

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

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

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