陳 剛 齊洪巖 李 維 張 融 鮮成鋼 王振林

(①中國石油大學(xué)(北京)非常規(guī)油氣科學(xué)技術(shù)研究院，北京 102249； ②中國石油新疆油田分公司勘探開發(fā)研究院，新疆克拉瑪依 834000)

0 引言

現(xiàn)今油氣勘探方向已逐步由構(gòu)造油氣藏轉(zhuǎn)向巖性、非常規(guī)等油氣藏，而這類油氣儲(chǔ)層多表現(xiàn)為縱向薄、橫向變化快的特征，這就要求現(xiàn)場采集到較寬頻、高頻地震記錄，且后期處理采用反Q濾波、反褶積等一系列高頻補(bǔ)償技術(shù)。然而，地震信號的高頻成分往往含較多噪聲，顯著降低地震資料的橫向連續(xù)性和信噪比，嚴(yán)重影響數(shù)據(jù)的分析和處理效果。去除噪聲、獲取高信噪比地震記錄一直是地震勘探要解決的首要問題[1-2]，目前地震數(shù)據(jù)處理中主要利用有效信號與噪聲的頻率或視速度差異壓制噪聲[3-4]，常見方法有中值濾波[5-6]、奇異值分解[7-8]、f-x域預(yù)測濾波和小波變換等[9-11]。

卡爾曼濾波于1960年由美國控制論專家Kalman[12]首次提出，是一種最優(yōu)化自回歸數(shù)據(jù)處理算法，即以最小均方差為估計(jì)準(zhǔn)則，采用狀態(tài)轉(zhuǎn)移反映系統(tǒng)變化的內(nèi)在規(guī)律，以“預(yù)估—測量—校正”的順序遞推，從被污染系統(tǒng)中估計(jì)出系統(tǒng)本來面目?？柭鼮V波算法原理表明其在有效信號恢復(fù)方面具有很大優(yōu)勢，已在目標(biāo)跟蹤、自動(dòng)控制、導(dǎo)航、氣象、空間物理、電力系統(tǒng)、醫(yī)學(xué)、通信等領(lǐng)域得到廣泛應(yīng)用[13-20]。在地震勘探領(lǐng)域，Crump[21]最先將其應(yīng)用于反褶積處理以得到反射系數(shù)序列。隨后，很多學(xué)者在此基礎(chǔ)上也將卡爾曼濾波用于地震資料反褶積處理，都取得很好效果。

董敏煜等[22]在Crump[21]和Mendel[23]研究成果的基礎(chǔ)上，進(jìn)一步提出了基于卡爾曼濾波估計(jì)時(shí)變子波的自適應(yīng)反褶積方法，顯著地提高了地震資料的分辨率、信噪比和同相軸連續(xù)性。鄧小英等[24]提出一種新的基于卡爾曼濾波的地震記錄同相軸跟蹤方法?；魢鴹澋萚25-26]通過在貝葉斯約束反演中建立卡爾曼濾波模式，實(shí)現(xiàn)縱向和橫向同時(shí)約束AVO反演，大大提高了反演結(jié)果的信噪比和橫向連續(xù)性。鄭笑雪等[27]通過卡爾曼濾波實(shí)現(xiàn)反演的橫向約束，進(jìn)而提高反演結(jié)果的信噪比和橫向連續(xù)性，同時(shí)表明卡爾曼濾波在提高地震資料信噪比和同相軸連續(xù)性方面確有明顯優(yōu)勢。

針對地震資料中存在的信噪比低和橫向連續(xù)性差的問題，本文基于卡爾曼濾波原理建立迭代循環(huán)組合濾波模式，將橫向(地震道)與縱向(時(shí)間)迭代循環(huán)濾波次數(shù)設(shè)定為4∶1，實(shí)現(xiàn)地震資料高效濾波，該方法在濾波效果和效率上都優(yōu)于一般平滑濾波、單一橫向或縱向迭代循環(huán)濾波。通過Marmousi2理論模型[28]和實(shí)際資料測試，處理結(jié)果驗(yàn)證了該方法的可靠性和實(shí)用性。

