全域正則化快速匹配追蹤稀疏地震反演方法

裴 松 印興耀 李 坤

(中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島 266580)

0 引言

地球物理學(xué)家希望通過觀測數(shù)據(jù)反推得到地下介質(zhì)模型參數(shù)，即建立觀測數(shù)據(jù)與地下介質(zhì)參數(shù)之間的關(guān)系。地震勘探中，地震記錄可表示為地震子波與反射系數(shù)的褶積，也可表述為子波矩陣與反射系數(shù)向量的乘積。地震反演在無噪情況下可利用簡單的矩陣求逆獲得目標(biāo)泛函的解。然而，地震記錄通常是包含噪聲的，這便需要利用最小二乘等方法獲得最優(yōu)解[1-2]。地震反演技術(shù)自提出以來便被廣泛應(yīng)用于烴類檢測、儲層預(yù)測等領(lǐng)域[3-6]。隨著地震勘探技術(shù)的發(fā)展與進步，對反演精度與分辨率提出了更高要求[7-8]。稀疏地震反演除了可獲得稀疏反射系數(shù)及“塊化”阻抗之外，還可有效提高分辨率、橫向連續(xù)性及層位邊界保真度[9-10]。

匹配追蹤(Matching Pursuit，MP)算法最早是由Mallat等[11]提出，用以進行信號分解及高精度時頻分析的一種方法。近年來，MP被廣泛應(yīng)用于信號稀疏表示、高精度時頻分析、地震信號分解與重構(gòu)[12]、噪聲與強反射壓制[13-16]、儲層預(yù)測[17]、薄層頂?shù)追瓷湎禂?shù)提取[18]及砂體尖滅點刻畫[19]等眾多領(lǐng)域。此外，MP是求解L0范數(shù)的經(jīng)典算法，可有效獲取目標(biāo)泛函的稀疏解[20]。然而，MP是一種貪婪的迭代算法，因此如何提高MP算法的計算效率一直是業(yè)內(nèi)的熱點問題。針對MP算法計算效率低等問題，眾多MP的發(fā)展算法被提出，三參數(shù)快速MP可以通過縮減冗余字典的規(guī)模提高計算效率[21]；動態(tài)MP通過改變時頻原子局部瞬時屬性提高信號分解速度[22-23]；正則化正交MP可在字典不完備的情況下實現(xiàn)有效的稀疏信號恢復(fù)，且該算法具備良好的魯棒性[24]。雖然在信號分析領(lǐng)域MP的發(fā)展算法層出不窮，但在地震反演領(lǐng)域，如何提高MP反演算法計算效率與穩(wěn)定性仍是熱點問題。

時間域地震反演具有穩(wěn)定性高、方法簡單可行的特點[25-28]，但其分辨率較低。頻率域反演通過構(gòu)建頻率域褶積模型，可以進行地震信號頻帶內(nèi)的頻率解耦，從而獲得高分辨率反演結(jié)果，但其穩(wěn)定性不足[20]。此外，地震反演是病態(tài)、多解的，反演多解性主要表現(xiàn)為多個地下參數(shù)模型都可產(chǎn)生相同的地震響應(yīng)[29-30]。在反演框架中添加如高斯分布等先驗信息是緩解反演多解性的有效方法之一，然而這種方式無法獲得反射系數(shù)稀疏解，也就無法獲得“塊化”阻抗信息[10]。

本文首先從MP冗余字典的構(gòu)建出發(fā)，結(jié)合時間域反演穩(wěn)定性好、頻率域反演分辨率高的特性，構(gòu)建時頻域冗余字典。對于MP算法，本文提出全域正則化快速MP算法，每次迭代同時選擇所有可能的匹配原子，并通過正則化方法依據(jù)能量篩選出最優(yōu)原子作為最終的匹配原子庫。在此基礎(chǔ)上，將初始模型約束引入反演框架，構(gòu)建初始模型約束的時頻域地震反演目標(biāo)泛函，并通過最小二乘求解。最后，利用模型及實際地震資料測試了本文方法的有效性及實用性。模型及實際資料測試結(jié)果表明：本文方法可以有效獲取反射系數(shù)的稀疏解。此外，相較于常規(guī)MP反演方法，本文方法的計算效率、反演結(jié)果的分辨率及邊界保真度均得到了有效提升。

1 理論與方法

1.1 初始模型約束的時頻域地震反演

在地震反演領(lǐng)域，經(jīng)典正演模型可表示為

S(t)=W·R+N

(1)

式中S(t)、W、R、N、t分別表示時間域地震信號、地震子波矩陣、反射稀疏序列、噪聲序列及時間。針對時間域穩(wěn)定性高、但分辨率有限的問題，聯(lián)合域褶積模型被提出[31]

S(f)=F(f)·R+N(f)

(2)

