潛江凹陷GH地區鹽巖強反射屏蔽下砂體識別方法研究

張水山,李銘華,羅 兵,張道洪,易小芹,張 磊,李 瓊,岳 林

(1.中國石化江漢油田分公司 勘探開發研究院,武漢 430223;2.成都理工大學 地球物理學院,成都 610059)

由于江漢盆地沉積在古氣候時期屬于強蒸發環境,潛江凹陷GH地區沉積鹽巖發育。潛江凹陷自西北向東南依次為砂泥巖區—過渡相區—鹽巖區,而GH地區屬于過渡相區,存在巖性組合復雜、砂巖線錯綜尖滅、厚鹽巖夾薄泥巖、膏泥巖中夾砂巖等情況。在地震剖面中鹽巖—泥巖波阻抗差異大,易形成強反射;砂—膏泥巖波阻抗差異小,形成弱反射,從而易被上覆強反射軸所掩蓋。這一現象在潛江凹陷地區廣泛存在,且該地區砂體具有數量多、面積小、分布散、儲層薄、砂體縱橫向變化快的特點,常規的地震資料很難滿足當前儲層解釋的要求,需要進行強反射剝離處理,以達到精細砂體解釋的要求。

對于強反射屏蔽弱有效信號這一問題,許多學者進行了研究。Mallat S.等[1-2]提出的基于匹配追蹤的稀疏分解方法,開創了信號自適應稀疏分解理論先河,隨后被地球物理學者引入到地震勘探中。Liu J.L.等[3-5]先后使用Ricker子波和Morlet小波構建原子庫,首先建立起了該方法與地震勘探的聯系,提出了動態最優搜索策略,并通過先驗參數約束搜索區間以減少搜索次數,提高了分解效率。張繁昌等[6-8]在Liu J.L.等[4]基礎上,采用動態匹配子波庫,實現了地震數據快速分解重構,并利用復數道分解對方法做了進一步改進,提高了分解速度與精度,對匹配子波進行施密特正交變換,降低了過完備匹配子波庫的冗余度。陳發宇等[9]采用地震有限頻帶約束初始頻率,提高了匹配速度。劉霞等[10]提出一種粒子群快速優化算法,用于快速搜索地震信號稀疏分解的最優匹配原子。王珺等[11]將遺傳算法和正交時頻原子相結合的匹配追蹤,實現了快速匹配追蹤,降低了分解度,提高了方法適用性。張雪冰[12]對地震信號的稀疏時頻分解做了系統的研究,給出了不同稀疏分解方法的適用條件,并通過引入衰減的Ricker子波,實現了較準確的Q值估計。鄧世廣等[13]運用并行算法,將匹配追蹤與偽Wigner-Ville分布結合,顯著提升了運算效率,得到了更高時頻分辨率的時頻譜。印興耀等[14]采用局部頻率對小波進行頻域先驗約束,加快了匹配的速度。郭志偉[15]研究了時頻分析在高精度地震資料處理重要環節中的許多應用。蔡涵鵬等[16]采用多參數構建Morlet匹配子波,快速精確地進行時頻譜分析,在烴類檢測方向取得了較好的效果。李海山等[17]采用低頻原子優化過完備字典,很好地去除了煤層強屏蔽。朱博華等[18]具體地闡述了匹配追蹤強反射去除時的原子構建參數的影響。張在金等[19]采用層位時間與解釋子波約束冗余字典,實現了快速去除強屏蔽,并提取到了儲層下方的低頻伴影異常。楊午陽等[20]以雷克子波作為原子,利用最大相關性估計分解方法提高了分解性能。何峰等[21]提出了基于測井聯合匹配追蹤的強反射去除方法,解決了原子子波難以估計的問題。常健強等[22]利用多道匹配追蹤子波分解與重構方法,對自來屯地區K21強反射層進行了剝離,提高了儲層預測精度。

匹配追蹤算法是一種信號稀疏分解算法,是壓縮感知領域中的信號重構算法,主要用于高維數據降維和信號稀疏表示。它可以通過迭代過程逐步獲取原始信號的稀疏表示,并且在迭代過程中可以進行實時處理,非常適合大規模數據處理和實時信號處理。本文通過研究地震匹配追蹤分解過程中的匹配速度敏感參數與匹配精度敏感參數,不依賴于冗余字典,而是對字典精度降采樣處理,進行多段強反射信號的匹配與重構,在此基礎上針對性地提出了基于雙精度降采樣字典的分步正交匹配追蹤強反射快速剝離改進算法。通過改進算法對實際資料進行強反射去除,有效去除了鹽巖強反射信號并保留了儲層弱反射信號,得到了更多構造信息。

1 改進正交匹配追蹤方法實現[12]

1.1 正交匹配追蹤原理

目前推廣于地震信號分解與處理領域,基于最小二乘問題的匹配追蹤算法是一種貪婪的迭代算法,其度量原子與信號匹配程度的數學工具為最小二乘解。