1 方法原理

1.1 橫向卡爾曼地震濾波

卡爾曼濾波分為預(yù)測和校正兩個(gè)步驟，預(yù)測是基于上一時(shí)刻狀態(tài)估計(jì)當(dāng)前時(shí)刻狀態(tài)，而校正則是綜合當(dāng)前時(shí)刻的估計(jì)狀態(tài)與觀測狀態(tài)，估計(jì)出最優(yōu)狀態(tài)[24]。本文依據(jù)卡爾曼濾波原理[12]建立卡爾曼地震道濾波模式。

目標(biāo)地質(zhì)體通常具有的橫向相似性導(dǎo)致相鄰地震道的高相似性，那么據(jù)相鄰地震道可估計(jì)本道最優(yōu)狀態(tài)。首先建立地震道之間的狀態(tài)轉(zhuǎn)換方程

(1)

再建立地震道的測量方程

(2)

根據(jù)上述狀態(tài)方程可優(yōu)化地震道數(shù)據(jù)Xk，對應(yīng)的Xk道的標(biāo)準(zhǔn)差Pk的更新公式為

(3)

式中：Pk、Pk-1分別是第k和第k-1道地震道的標(biāo)準(zhǔn)差；Rk為Pk更新過程中的誤差標(biāo)準(zhǔn)差。

有了地震道狀態(tài)的預(yù)測結(jié)果和地震道的測量值，就可對Xk進(jìn)行最優(yōu)化估算

(4)

(5)

式中Qk為系統(tǒng)過程的標(biāo)準(zhǔn)差。

到此，實(shí)施了對Xk的最優(yōu)估計(jì)。為了進(jìn)一步對Xk+1進(jìn)行最優(yōu)估計(jì)，需得到更新后的Pk

Pk=Pk(1-gk)

(6)

基于卡爾曼濾波原理建立的以上濾波模式，在確定P0、Rk和Qk后就可對地震資料進(jìn)行橫向?yàn)V波。這里，P0取開始濾波的第1道的標(biāo)準(zhǔn)差，而針對沉積巖層，相鄰地震道之間的差異可近似看成平穩(wěn)的隨機(jī)高斯白噪聲，故Pk、Qk可分別由下式求取

(7)

Qk=Ak[std(Xk-Xk-1)+std(Xk+1-Xk)]

(8)

式中std(·)表示取標(biāo)準(zhǔn)差。

1.2 縱向卡爾曼地震濾波

類似于上述橫向?yàn)V波。依據(jù)卡爾曼濾波模式，首先建立地震道相鄰時(shí)刻之間的狀態(tài)轉(zhuǎn)換方程

(9)

再建立地震數(shù)據(jù)t時(shí)刻的測量方程

(10)

有了狀態(tài)方程，就可優(yōu)化t時(shí)刻的地震數(shù)據(jù)Xt，那么對應(yīng)Xt的標(biāo)準(zhǔn)差Pt的更新公式為

(11)

式中：Pt、Pt-1分別是t和t-1時(shí)刻地震數(shù)據(jù)的標(biāo)準(zhǔn)差；Rt為Pt更新過程中的誤差標(biāo)準(zhǔn)差。

基于t時(shí)刻地震數(shù)據(jù)的預(yù)測結(jié)果和對應(yīng)的測量值，就可對Xt進(jìn)行最優(yōu)化估算

(12)

(13)

式中Qt為系統(tǒng)過程的標(biāo)準(zhǔn)差。

到此，完成了對t時(shí)刻的Xt的最優(yōu)估計(jì)。為了對Xt+1進(jìn)行最優(yōu)估計(jì)，需得到更新的Pt

Pt=Pt(1-gt)

(14)

基于卡爾曼濾波原理建立的以上濾波模式，在確定P0、Rt和Qt后就可對地震資料進(jìn)行縱向?yàn)V波。此處P0取濾波開始時(shí)刻的標(biāo)準(zhǔn)差，而針對沉積巖層，地震相鄰時(shí)刻之間的差異可近似看成平穩(wěn)的隨機(jī)高斯白噪聲，因此Rt、Qt分別由下式求取