式中S(f)、F(f)、N(f)分別為地震信號頻譜、頻率域核矩陣及噪聲頻譜，f表示頻率。具體地，式(2)中時頻域核矩陣的公式為

F(f)=

(3)

式中：fi為所選頻率,i=1,2,…,n；地震信號S(t)的采樣點數(shù)為m;W(2πfi)為子波頻率域形式。一般地，地震信號的有效頻帶為10～60Hz，而對于寬頻地震信號，該范圍的上、下限應(yīng)該由地震信號分析獲得。本文所用為常規(guī)地震信號，因此頻帶設(shè)定為10～60Hz。由于信號頻譜是由虛部及實部兩部分組成，因此可將時頻域聯(lián)合褶積模型改寫為

(4)

式中：下標(biāo)r代表信號實部；下標(biāo)i表示信號虛部。將式(4)改寫為簡單的矩陣形式

O=DR(t)+n

(5)

為了提高反演結(jié)果的橫向連續(xù)性，降低反演多解性，縱波阻抗初始模型約束被引入至反演目標(biāo)泛函。需要注意的是，初始模型并非低頻模型?？v波阻抗低頻模型為

(6)

式中Plow為m×1階向量?？v波阻抗初始模型為

(7)

聯(lián)合時頻域褶積模型與初始模型約束，可得

(8)

式中：C為m×m階下三角陣，也即積分矩陣；Pini為m×1階向量。

1.2 全域正則化快速匹配追蹤

MP是一種貪婪的迭代尋優(yōu)方法，通過遍歷冗余字典逐次尋找最優(yōu)解[32]。MP主要包含以下幾個步驟：①依據(jù)信號本身構(gòu)建相應(yīng)的冗余字典；②求取信號與冗余字典的內(nèi)積，也即投影；③尋找內(nèi)積最大值的位置并在字典中搜索對應(yīng)的最佳匹配原子；④利用原信號減去最佳匹配原子獲得殘差；⑤將殘差作為新信號重復(fù)①～④，直至殘差小于設(shè)定閾值。對于地震信號，經(jīng)過次迭代后的地震信號可表示為

(9)

式中：AIte、χIte、d分別為經(jīng)第Ite次迭代所得最佳匹配原子的振幅、最佳匹配原子及殘差；k為迭代次數(shù)。

以下將著重介紹全域正則化快速匹配追蹤算法。首先，計算信號與冗余字典的內(nèi)積

Pro=〈Sig,H〉

(10)

式中Sig與H分別表示時頻域與縱波阻抗初始模型計算的對應(yīng)信號與冗余字典，相應(yīng)表達(dá)式為

(11)

式中：Sig為(2n+m)×1階向量；H為(2n+m)×m階矩陣；Pro為m×1階向量；α1與α2分別為時頻域地震反演與初始模型約束的權(quán)重。α1一般取1即可；α2越大，反演結(jié)果越向初始模型靠攏。與常規(guī)多原子匹配追蹤不同，本文方法在每次迭代中同時選取所有極大值位置點χ作為匹配原子備選庫

χ=localmax[Pro]

(12)

在此基礎(chǔ)上，需要對備選原子庫進行正則化篩選，在所有子集中依據(jù)能量篩選出能量最大的匹配原子庫

J0∈χs.t. |ui|≤2|uj|i,j∈J0

(13)

式中u、i和j的設(shè)定及具體正則化過程如圖1所示。式(13)所篩選出的原子庫即為在原子庫χ的所有子集中具有最大能量的子集J0，因此在某次迭代后最終的原子庫為Loc=Loc∪J0(Loc在初次迭代前可設(shè)置為空集φ，詳情可參考流程圖中步驟二)。相較于常規(guī)快速匹配追蹤反演中人為采用最大內(nèi)積點的百分比函數(shù)進行計算，基于正則化的全域快速匹配追蹤算法可更好地篩選出全域字典中所有的符合條件的原子，在提高計算效率的同時有效提高了反演的準(zhǔn)確性。在獲得匹配原子庫后，可以利用最小二乘來獲取對應(yīng)振幅

(14)

此時，殘差信號可表示為

(15)

得到殘差后，令Sig=Residual，并重新進行正則化篩選與反演，直至迭代次數(shù)達(dá)到設(shè)定值或殘差小于閾值。圖1給出了本文方法的詳細(xì)流程。

圖1 全域正則化快速匹配追蹤稀疏地震反演方法流程圖

1.3 一維模型試算

為了驗證本文方法的可行性，本文從測井資料中抽取了縱波速度與密度來合成縱波阻抗，并利用縱波阻抗計算得到反射系數(shù)。采用主頻為30Hz的雷克子波與反射系數(shù)褶積得到了如圖2a所示的合成地震記錄，其采樣間隔為2ms，采樣點數(shù)為250。若要獲得稀疏解，需要降低迭代次數(shù)。