設原始信號為x,信號長度為L,則R(0)=x。以第k次迭代為例,其迭代過程相當于在冗余字典中尋找與殘差信號最匹配的原子W(k)。即:

(1)

(2)

其中:W(k)為Dpos(k),是該字典中的第pos(k)處的原子。pos(k)采用下式計算:

pos(k)=jmin((Wj-R(k-1))2),j=1,2,3,…,M

(3)

但在求解每個原子的振幅系數時采用下式計算:

(4)

上式用矩陣可以表示為:

(5)

其中Φ(k)為擴充的匹配原子庫,由k個原子組成,表示為:

Φ(k)=[Φ(k-1),W(k)]

(6)

A(k)為Φ(k)中的每個原子通過式求出的振幅向量,x為地震信號x(t)的向量形式,那么前k次迭代的子成分信號之和構成重構信號,為:

S(k)=Φ(k)A(k)

(7)

殘差信號為:

R(k)=x-Φ(k)A(k)

(8)

由于正交匹配追蹤要求所有原子Φ(k)均與剩余信號R(k)正交,即:

Φ(k)TR(k)=0

(9)

結合式(8)與式(9)可得,振幅系數向量的近似解為:

A(k)= (Φ(k)TΦ(k))-1Φ(k)Tx

(10)

值得注意的是第k-1次迭代與第k次的前k-1個原子相同,但振幅系數不相等。圖1為正交匹配追蹤的前兩次迭代示意圖。

Φ(k-1)(j)=Φ(k)(j),j=1,2,3,…,k-1

(11)

(12)

x=Φ(k)A(k)+R(k)

(13)

從圖1中可以看出,在第二次迭代搜索出最佳原子后,將第二次匹配到的原子擴充到匹配原子矩陣里,此時匹配原子矩陣包含前兩次迭代的最佳原子,通過匹配原子矩陣求取振幅向量,使x在W(1),W(2)所張開的平面Φ上投影最大,即S(2)最大,殘差R(2)最小,第二次迭代后的殘差R(2)⊥S(2)。

圖1 正交匹配追蹤算法前兩次迭代示意圖Fig.1 Schematic diagram of the first two iterations of the orthogonal matching tracking algorithm

1.2 方法改進

應用匹配追蹤分解地震信號時,首先要構造合適的過完備時頻字典。過完備時頻字典通過傳播時間、主頻、相位等具有地球物理意義的參數來定義,并通過各種參數的精細組合,來近似的獲得所有的地震子波原子。此時,字典中原子的數量遠大于原子信號的長度,稱為過完備字典或冗余字典。首先針對研究區GH的疊后地震數據體做了子波分析,如圖2所示。

圖2 子波分析Fig.2 Wavelet analysis

從地震資料中提取到如圖2所示子波,為一主頻31.5 Hz的零相位子波,長度約為60 ms。故本文運用Ricker子波來表征研究區地震子波。定義時頻字典需許多參數,如時移、頻率、相位、振幅等,這些參數既是構建稀疏字典的參數,當完全分解了地震信號后,又相當于反演了地下介質反射特征參數。理想的分解是所匹配的子波等于地下某界面所反射的子波;所匹配的原子參數等于地下某界面的反射參數。這種情況下,當匹配到強反射子波時,我們認為其就是地下強反射界面的反射子波,將該子成分信號從地震信號中減去,則完成了強反射軸的去除。那么,如何快速而準確的匹配成了關鍵。

常規的提升匹配精度的方法是擴大冗余字典的采樣率,這種做法是最直接的提升匹配精度的方法,但卻會極大的降低匹配速度。本文采用兩步字典法來既快速又精確地求解最小二乘值的最優解。標準的Ricker子波為:

r(t,f)=e-(πft)2(1-2(πft)2)

(14)

其中,t表示時移,f是子波主頻,子波原子采樣點長度L=20。首先用時移間隔為dt=0.001 s,頻率區間取30～40 Hz,間隔df=1 Hz,生成字典DⅠ。字典DⅠ按下式所示排列,其中每一個子波原子為一個列向量。

DⅠ=[r(t1,f1)r(t2,f1)r(t3,f1)…r(tLt,f1)
r(t1,f2)r(t2,f2)r(t3,f2)…r(tLt,fNf)]

(15)

其中:Lt=31為最大時移點數,Nf=11為最大頻率偏移點數,故字典DⅠ為一L×Lt×Nf的矩陣。矩陣寬M=Lt×Nf=341,似乎M不滿足遠大于L,這是因為本文采用了兩步字典法求解匹配原子,在后面構建DⅡ字典時加以解釋。

用字典DⅠ對正演剖面進行匹配追蹤,(正演剖面與實際剖面均為1ms采樣,故在與20點的子波原子匹配計算內積時會進行降采樣處理),最小二乘系數曲線如圖3所示。

圖3 最小二乘系數Fig.3 Least squares coefficient