(15)

Qt=At[std(Xt-Xt-1)+std(Xt+1-Xt)]

(16)

1.3 橫縱向組合的卡爾曼地震資料濾波

倘若只考慮橫向?yàn)V波，不僅會(huì)一定程度拉平地震同相軸，而且會(huì)在增加道與道之間相似程度的同時(shí)降低地震的橫向分辨率，導(dǎo)致地質(zhì)體橫向邊界不清晰。為此，本文在對地震剖面進(jìn)行橫向卡爾曼濾波后再做縱向卡爾曼濾波以約束地震同相軸的縱向真實(shí)位置。具體方法為： 實(shí)施一次上述橫向?yàn)V波模式稱為一次迭代循環(huán)濾波； 更新地震剖面，針對更新剖面進(jìn)行下次迭代循環(huán)濾波，繼續(xù)對地震剖面橫向進(jìn)行4N(N為整數(shù))次迭代循環(huán)濾波； 再從時(shí)間(縱向)上也進(jìn)行N次迭代循環(huán)濾波；如此以橫縱向迭代循環(huán)濾波次數(shù)比為4∶1進(jìn)行M輪(次)大循環(huán)濾波，這樣就可實(shí)現(xiàn)地震資料的橫縱向?yàn)V波。顯然，總濾波次數(shù)為5MN次。

2 方法測試

2.1 模型測試

為測試本文所提濾波方法的有效性和可靠性， 從Marmousi2模型(圖1a)中選用含有油、氣和水的砂巖部分模型(圖1b)。將由Zoeppritz方程計(jì)算的反射系數(shù)與主頻為35Hz的雷克子波進(jìn)行褶積，得到圖2a所示合成記錄剖面，可見剖面上氣藏和油藏的反射特征和橫向邊界較清楚； 再對該數(shù)據(jù)加入信噪比是1∶2的高斯隨機(jī)噪聲(圖2b)，可見加噪后剖面上油藏反射特征幾乎被噪聲淹沒。圖3a和圖3b分別為二維中值濾波1次和100次后結(jié)果，其噪聲得到一定壓制，但整體信噪比依然較低，且濾波次數(shù)增加對改善噪聲壓制效果并不明顯。圖4a和圖4b分別是單一橫向卡爾曼濾波100次和200次后結(jié)果，剖面同相軸、油藏和氣藏得到一定的恢復(fù)，噪聲被壓制，但出現(xiàn)較多水平層狀反射。圖5a和圖5b分別為單一縱向卡爾曼濾波25次和50次剖面顯示，可見縱向?yàn)V波去噪效果不佳，橫向連續(xù)性改善甚微，且濾波超過50次后同相軸出現(xiàn)失真。圖6a和圖6b分別為本文方法濾波50次和100次結(jié)果，可見該方法較好地恢復(fù)了有效信號并合理壓制了噪聲，同相軸反射特征也得到較好恢復(fù)； 同時(shí)，本文方法濾波效果和效率都優(yōu)于單一橫向或縱向?yàn)V波方法。

圖1 Marmousi2模型(a)及選取的(綠框)測試部分(b)[1]

圖2 原始合成(a)及加入1∶2高斯隨機(jī)噪聲(b)剖面

圖3 二維中值濾波1次(a)和100次(b)后剖面對比

圖4 單一橫向卡爾曼濾波100次(a)和200次(b)后剖面對比

圖5 單一縱向卡爾曼濾波25次(a)和50次(b)后剖面對比

2.2 實(shí)際資料測試

選用新疆吉木薩爾蘆凹陷實(shí)際地震資料(圖7)測試本文方法的實(shí)用性。針對該區(qū)二疊系蘆草溝組頁巖油儲(chǔ)層巖性復(fù)雜、縱向薄且目的層小斷裂發(fā)育導(dǎo)致水平井鉆遇率低的問題，重新采集了高密度高次覆蓋寬方位地震三維資料。

