安圣培,劉 韜,胡天躍

(1.中國石油化工股份有限公司石油勘探開發(fā)研究院,北京100083;2.北京大學地球與空間科學學院,北京100871)

目前,中國西部地區(qū)已經成為油氣勘探的重點地區(qū)之一,該地區(qū)地表條件相對復雜,如山前帶,戈壁和沙漠,近地表非成巖低速介質區(qū)往往存在近地表高程、厚度和速度的空間變化,靜校正問題突出,從而影響疊加成像效果以及地表風化層速度異常的檢測。

折射類方法在解決靜校正問題中發(fā)揮了重要作用,如用于基準面靜校正的折射靜校正方法[1-2]和層析靜校正方法[3-4],以及基于統(tǒng)計相關實現的剩余靜校正方法[5-6]。這些方法的有效性都依賴于初至拾取的準確度,但復雜地表條件往往伴隨著強背景噪聲,降低了初至的信噪比,導致傳統(tǒng)的自動初至拾取的準確度下降,而手動初至拾取在面對大量地震數據時會產生極大的人力和時間成本。針對這一問題,有學者提出使用濾波類方法或地震干涉類方法提高初至信噪比,以保證自動拾取初至的準確性[7-9]。

在地震干涉類方法中,MIKESELL等[10]和ZHANG等[11]使用基于互相關的地震干涉法,計算兩點間的時延,并使用時延信息進行近地表速度建模,再基于速度模型計算剩余靜校正量。時延估計使用了互相關的方法,但互相關容易受高斯噪聲的影響,尤其是相干的高斯噪聲[12]。為解決這一問題,AN等[8]提出使用高階累積量替代互相關以提高地震干涉法的抗噪性,并由此發(fā)展了累積量相干積累法用于加強初至波信號。高斯噪聲的高階累積量為0,可有效提取高斯噪聲中的非高斯信號,更加符合實際地震數據的情況。高階累積量同時提供信號的幅度和相位信息,可用于描述信號系統(tǒng)的非線性特征?；谝陨咸攸c,高階累積量已經被廣泛應用于子波特征分析、噪聲衰減以及初至拾取等地震數據處理領域[13-15]。

在上述研究的基礎上,本文提出了基于高階累積量時間延遲估計的自動初至剩余靜校正方法,該方法使用高階累積量和相干疊加的方法,估計兩點間的初至波時延,再使用不同時延信息的組合得到炮檢點的剩余靜校正量。本文方法的優(yōu)勢在于:①互相關使用兩道數據進行時延估計,只能利用二維地震信息,而高階累積量使用多道數據,可以利用三維地震數據信息增強有效信號;②互相關是基于噪聲為高斯分布的假設,但實際地震數據很可能不滿足這一假設,導致道間時延估計的準確度不高,而高階累積量并不要求噪聲為高斯分布,更加符合實際地震數據的情況,對背景噪聲的壓制效果更優(yōu);③再結合地震干涉法中相干疊加的思路進一步提高信號的有效疊加次數,增強信噪比,保證本文方法獲得穩(wěn)定且可靠的道間時延估計結果;④本文方法采用完全數據驅動,相比基于初至拾取的剩余靜校正方法,在面對海量低信噪比地震數據時,可以有效降低初至拾取所需的高昂時間和人力成本。

1 方法原理

本文方法主要包括兩個部分:①使用高階累積量估計兩點間的初至波時延,并使用相干疊加的方法增強信噪比,保證時延估計的準確性;②使用時延信息的組合得到炮檢點的剩余靜校正量。

1.1 高階累積量時延估計

本文方法首先使用高階累積量估計兩點間的時延。二階統(tǒng)計工具,如互相關,常被用于時延估計,但互相關對高斯噪聲敏感,因而人們開始使用高階統(tǒng)計工具進行時延估計,如高階累積量。在實際問題中,加性噪聲往往具有高斯性,因此,使用高階累積量進行時延估計更加合理,因為高斯過程的高階累積量為0[16]。對零均值的平穩(wěn)隨機過程{x0(t),x1(t),…,xk-1(t)},其k階累積量由{x0(t),x1(t+τ1),…,xk-1(t+τk-1)}的聯合k階累積量表示如下[16]:

Ck(τ1,τ2,…,τk-1)=Cum{x0(t),x1(t+τ1),…,xk-1(t+τk-1)}

(1)

其中,Cum{·}是累積量算子,Ck(·)表示k階累積量函數,τ1,τ2,…,τk-1分別代表x0與x1,x2,…,xk-1之間的時延。常用的二階、三階和四階累積量函數的表達式如下[16]:

C2(τ)=E{x0(t),x1(t+τ)}

(2)

C3(τ1,τ2)=E{x0(t),x1(t+τ1),x2(t+τ2)}

(3)

C4(τ1,τ2,τ3)=E{x0(t),x1(t+τ1),x2(t+τ2),

x3(t+τ3)}-C2(τ1)C2(τ2-τ3)-C2(τ2)·

C2(τ3-τ1)-C2(τ3)C2(τ1-τ2)

(4)

其中,E{·}是期望算子。公式(2)表示協(xié)方差,在零均值的假設下等價于互相關。

基于高階累積量進行兩點間初至波時延估計的具體物理過程如下。如圖1所示,對于炮點S激發(fā),檢波點Rj接收的地震記錄Uj(t),假設它是由初至波Gj(t)和高斯隨機噪聲w(t)組成的,其表達式為:

Uj(t)=Gj(t)+w(t)

(5)

使用炮點S在檢波點Rj,Rj+1,…,Rj+n接收的n+1道地震記錄,計算其對應的n+1階累積量的公式如下:

Cj;j+1,j+2,…,j+n(τ1,τ2,…,τn)=Cum{Uj(t),

Uj+1(t+τ1),…,Uj+n(t+τn)}=Cum{Gj(t),

Gj+1(t+τ1),…,Gj+n(t+τn)}

(6)

其中,在階數≥3的高階累積量中,w(t)的響應為0。n+1階累積量包含n個時延τp(p=1,2,…,n),τp表示檢波點Rj與Rj+p之間的初至波時延。圖1中的V代表虛源位置,所謂虛源指的是假想的震源,并非真實存在的震源。S到Rj的旅行時與S到Rj+1的旅行時共用了S到V的旅行時,對應圖1中的灰色粗線,共用旅行時在估計時延的過程中被相互抵消,因此,τp只反映了兩個檢波點之間的時延,對應圖1中的紅線。高階累積量方法使用了n+1道信號進行時延估計,而互相關方法只使用兩道進行時延估計,因此,高階累積量方法的抗噪性更好。n+1階累積量Cj;j+1,j+2,…,j+n(τ1,τ2,…,τn)是n維函數,理論上可以直接求取n個時延,但隨著n的增大,存儲量和運算量都會顯著增加。因此,本文方法每次只計算一個時延,實現方式如下:

圖1 使用高階累積量估計兩點間時延示意

(7)

其中,τ是第j+1個檢波點與第j個檢波點之間的初至波時延。公式(7)只以兩道地震記錄為輸入,理論上仍可以得到n階累積量,但由于Gj(t)與自身的時延為0,因此,不再需要計算累積量中的除τ以外的其它時延,大大節(jié)省了運算量,并且仍然利用了高階累積量的抗噪性優(yōu)勢。

1.2 相干疊加