圖2為利用本文方法獲得的稀疏解，由圖可知，本方法可以很好地進行稀疏反射系數(shù)反演，且稀疏度與反射系數(shù)稀疏度較為一致，塊化縱波阻抗的反演結(jié)果也與模型值高度相似。需說明的是本文中所有一維測試結(jié)果中的黑色、紅色及藍(lán)色曲線均分別代表模型值、反演值和初始模型值。

圖2 無噪地震記錄稀疏反演結(jié)果

為了驗證本文方法的有效性與實用性，分別利用本文方法與常規(guī)方法(快速匹配追蹤地震反演方法)計算得到了如圖3、圖4所示的非稀疏反演結(jié)果。除迭代次數(shù)外，其余參數(shù)設(shè)置相同。其中，本文方法僅需15次迭代即可得到非稀疏反演結(jié)果(圖3c)，而常規(guī)方法則需要50次迭代得到圖4c所示的非稀疏反演結(jié)果。此外，本文方法與傳統(tǒng)方法用時分別為0.428s和1.265s。圖3、圖4中縱波阻抗反演結(jié)果與真實縱波阻抗相關(guān)系數(shù)分別為0.98和0.96。因此，本文方法反演所得結(jié)果在細(xì)節(jié)處與模型的匹配度稍高于常規(guī)方法所得反演結(jié)果。由對比結(jié)果可知，本文方法在有效提高計算效率的同時還可獲得更準(zhǔn)確的反演結(jié)果。

圖3 本文方法無噪地震記錄非稀疏反演結(jié)果

圖4 常規(guī)方法無噪地震記錄非稀疏反演結(jié)果

為測試本文方法的魯棒性，在合成地震記錄中加入了不同的隨機高斯噪聲，圖5a、圖6a、圖7a分別展示了信噪比(SNR)分別為5、2、1的含噪地震記錄。由圖可知，雖然隨著信噪比的降低，反演值與模型值有所差異，但即使在信噪比為1的情況下，反演結(jié)果與模型值的相對變化吻合度依然較高。

圖5 SNR=5的反演結(jié)果

圖7 SNR=1的反演結(jié)果

此外，筆者還對本文所提方法進行了質(zhì)控測試(圖8)。所謂質(zhì)控，即固定某些參數(shù)，只改變某個或部分參數(shù)來測試該參數(shù)對反演結(jié)果的影響。圖8a為信噪比為1的地震記錄，后續(xù)質(zhì)控測試都是在此地震記錄之上進行的。首先，固定α1與α2，改變迭代次數(shù)來測試迭代次數(shù)(同圖1“初始設(shè)定”中的迭代次數(shù))對反演結(jié)果的影響。圖8b～圖8e分別展示了迭代次數(shù)為30、20、10、9次的縱波阻抗反演結(jié)果。由圖可知：對于含有較強噪聲的地震記錄，隨著迭代次數(shù)的降低，反演結(jié)果逐漸趨于穩(wěn)定。在此基礎(chǔ)上，固定迭代次數(shù)(9次)不變，通過改變初始模型權(quán)重參數(shù)(α2)來測試初始模型約束對于反演結(jié)果的影響。圖8f～圖8i的迭代次數(shù)都為9次，而初始模型權(quán)重參數(shù)分別為0.6、0.7、0.9、1.2。由圖8可知，當(dāng)減小迭代次數(shù)后，增加初始模型權(quán)重可以令反演結(jié)果靠近初始模型，從而進一步穩(wěn)定反演結(jié)果。因此，當(dāng)?shù)卣鹦盘栔泻羞^強噪聲(SNR<2)時，需要減小迭代次數(shù)同時提高初始模型約束項的權(quán)重(α2)獲取穩(wěn)定反演結(jié)果。綜上，無噪及含噪地震信號的測試都說明了本文方法的可行性和有效性。

圖8 質(zhì)控測試

2 二維模型試算

逆掩斷層模型是測試稀疏反演方法邊界保真度的有效模型。因此，為了進一步驗證方法有效性，本文利用Marmousi逆掩斷層二維模型對全域正則化快速匹配追蹤算法進行測試。針對合成地震記錄圖9a及原始縱波阻抗(圖9b)，利用式(7)可獲得對應(yīng)的初始模型圖9c。分別采用常規(guī)快速匹配追蹤算法及本文方法對該逆掩斷層模型進行有效性測試。圖9d與圖9e分別為常規(guī)方法與本文方法的測試結(jié)果。由對比可知，利用全域正則化快速匹配追蹤方法得到的縱波阻抗反演結(jié)果相較于常規(guī)方法得到的縱波阻抗反演結(jié)果具有更好的層位邊界保真度。尤其在黑色箭頭與黑色方框所示區(qū)域，可見本文方法所得結(jié)果具有清晰的邊界與更好的橫向連續(xù)性。