根據字典DⅠ排列結構(公式15),可以看出最小二乘曲線對頻率并不敏感,但對時移非常敏感。且由于我們的字典原子主頻較低,可以認為在較小的區間內,最小曲線單調性變化不大。我們放大圖3中紅色箭頭處的最小二乘系數最小值附近曲線,如圖4所示。

從圖4中可以看出,使用字典DⅠ進行匹配,最小值在第9點,此時時移間隔過大,求解精度不夠。明顯第9點曲線左邊的梯度小于右邊的梯度。故可以通過此點確定精確值所在區間為(8,10),此為1步匹配法。進一步地,在區間(8,10),增加100倍時移采樣率,變為[8：0.01：10],得到臨時的局部精細字典DⅡ,故本文中實際字典長度M為101×Lt×Nf,滿足遠大于L的要求。且在逐點內積的過程中,避免了大量在非最小值區間的計算,使得運算過程中時間復雜度更低。其最小二乘曲線如圖 5中紅色曲線所示。

圖4 最小二乘方法最小值附近的曲線Fig.4 Curve near the minimum value of the least square method

圖5 改進前后最小二乘曲線對比Fig.5 Comparison of least squares curves before and after improvement

從圖5中可以看出,增加采樣率后,最小二乘曲線更加平滑,且出現了比點9處更小的值,其位置如圖中紅色點所示,其橫坐標8.75為所匹配出的時移解。所匹配出的原子波形對比如圖 6。

從圖6中可以看出,字典DⅡ匹配到的原子相對字典DⅠ匹配到的原子有一定偏置,匹配到的原子更加準確。

圖7為各道的偏置曲線,可以看出算法對各道的原子均有修正。兩步法求解原子相比直接將字典DⅠ擴大100采樣率更高效,計算速度更快,滿足實際應用的需求。但此時方法仍存在一定的缺點,雖然此時使用的低采樣率匹配強反射精度更高,但使用低采樣字典生成重構信號時會出現波形被截斷。如圖8所示的90-100道,子波波形被截斷,波形的旁瓣被“削去”,這會導致去除強反射后的剩余剖面不自然、強反射去除不干凈等問題,甚至在重建剖面產生假象(圖中箭頭處)。通過提高重構子波采樣點為41,解決了這一問題,如圖9。使用41點的原子子波去除之后的剖面上,剩余的信號同相軸更連續,剩余剖面上沒有出現因信號截斷而產生的假象,整體比20點短時窗原子子波去除效果更好。結果表明,通過在匹配時對原信號降采樣,可以減小高頻成分對匹配精度的影響,提高算法穩定性與匹配速度。在恢復信號時對子成分信號恢復采樣,能減小處理帶來的假象。

圖6 改進前后匹配原子波形對比Fig.6 Comparison of matched atomic waveforms before and after improvement

圖7 相較于改進前的偏置曲線Fig.7 Compared with the bias curve before the improvement

圖8 20點子波匹配的強反射Fig.8 Strong reflection of 20 dot wave matching

圖9 41點子波匹配的強反射Fig.9 Strong reflection of 41 dot wave matching

2 應用實例

圖10 Gu35井巖性解釋成果Fig.10 Lithology interpretation results of Well Gu35

圖11 實際地震剖面Fig.11 Actual seismic section

基于該方法對研究區進行強反射去除處理,如圖15所示,圖中紅色曲線為沿有效信號波峰位置。

圖12 原始剖面Fig.12 Original section

圖13 高分辨率處理后剖面(45 Hz)Fig.13 High-resolution profile (45 Hz)

圖14 鹽巖強反射剝離處理后剖面Fig.14 Section of salt rock after strong reflection stripping treatment

從圖16對比結果可以看出,處理前強反射屏蔽了砂巖弱信號,導致砂體平面特征不明顯。通過本文方法處理后,成功顯示出砂巖弱反射信號,提取的砂體平面特征明顯,與鉆井鉆遇砂巖位置對應良好,預測吻合率達到了84.60%。

沿解釋的砂體界面提取原始數據與高分辨率及強反射剝離后數據的均方根屬性,得到砂體平面特征如圖16。

圖16 處理前后均方根振幅屬性對比Fig.16 Comparison of RMS amplitude attributes before and after treatment

圖17 強反射剝離后砂體流體檢測平面屬性Fig.17 Fluid detection plane properties of sand body after strong reflection stripping

3 結 論

a.本文基于最小二乘優化的正交匹配追蹤算法進行研究,對GH地區鹽巖強反射屏蔽儲層有效信號的難題進行分析研究。利用地震資料,詳細分析了原子字典構建參數:采樣率、時移、頻率、相位、振幅對匹配速度與匹配精度的影響。提出基于雙精度降采樣字典的分步正交匹配追蹤強反射快速剝離改進算法,減小了在非最優解區間的匹配運算,兼顧了匹配速度與精度,提高了算法穩定性。