為了保證兩點間時延估計的準確性,本文進一步使用相干疊加方法,增強兩點間初至波信噪比,具體過程如下。如圖2所示,對炮點S1在檢波點R1和R2的初至波記錄進行高階累積量計算,共用的S1到虛源V的旅行時在時延估計過程中被抵消,因此,得到了R1和R2間的初至波響應(圖2b中紅線所示),圖中實線表示旅行時為正,虛線表示旅行時為負。使用其它炮點Si(i=2,3,…,N)重復上述過程,炮點Si到V的旅行時仍共用,在時延估計時被抵消,因此,使用不同炮點得到的高階累積量結果中包含了相同的R1到R2的道間時差(圖2b)。對圖2b中的所有結果進行疊加,實現了初至波響應的同相位疊加,而非相干的隨機噪聲由于具有不同的時延,不能實現同相位疊加,疊加能量將小于初至波。至此,實現了使用相干疊加增強兩點間初至波響應的信噪比。

圖2 相干疊加示意a 對炮點在檢波點R1和R2的初至波記錄進行高階累積量計算; b 不同炮點得到的高階累積量結果; c 圖2b的疊加結果

1.3 時延組合求取剩余靜校正量

(8)

(9)

(10)

(11)

圖3 使用時延信息求取剩余靜校正量示意a 使用檢波點左側炮點的地震記錄進行檢波點Rj和Rj+1間的時延估計; b 使用檢波點右側炮點的地震記錄進行檢波點Rj+1和Rj間的時延估計; c 使用a和b中得到的時延信息進行相減得到相鄰檢波點間的相對剩余靜校正量

ΔTj,j+1和ΔTj+1,j都包含了傳播至Rj和Rj+1的兩組上行旅行時和沿高速層頂傳播的旅行時。假設地表高程和近地表低速層厚度引起的長波長靜校正量已經被消除,則ΔTj,j+1和ΔTj+1,j包含的上、下行旅行時相同;同時高速層頂界面的起伏在長波長靜校正中被修正,ΔTj,j+1和ΔTj+1,j包含的沿高速層頂的傳播距離也近似相等,由于Rj和Rj+1為相鄰檢波點,在相對小的空間范圍內,可以認為高速層頂的速度變化不大,因此ΔTj,j+1和ΔTj+1,j包含的沿高速層頂傳播的旅行時也近似相等,ΔTj,j+1和ΔTj+1,j近似相等,ZHANG等[11]給出過類似的證明。因此,公式(10)和公式(11)相減得到:

(12)

1.4 方法流程

本文方法的實現流程如圖4所示。如果數據存在長波長靜校正問題,需要先進行高程靜校正或野外靜校正;針對數據可能存在的異常高頻噪聲,使用高通濾波進行噪聲壓制;使用高階累積量和相干疊加方法得到兩點間的初至時延;最后,使用時延信息的組合得到炮檢點的剩余靜校正量。

圖4 本文方法流程

2 合成數據實例

2.1 高階累積量時延估計實例

本節(jié)通過合成數據實例,對比了高階累積量和互相關方法應用于時延估計的抗噪性。圖5a和圖5b展示了僅包含一組折射波信號的兩個單炮記錄,折射波視速度為2500m/s,檢波點間距為15m,兩個炮點相距60m,雷克子波主頻為60Hz。圖5c和圖5d分別為對圖5a和圖5b加入相干高斯噪聲的單炮記錄,其信噪比為-9dB。圖5a和圖5b的互相關結果如圖5e所示。圖5f、圖5g和圖5h分別為對含噪數據(圖5c和圖5d)使用互相關、三階累積量和四階累積量得到的結果。由于互相關對相干噪聲敏感,圖5f的互相關結果中除了包含折射波的響應,還混入了相干噪聲的響應,信噪比為3dB,圖5g和圖5h中的相干噪聲被有效壓制,殘余噪聲也不具有相同的時延,因此無法實現相干疊加,其信噪比分別達到了9dB和15dB。圖5i、圖5j、圖5k和圖5l分別為對圖5e、圖5f、圖5g和圖5h的結果沿檢波點疊加得到的疊加單道。圖5e使用不同檢波點得到的折射波響應具有相同的時延,因此在疊加過程中實現了相干疊加,疊加結果(圖5i)中的峰值位置與理論時延(紅色虛線標示)一致。圖5j中相干噪聲響應的峰值大于折射波響應的峰值,因此難以準確得到折射波響應的時延估計。圖5k和圖5l中包含的折射波響應的峰值位置與理論時延一致。該實例證明了三階和四階的高階累積量方法壓制相干白噪聲的效果優(yōu)于互相關方法,同時,高階累積量的階數越高,得到的折射波響應的信噪比也越高。