圖9 二維模型測試結(jié)果

本文還針對二維模型進行了抗噪性測試，圖10a與圖10b分別展示了SNR=1的地震剖面與反演結(jié)果。由圖可知，反演結(jié)果的橫向連續(xù)性有些下降(圖10b中箭頭所示)，但其邊界仍然較為清晰，且反演結(jié)果依然與模型保持了高度一致。

圖10 抗噪性測試結(jié)果

此外，仍然對二維模型反演的用時進行了記錄。常規(guī)方法與本文方法用時分別為118s和35s。因此，本文方法的計算效率要高于傳統(tǒng)快速匹配追蹤算法，該測試同樣驗證了全域正則化快速匹配追蹤稀疏地震反演方法的有效性。

3 實際資料應(yīng)用

在模型測試的基礎(chǔ)上，還需要利用實際資料來對本文方法的實用性進行測試。從中國東部M工區(qū)抽取了一個由18600道組成的三維地震數(shù)據(jù)體，測試本文方法的效果，并抽出兩個剖面進行測試結(jié)果展示。從圖11a所示的地震剖面可知，該工區(qū)地層橫向連續(xù)性較差，且存在多處弱反射區(qū)域。此外，快速的縱向變化及斷層的發(fā)育都給地震反演方法帶來了一定的困難與挑戰(zhàn)。

為了進行更加精確的反演，需利用高頻濾波后的井?dāng)?shù)據(jù)建立縱波阻抗低頻模型。圖11b為縱波阻抗低頻模型，利用式(7)即可獲得相應(yīng)的縱波阻抗初始模型。圖11c與圖11d為采用全域正則化快速匹配追蹤稀疏地震反演方法獲得的稀疏反射系數(shù)反演結(jié)果與塊化縱波阻抗反演結(jié)果。由圖11c可以看出，稀疏反射系數(shù)反演結(jié)果具有良好的橫向連續(xù)性及稀疏度，其分辨率相對較高。由圖11d可知，塊化縱波阻抗反演結(jié)果的層位邊界較為清晰，可以有效對地層邊界進行解釋與刻畫。為了說明方法的準(zhǔn)確性，從三維反演結(jié)果中抽取一個剖面與井進行對比(圖12)。由于井?dāng)?shù)據(jù)頻率較高，對實際井?dāng)?shù)據(jù)進行高頻濾波，濾掉60Hz以上的頻率成分，以與地震數(shù)據(jù)保持一致。由圖12可知，利用本文方法所得的縱波阻抗反演結(jié)果與井保持了高度一致，也說明了本文方法的有效性與實用性。

圖11 三維地震數(shù)據(jù)反演結(jié)果

圖12 井震對比

此外，為了進一步驗證方法有效性，分別利用本文方法與常規(guī)方法對該工區(qū)另一剖面(圖13a)進行測試。圖13b與圖13c分別為本文方法及常規(guī)方法(快速匹配追蹤地震反演方法)縱波阻抗反演結(jié)果。由圖可知，本文方法反演結(jié)果橫向連續(xù)性更高(白色箭頭所指區(qū)域)。本文方法與常規(guī)方法的迭代次數(shù)分別為8次(用時3.8s)與30次(用時10.2s)。因此，相較于常規(guī)方法，本文方法可在更少的迭代次數(shù)及耗時的情況下獲得更穩(wěn)定精確的反演結(jié)果。

圖13 不同方法縱波阻抗反演結(jié)果對比

4 結(jié)論

(1)全域正則化快速匹配追蹤稀疏地震反演方法首先通過時頻域褶積模型與縱波阻抗約束構(gòu)建匹配追蹤冗余字典。在此基礎(chǔ)上利用對全域信號進行整體分析得出所有可能的匹配原子構(gòu)建備選原子庫，最后利用正則化方法篩選出能量最大子集構(gòu)成最終的匹配原子庫進行反演。

(2)在利用全域正則化快速匹配追蹤稀疏地震反演方法時，需要控制迭代次數(shù)。當(dāng)?shù)卣鹳Y料包含較強噪聲時，需要降低迭代次數(shù)、提高初始模型約束權(quán)重來獲得穩(wěn)定解。相反，若地震資料信噪比較高，則可以提高迭代次數(shù)獲得較為精確的反演結(jié)果。

(3)模型與實際資料測試表明，本文方法在保證反演準(zhǔn)確性與魯棒性的基礎(chǔ)上有效地提高了匹配追蹤稀疏反演方法的計算效率。通過對比由本文方法與常規(guī)方法獲得的縱波阻抗反演結(jié)果可知：相較于常規(guī)方法，本文方法反演結(jié)果的橫向連續(xù)性及邊界保真度均有所提高。