圖6 橫縱向同時(shí)濾波50次(a)和100次(b)后剖面對比

圖7 實(shí)際原始地震剖面(a)及本文方法濾波后剖面(b)

圖7a為OVT域疊前時(shí)間偏移成果數(shù)據(jù)?？梢娬w斷裂特征清楚，合成記錄與地震記錄對應(yīng)性較好，整體相關(guān)系數(shù)為70.5%(圖8a)。但為提高資料的薄層識別能力，保留較多高頻信息，高頻段的噪聲卻影響目的層小斷裂識別，由于目的層斷距大多小于10m，斷層的剖面地震相特征主要以同相軸彎曲、撓曲為主，因此曲率屬性在識別小斷裂方面具有優(yōu)勢。圖9a中目的層斷層多以延伸短的正斷層為主，但由于噪聲干擾，斷裂之間組合特征，尤其是靠資料邊界的斷層特征不清楚(圖中紅色圓圈內(nèi))。通過本文濾波方法后，資料信噪比及橫向連續(xù)性得到明顯改善，剖面上斷層特征清晰(圖7b)，濾波后地震記錄與原始記錄在縱、橫向上振幅強(qiáng)弱關(guān)系以及波組關(guān)系都保持較好的一致性，且合成記錄與原地震記錄對應(yīng)性仍然較好，整體相關(guān)系數(shù)保持在大約71.6%(圖8b)。

圖8 實(shí)際原始地震剖面(a)及本文方法濾波后剖面(b) 左：W1井旁道； 中：井?dāng)?shù)據(jù)合成道； 右左圖與中圖的相關(guān)

在濾波地震資料上提取目的層的曲率屬性，發(fā)現(xiàn)小斷層及其組合的特征清楚，小斷裂主要呈南北向條帶狀分布，走向以北東—南西向?yàn)橹?，同時(shí)在濾波后平面上可較清楚地展示一些與其正交的延續(xù)相對較長的斷裂，原始資料部分未識別出的斷層，在濾波后屬性平面上可清晰展示且去除了絕大部分的干擾信息(圖9b)； 同時(shí)從濾波后地震頻譜看，該方法與傳統(tǒng)的頻率濾波方法不同，只是消除高頻段的噪聲成分，保護(hù)了高頻有效信號(圖10)，還對10～30Hz頻帶缺失的部分頻率成分起到一定改善作用(圖10藍(lán)圈)。這樣就為工區(qū)目的層微構(gòu)造精細(xì)解釋及后續(xù)水平井位部署和壓裂方案設(shè)計(jì)提供了有利的保障。

圖10 本文方法濾波前(a)、后(b)頻譜對比圖

3 結(jié)論

地震資料經(jīng)過高頻補(bǔ)償處理后往往會(huì)表現(xiàn)出信噪比降低和橫向連續(xù)性變差等特點(diǎn)，本文鑒于卡爾曼濾波在恢復(fù)有效信號方面的優(yōu)勢，提出基于卡爾曼濾波原理對地震資料橫縱向同時(shí)濾波的方法，目的是提高地震資料信噪比和橫向連續(xù)性。通過方法研究并針對模型和實(shí)際資料測試，取得以下認(rèn)識：

(1)對地震資料做橫向卡爾曼濾波，所取得的信噪比和連續(xù)性的改善遠(yuǎn)優(yōu)于單一縱向?yàn)V波；

(2)經(jīng)過測試發(fā)現(xiàn)，基于卡爾曼濾波原理的地震資料橫、縱向以迭代循環(huán)濾波次數(shù)4∶1的形式進(jìn)行橫縱向組合濾波，在改善地震資料的信噪比和連續(xù)性的效果和效率兩方面都優(yōu)于單一橫向或縱向及二維中值濾波；

(3)實(shí)際資料測試表明，本文方法能顯著提高地震資料信噪比及同相軸的連續(xù)性，且不會(huì)造成高頻成分缺失，可有效確保構(gòu)造解釋精度。