圖5 合成數據的測試結果a,b 僅包含一組折射波信號的單炮記錄; c,d 加入相干高斯噪聲的單炮記錄; e 圖5a和圖5b的互相關結果; f 圖5c和圖5d的互相關結果; g 圖5c和圖5d的三階累積量結果; h 圖5c和圖5d的四階累積量結果; i 對圖5e疊加得到的疊加單道; j 對圖5f疊加得到的疊加單道; k 對圖5g 疊加得到的疊加單道; l 對圖5h疊加得到的疊加單道

2.2 求取剩余靜校正量實例

采用聲波有限差分正演得到合成數據,所用模型參數見表1,炮距為40m,道距為10m,炮點數100,檢波點數400。剩余靜校正量由隨機序列產生,滿足均值為0和標準差為8ms的正太分布,正演得到的單炮記錄如圖6a所示,由于剩余靜校正量的存在,圖6a中同相軸的連續(xù)性較差。將采用本文方法得到的剩余靜校正量應用到單炮記錄中,得到的結果如圖6b 所示。由圖6b可見,單炮記錄同相軸的連續(xù)性顯著提高。圖7給出了對不含噪合成數據采用本文方法得到的剩余靜校正量。其中,灰線表示真實靜校正量,紅線表示計算得到的靜校正量。由圖7可見,本文方法計算得到的剩余靜校正量與真實值完全吻合。

圖6 對不含噪合成數據(a)應用本文方法進行剩余靜校正處理后的結果(b)

圖7 對不含噪合成數據應用本文方法得到的剩余靜校正量

表1 模型參數

對圖6a數據加入相干高斯噪聲,得到信噪比為-7dB的含噪數據(圖8a)。將采用本文方法得到的剩余靜校正量應用到含噪單炮記錄中,得到如圖8b所示的結果。由圖8b可見,反射波同相軸連續(xù)性明顯提高。圖8c和圖8d分別是圖8a和圖8b方框標示區(qū)域的局部放大圖。圖9展示了對含噪合成數據使用互相關和本文方法得到的剩余靜校正量。由圖9可以看出,本文方法得到的剩余靜校正量與真實值基本一致,互相關方法受相干噪聲影響,其結果偏離真實值較多。該實例驗證了本文方法求取剩余靜校正量的有效性,同時,相比互相關方法,本文方法不易受相干噪聲影響,抗噪性更好。

圖8 對含噪合成數據應用本文方法進行剩余靜校正處理后的結果a 未做剩余靜校正的單炮記錄; b 采用本文方法進行剩余靜校正處理后的單炮記錄; c 圖8a方框區(qū)域局部放大結果; d 圖8b方框區(qū)域局部放大結果

圖9 對含噪合成數據使用互相關和本文方法得到的剩余靜校正量

3 實際資料應用

選取中國西部某地區(qū)的實際地震資料驗證本文方法的有效性。該工區(qū)地表存在較強的近地表速度異常,剩余靜校正問題突出。對比了本文方法和某商業(yè)軟件中的基于初至拾取的初至波剩余靜校正模塊的結果,該模塊采用的是廣義互換的方法?；诔踔潦叭〉氖Ｓ囔o校正方法要求拾取全部道集的初至到時,但在遠偏移距處,受該工區(qū)復雜近地表影響,初至信噪比極低,不得不進行耗時耗力的手動初至拾取,而本文方法可以避免初至拾取的工作量。圖10a和圖10b對比了本文方法和基于手動初至拾取方法得到的檢波點和炮點的剩余靜校正量,結果基本吻合,僅在檢波點號為900附近位置存在微弱差異。將剩余靜校正量應用到道集數據,選取剩余靜校正問題比較突出的原始單炮(圖11a)進行比較,該單炮初至波同相軸明顯不連續(xù),圖11b和圖11c分別是應用本文方法和基于初至拾取方法的結果,兩種方法都明顯改善了初至波同相軸的連續(xù)性,圖11d至圖11f分別是圖11a至圖11c黑色方框處的局部放大圖,反射波同相軸的連續(xù)性也明顯改善。圖12對比了本文方法和基于手動初至拾取方法得到的疊加剖面。圖12a為未進行剩余靜校正的疊加剖面。由圖12a可見,該疊加剖面同相軸連續(xù)性較差。采用本文方法得到的疊加剖面(圖12b)同相軸的連續(xù)性顯著改善,并且與基于手動初至拾取的剩余靜校正方法相比,本文方法避免了手動初至拾取的繁重工作量。

圖10 本文方法和基于手動初至拾取方法應用于實際數據得到的剩余靜校正量對比a 檢波點剩余靜校正量對比; b 炮點剩余靜校正量對比

圖11 本文方法和基于手動初至拾取方法進行剩余靜校正的單炮記錄效果對比a 未進行剩余靜校正的單炮記錄; b 采用本文方法進行剩余靜校正后的單炮記錄; c 基于手動初至拾取方法進行剩余靜校正后的單炮記錄; d 圖11a黑色方框區(qū)域的局部放大結果; e 圖11b黑色方框區(qū)域的局部放大結果; f 圖11c黑色方框區(qū)域的局部放大結果

圖12 本文方法和基于手動初至拾取方法進行剩余靜校正后的疊加剖面效果對比a 未進行剩余靜校正的疊加剖面; b 采用本文方法進行剩余靜校正后的疊加剖面; c 基于手動初至拾取方法進行剩余靜校正后的疊加剖面

4 結論與認識

本文發(fā)展了基于高階累積量時間延遲估計的自動剩余靜校正方法,采用高階累積量得到兩點間的初至波時延;使用相干疊加提高信噪比,保證時延估計的準確性;再使用時延信息的組合得到炮檢點的剩余靜校正量。與互相關類方法相比,使用高階累積量進行時延估計可有效利用三維地震數據信息增強有效信號,高階累積量不再要求噪聲為高斯分布,與實際地震資料的情況更加吻合;采用相干疊加方法提高信號有效疊加次數和信噪比,保證本文方法獲得穩(wěn)定且可靠的道間時延估計結果;本文方法采用完全數據驅動,與基于初至拾取的剩余靜校正方法相比,可有效避免初至拾取所需的高昂時間和人力成本。模型試驗和實際資料的應用結果表明,本文方法可以明顯改善疊加剖面中反射波同相軸的連續(xù)性,為后續(xù)地震數據處理提供高品質資料。

本文方法涉及到高階累積量階數的選取,選取的階數越高,噪聲壓制效果越好,但運算量也同步增加,本文選取了常用的四階累積量,在具體應用于實際地震資料時,可事先通過參數測試確定適合的累積量階數,實現高效穩(wěn)定的剩余靜校正量求取。另外,本文方法要求觀測系統(tǒng)規(guī)則化,以實現兩點間初至波響應的相干疊加,針對不規(guī)則觀測系統(tǒng),可事先使用空間插值技術重構得到規(guī)則化道集。本文方法實現了自動剩余靜校正,未來可結合人工智能類方法實現流程參數的自動優(yōu)化選取,進一步提高自動化程度以減少人工成本,更加高效地應用于實際工區(qū)的地震資料處理。本文方法可進一步與深度學習結合,改進傳統(tǒng)的基于特征值的時延估計方法,提高道間時延估計的準確率和對更加復雜噪聲環(huán)境的適用性,并同時降低特征提取和篩選的人工成本,提高計算效率。