張亞南，戴永壽，陳 健，魏玉琴，丁進(jìn)杰，張漫漫

中國石油大學(xué)（華東）信息與控制工程學(xué)院，山東東營 257061

1 引 言

地震子波提取是地震資料處理和地震解釋的重要組成部分，但子波提取過程中子波相位往往估計不準(zhǔn)確，嚴(yán)重影響地震反演處理和地震解釋.針對地震子波相位提取的問題，相關(guān)領(lǐng)域的學(xué)者做了大量研究.唐斌等［1］將子波相位分成最小相位和最大相位兩個部分，通過高階累積量方法將二者恢復(fù)，從而提取子波相位.伊振林等［2］利用延遲因子構(gòu)建混合相位濾波器進(jìn)行混合相位反褶積，通過最小熵準(zhǔn)則和Lp范數(shù)約束反褶積結(jié)果，進(jìn)而獲取最佳濾波器和子波相位.崔慶輝等［3］通過地震道自相關(guān)提取子波振幅，利用雙譜法提取子波相位，最終重構(gòu)出混合相位的地震子波.但在實際應(yīng)用過程中，以上方法很難得到精確結(jié)果，以至于造成反演結(jié)果的失真.為改善反演結(jié)果，國內(nèi)外科研人員提出了多種相位校正方法.Levy等［4－5］假設(shè)提取的子波與真實子波的差異為不依賴于頻率的常數(shù)，采用常相位校正方法對反演結(jié)果進(jìn)行校正，以消除子波相位估計不準(zhǔn)的影響.Wang［6－7］提出了一種基于波場下延理論的常Q反濾波方法，通過兩項指數(shù)對大地Q濾波的相位和振幅進(jìn)行校正，并給出其穩(wěn)定性的證明.到目前為止，人們提出的各種相位校正方法并不能完全有效地消除反演結(jié)果中子波相位估計不準(zhǔn)的影響，因此子波相位對反演影響的研究將成為今后的攻關(guān)課題和研究熱點.Yuan等［8］針對滑動平均（Moving average，MA）模型描述下地震子波相位估計不準(zhǔn)對反演結(jié)果影響的問題進(jìn)行了討論，并給出了多種方法評價反演結(jié)果.但由于其模型所含參數(shù)較多，造成多種相位可能性；且所用子波的長度過短，與實際子波相比具有一定差異［9］.因此，本文在ARMA模型描述子波的基礎(chǔ)上，忽略噪聲等因素的影響，重點探討子波相位對反射系數(shù)序列反演結(jié)果的影響，并通過理論分析和仿真實驗研究其影響規(guī)律.

本文首先采用在z域?qū)ΨQ映射ARMA模型零極點的方式構(gòu)造出一系列相同振幅譜、不同相位譜的子波，進(jìn)而對人工合成地震記錄進(jìn)行反射系數(shù)序列反演并對結(jié)果進(jìn)行分析，通過仿真實驗和實際地震資料處理來驗證理論分析結(jié)果.本文采用兩種方法對反演結(jié)果進(jìn)行評價，以確定出真實或準(zhǔn)確的反射系數(shù)序列.

2 構(gòu)造不同相位譜的地震子波

戴永壽等［10］將ARMA模型引入到地震子波的估計當(dāng)中，ARMA模型相對于MA模型具有參數(shù)吝嗇的性質(zhì)，用較少的參數(shù)即可描述一個精確的子波.Robinson褶積模型［11］可用ARMA模型表示為

其中，x（n）為地震記錄；r（n）為反射系數(shù)序列；子波AR部分的參數(shù)ai和MA部分的參數(shù)bk均為實數(shù)；p為AR階數(shù)，q為MA階數(shù).

假設(shè)ARMA描述下地震子波z域變換的系統(tǒng)函數(shù)W（z）無零極點對消，則

其中w（m）為地震子波的時間域序列，W（z）為其z域表示.由于ai和bk均為實系數(shù)，因此系統(tǒng)零點ck與極點di均為實根或者共軛復(fù)根［12］，A為增益常數(shù).當(dāng)z＝ejω時，式（2）中子波的z域表示轉(zhuǎn)換為頻率響應(yīng)，即

3 反射系數(shù)序列反演

在已經(jīng)構(gòu)造出一系列相同振幅譜、不同相位譜子波的基礎(chǔ)上研究子波相位不準(zhǔn)對反射系數(shù)反演的影響，首先需要找尋一個合適的反演方法.現(xiàn)有的反演方法已經(jīng)十分成熟，如最小二乘反褶積、最小熵反褶積［13］、稀疏脈沖反褶積［14］、多分辨率地震信號反褶積［15］、基追蹤反演［16］等.由于ARMA模型描述下的子波在時域表現(xiàn)為無限長脈沖響應(yīng)序列，時域反褶積處理不可避免地要對子波進(jìn)行截斷.為減少截斷誤差對反演過程影響，我們選擇譜除法在頻域進(jìn)行反褶積處理，處理結(jié)束后將頻域結(jié)果映射到時域進(jìn)行比較.

將Robinson褶積模型經(jīng)傅里葉變換可得

其中X（ejω）、W（ejω）和R（ejω）分別為地震記錄、地震子波和反射系數(shù)序列的頻域表示，則反射系數(shù)序列的估計可表示為

由式（14）可知，子波相位估計不準(zhǔn)的反褶積結(jié)果中會殘留一個純相位濾波器，即反射系數(shù)序列反演的結(jié)果為真實反射系數(shù)序列與一個純相位濾波器的褶積，且此純相位濾波器的相位譜為真實子波和估計子波相位譜之差.估計的反射系數(shù)序列應(yīng)當(dāng)和真實反射系數(shù)序列具有相同的振幅譜和能量，僅相位譜存在差異.由于純相位濾波器的影響，估計的反射系數(shù)序列出現(xiàn)了相位偏移，不符合真實反射系數(shù)序列的分布規(guī)律，有可能造成假象的出現(xiàn)和分辨率的降低［8］.

4 數(shù)據(jù)仿真驗證

為了驗證ARMA模型描述下子波相位估計不準(zhǔn)對反射系數(shù)反演結(jié)果的影響，本文使用具有四個反射系數(shù)脈沖的稀疏模型進(jìn)行實驗，其長度為1000ms，采樣頻率為1kHz，如圖1所示.其中稀疏反射系數(shù)脈沖的位置和大小分別為（200ms，1.0）、（300ms，－0.8）、（600ms，0.8）和（800ms，－1.0）.

本實驗中采用因果混合相位模型的真實子波，其ARMA描述下的差分形式為

分別求取式（16）中AR和MA部分的根可得真實子波的極點：α1＝0.8643＋0.1666j，α2＝0.8643－0.1666j，α3＝0.3107＋0.4177j，α4＝0.3107－0.4177j；真實子波的零點：β1＝1.2015，β2＝－0.2008＋0.8013j，β3＝－0.2008－0.8013j.子波的時域和頻域表示如圖2所示，其中時域波形進(jìn)行了截斷處理.文中相位譜中令相位角φ∈（－π，π］.

將真實子波與反射系數(shù)序列作褶積處理合成地震記錄，合成的地震記錄如圖3所示.由于所采用的子波為ARMA模型描述，在時間域具有很強的延續(xù)性，所以圖3a中會出現(xiàn)較長的拖尾現(xiàn)象.圖3b與圖1相比，合成地震記錄的能量大部分集中在了低頻區(qū)域.

按照第2節(jié)中所介紹的方法構(gòu)造出一系列相同振幅譜、不同相位譜的子波，其z域中零極點的分布如圖4所示.其中“○”代表零點，“×”代表極點.第一行（4a—4d）、第 二 行 和 第 三 行（4e—4l）、第 四 行（4m—4p）分別為因果、混合因果、反因果子波的零極點分布，第一列、第二列和第三列、第四列分別為最小相位、混合相位、最大相位子波的零極點分布.圖4b為真實的因果混合相位子波在z域的零極點分布.

圖5所示為相同振幅譜不同相位譜的地震子波的時域序列和相位譜.由圖5可知，時間域各道子波具有相同振幅譜，由于相位譜的不同，在時域中表現(xiàn)出了不同的性質(zhì).圖5b中的每道波形為圖5a中相應(yīng)地震子波的相位譜，相位角φ∈（－π，π］.由圖5a可以看出，1—4道子波為因果子波，其時間域響應(yīng)均在零點之后；5—12道子波為混合因果子波，在零點前后均有時間域響應(yīng)；13—16道子波為反因果子波，其時間域響應(yīng)均在零點之前，其中第2道為真實子波.在實際的地震勘探過程中，子波應(yīng)當(dāng)是因果混合相位的，但由于檢波器等原因會引入少量的反因果成分，故在實際處理過程中不能單純地將子波看作僅含因果成分.由圖4和圖5a可以看出，在因果穩(wěn)定系統(tǒng)的情況下，零點均在單位圓內(nèi)為最小相位，內(nèi)外均有為混合相位，均在單位圓外為最大相位.當(dāng)系統(tǒng)為反因果時則與因果系統(tǒng)完全相反，混合因果時需要分情況討論，在此不再贅述.

應(yīng)用構(gòu)造出的相同振幅譜不同相位譜子波對合成的地震記錄進(jìn)行反射系數(shù)序列反演，即反褶積處理.16個不同相位譜子波反褶積處理后的結(jié)果如圖6所示，其中第2道為真實子波反演結(jié)果.

圖3 合成的地震記錄（a）時域波形圖和（b）振幅譜Fig.3 The time domain waveform（a）and amplitude spectrum（b）of synthetic seismogram

由圖6可以看出，僅真實子波能夠完全恢復(fù)真實的反射系數(shù)序列，其余相位不準(zhǔn)的子波反演結(jié)果均在時域呈現(xiàn)出反射系數(shù)拓展的形式，不滿足真實反射系數(shù)序列的稀疏性質(zhì).如果本實驗結(jié)果滿足第3節(jié)中的結(jié)論——在子波相位估計不準(zhǔn)的情況下反射系數(shù)序列反演的結(jié)果為反射系數(shù)序列與一個純相位濾波器的褶積，則對所有估計子波進(jìn)行反射系數(shù)序列反演的結(jié)果應(yīng)當(dāng)具有相同的振幅譜，即相同的能量.根據(jù)能量在時域和頻域等價原則，計算圖6中所有反射系數(shù)序列的時域和頻域能量，其結(jié)果如表1所示，第2道為真實子波反演的能量值.

由表1可知，不同相位子波的反演結(jié)果所含能量略有差異，其根本原因在于采用的ARMA模型為無限長脈沖響應(yīng)系統(tǒng)，由于實驗用反射系數(shù)序列的長度為1000ms，在實際處理過程中進(jìn)行了截斷處理，從而產(chǎn)生截斷誤差；數(shù)值計算存在精度誤差.但估計子波與真實子波反演結(jié)果的能量值最大差異不超過0.031%，屬可接受范圍，可認(rèn)為估計子波的反演結(jié)果具有相同的能量，可初步驗證文中第3節(jié)中的結(jié)論.

為進(jìn)一步驗證第3節(jié)中的結(jié)論，我們將求取出純相位濾波器.由式（14）可知，純相位濾波器P（ejω）＝ejφW（ω）－jφ~W（ω）由以下公式得到

若上式中分母中所估計的子波頻率響應(yīng)的某些頻率為零值時，將其用極小實數(shù)值代替以避免計算錯誤.

純相位濾波器應(yīng)當(dāng)具有相同的單位振幅譜，并且在時域和頻域應(yīng)當(dāng)具有單位能量.圖7所示為由式（17）所計算出的純相位濾波器，每道中各個頻率的幅值均為1，圖7a為其時間域波形，圖7b為其相位譜.對圖7中所示的16個純相位濾波器進(jìn)行時域和頻域能量的求取，結(jié)果如表2所示.

表1 不同子波反演結(jié)果的能量值Table 1 The energy inversion results with different wavelets

表2 不同純相位濾波器的能量值Table 2 The energy of different pure－phase filters

由表2可以看出，數(shù)值計算誤差和截斷誤差在0.03%以內(nèi)，屬可接受范圍，可認(rèn)為所求得的結(jié)果在時域和頻域均具有相同的單位能量.由圖7和表2可以驗證由式（17）計算的結(jié)果為純相位濾波器.

令真實反射系數(shù)序列激勵各純相位濾波器，所得結(jié)果如圖8所示.由圖8所得結(jié)果與圖6反射系數(shù)序列反演結(jié)果比較，二者十分相似，為進(jìn)一步驗證本文結(jié)果，將圖6與圖8數(shù)據(jù)根據(jù)相同道號求相似度，其結(jié)果如表3所示.

由表3可以看出，反射系數(shù)序列反演的結(jié)果與真實反射系數(shù)序列激勵純相位濾波器輸出結(jié)果的相似度非常高，其中的誤差可以認(rèn)為是數(shù)據(jù)截斷誤差和數(shù)值運算精度誤差.由此可以認(rèn)為在子波相位估計不準(zhǔn)的情況下反射系數(shù)序列反演的結(jié)果為反射系數(shù)序列與一個純相位濾波器的褶積.

為驗證純相位濾波器的相位譜為真實子波與估計子波的相位譜之差，將圖5中真實子波與各道子波的相位譜做相減運算，其結(jié)果如圖9所示.

由圖9和圖7b對比可看出二者十分相似，為進(jìn)一步驗證二者的相似性，將圖9和圖7b中的數(shù)據(jù)按相同道求取相似度，其結(jié)果如表4所示.由表4可以看出純相位濾波器的相位譜與真實子波－估計子波的相位譜之差相一致，可以認(rèn)為，純相位濾波器的相位譜為原始子波與估計子波的相位譜之差.

以上實驗結(jié)果表明，子波相位估計不準(zhǔn)對反射系數(shù)序列反演具有影響，反演結(jié)果為真實反射系數(shù)序列與一個純相位濾波器的褶積，且此純相位濾波器的相位譜為真實子波和估計子波的相位譜之差.

圖8 反射系數(shù)序列激勵不同純相位濾波器的輸出結(jié)果Fig.8 The output of reflection coefficient sequences incentive different pure－phase filters

表3 反射系數(shù)序列反演結(jié)果與純相位濾波結(jié)果的相似度Table 3 Similarity between inversion results of reflection coefficient sequences and pure－phase filter results

表4 純相位濾波器的相位譜與原始－估計子波相位譜之差的相似度Table 4 Phase spectrums similarity between pure－phase filters and the difference of real－estimated wavelets

圖9 真實子波與估計子波相位譜之差Fig.9 The phase spectrum difference between real wavelet and estimated wavelets

5 反射系數(shù)序列反演結(jié)果的評價

在實際地震信號處理過程中，反射系數(shù)序列和真實子波均為未知，真實子波的振幅譜較為容易求得，但其相位譜往往估計不準(zhǔn)，造成反射系數(shù)序列反演結(jié)果的相位偏移.假設(shè)真實反射系數(shù)序列是稀疏脈沖序列，時域能量集中在數(shù)個脈沖內(nèi)；子波相位不準(zhǔn)的反演結(jié)果雖然具有相同的能量，但由于改變了其相位譜，時域能量則會分散.因此應(yīng)用如下評價方法對反演結(jié)果進(jìn)行評價，以確定出真實或準(zhǔn)確的反射系數(shù)序列.

圖10 反射系數(shù)序列的（a）變分評價結(jié)果和（b）豐度評價結(jié)果Fig.10 Evaluation results of reflection coefficient sequences with Variation（a）and Kurtosis（b）

其中式（18）和式（19）分別稱為變分［17］和豐度［18］方法，應(yīng)用二者對本文的反射系數(shù)序列反演結(jié)果進(jìn)行評價，得到如圖10所示結(jié)果.

由圖10可以看出，變分方法在第2道取最小值，豐度方法在第2道取最大值，表明二者均可以確定出準(zhǔn)確的反射系數(shù)序列.由于真實反射系數(shù)序列包含在反演結(jié)果中，則確定出的準(zhǔn)確反射系數(shù)序列即為真實反射系數(shù)序列.

為驗證以上方法應(yīng)用的普遍性，假設(shè)真實子波為因果最小相位子波（第1道）、混合因果混合相位子波（第11道）、反因果最大相位子波（第16道），其反演結(jié)果與評價結(jié)果如圖11所示.

圖11表明，在真實子波為任意相位和任意因果性的條件下，均可以通過變分和豐度方法對反演結(jié)果進(jìn)行評價，從而確定出真實或準(zhǔn)確的反射系數(shù)序列.同樣，采用驗證時域能量集中性的其它評價方法也可得到相應(yīng)的結(jié)果，在此不作一一表述.

6 實際地震資料處理

為驗證子波相位不準(zhǔn)對反射系數(shù)序列反演的影響規(guī)律，我們將對實際地震資料和測井所得反射系數(shù)序列進(jìn)行處理.圖12b為勝利某區(qū)塊的地震記錄，地震記錄的采樣速率為2ms，井位于第289道，圖12a為根據(jù)該井的測井?dāng)?shù)據(jù)所轉(zhuǎn)換的時域反射系數(shù)序列.

對圖12所示的實際地震資料進(jìn)行地震子波提取.采用張廣智等［19］提出的井旁道地震子波精細(xì)提取方法結(jié)合測井資料提取井旁道子波，采用原始的統(tǒng)計性子波提取方法提取最小相位的地震子波，二者具有相同的振幅譜.子波提取結(jié)果如圖13所示.

結(jié)合井旁道測井資料所提取的地震子波為混合因果混合相位的，而傳統(tǒng)方法所提取的地震子波為因果最小相位的，如圖13所示.二者具有相同的振幅譜，僅相位譜存在差異，為驗證本文所提理論的有效性，對二者進(jìn)行相位譜差值化處理，以構(gòu)造純相位濾波器.構(gòu)造出的純相位濾波器如圖14所示.

應(yīng)用提取出的最小相位地震子波對實際地震記錄作反射系數(shù)反演處理，將反演結(jié)果與實際反射系數(shù)序列激勵純相位濾波器的輸出進(jìn)行比較，若二者相一致則可證明本文所提理論的有效性.測井所得實際反射系數(shù)序列與反演處理結(jié)果如圖15所示.

對比圖15a和圖15b可知，當(dāng)子波相位估計不準(zhǔn)時，得到的反演結(jié)果與真實的反射系數(shù)序列相差甚遠(yuǎn)，不能夠正確反映地層的實際信息；由于提取的子波振幅譜是準(zhǔn)確的，故地震記錄中的低頻成分經(jīng)反演后壓制的較為理想.對比圖15b和圖15c可以看出，反射系數(shù)序列激勵圖14中的純相位濾波器的輸出與最小相位子波反演的結(jié)果較為相似.對最小相位子波反演結(jié)果和反射系數(shù)序列純相位濾波結(jié)果計算相似度，其相似度為73.74%，二者具有較好的相似性.考慮到測井資料和地震道中的噪聲、深時轉(zhuǎn)換的誤差等原因，給子波反演結(jié)果和反射系數(shù)序列求取結(jié)果帶來一定誤差，二者相似度在可接受范圍內(nèi)，從而表明實際數(shù)據(jù)處理結(jié)果可以驗證文中結(jié)純相位濾波器的褶積，此濾波器的相位譜為真實子波和估計子波相位譜之差.由于純相位濾波器的作用，反褶積結(jié)果破壞了真實反射系數(shù)序列的稀疏性（時域能量集中性）.通過兩種方法對反演結(jié)果進(jìn)行了評價，評價結(jié)果表明二者均可確定出真實或準(zhǔn)確的反射系數(shù)序列.本文針對子波相位估計不準(zhǔn)對反射系數(shù)序列反演的影響進(jìn)行了理論分析、實驗仿真和實際地震資料處理，仿真和實際處理結(jié)果可以驗證理論結(jié)果，為后續(xù)進(jìn)一步提高反射系數(shù)序列反演結(jié)果精度指明了研究方向.

表5 實際資料處理的變分和豐度評價結(jié)果Table 5 Variation and Kurtosis evaluation results of actual seismic data

致 謝 感謝中國石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院印興耀教授和張廣智教授的大力支持，感謝中國石油大學(xué)（北京）袁三一老師的指導(dǎo)與建議.

（

）

［1］ 唐斌，尹成.基于高階統(tǒng)計的非最小相位地震子波恢復(fù).地球物理學(xué)報，2001，44（3）：404－410.

Tang B，Yin C.Non－minimum phase seismic wavelet reconstruction based on higher order statistics.ChineseJ.Geophys.（in Chinese），2001，44（3）：404－410.

［2］ 伊振林，王潤秋.一種新的混合相位反褶積方法.石油地球物理勘探，2006，41（3）：266－270.

Yi Z L，Wang R Q.A new mixed phase deconvolution method.OilGeophysicalProspecting（in Chinese），2006，41（3）：266－270.

［3］ 崔慶輝，芮擁軍，尚新民等.混合相位地震子波提取及應(yīng)用.石油物探，2011，50（5）：481－486.

Cui Q H，Rui Y J，Shang X M，et al.Mixed－phase seismic wavelet extraction and its application.GeophysicalProspecting forPetroleum（in Chinese），2011，50（5）：481－486.

［4］ Levy S，Oldenburg D W.The deconvolution of phase－shifted wavelets.Geophysics，1982，47（9）：1285－1294.

［5］ Levy S，Oldenburg D W.Automatic phase correction of common－midpoint stacked data.Geophysics，1987，52（1）：51－59.

［6］ Wang Y H.A stable and efficient approach of inverseQfiltering.Geophysics，2002，67（2）：657－663.

［7］ Wang Y H.Inverse Q－filter for seismic resolution enhancement.Geophysics，2006，71（3）：V51－V60.

［8］ Yuan S Y，Wang S X.Influence of inaccurate wavelet phase estimation on seismic inversion.AppliedGeophysics，2011，8（1）：48－59.

［9］ 梁光河.地震子波提取方法研究.石油物探，1998，37（1）：31－39.

Liang G H.On the methods of seismic wavelet extraction.GeophysicalProspectingforPetroleum（in Chinese），1998，37（1）：3l－39.

［10］ 戴永壽，王俊嶺，王偉偉等.基于高階累積量ARMA模型線性非線性結(jié)合的地震子波提取方法研究.地球物理學(xué)報，2008，51（6）：1851－1859.

Dai Y S，Wang J L，Wang W W，et al.Seismic wavelet extraction via cumulant－based ARMA model approach with linear and nonlinear combination.ChineseJ.Geophys.（in Chinese），2008，51（6）：1851－1859.

［11］ Robinson E A.Predictive decomposition of time series with application to seismic exploration.Geophysics，1967，32（3）：418－484.

［12］ 吳鎮(zhèn)揚.數(shù)字信號處理.北京：高等教育出版社，2004.

Wu Z Y.Digital Signal Processing（in Chinese）.Beijing：Higher Education Press，2004.

［13］ Wiggins R A.Minimum entropy deconvolution.Geoexploration，1978，16（1－2）：21－35.

［14］ Sacchi M D.Re－weighting strategies in seismic deconvolution.GeophysicalJournalInternational，1997，129：651－656.

［15］ 章珂，李衍達(dá)，劉貴忠等.多分辨率地震信號反褶積.地球物理學(xué)報，1999，42（4）：529－535.

Zhang K，Li Y D，Liu G Z，et al.Multiresolution seismic signal deconvolution.ChineseJ.Geophys.（in Chinese），1999，42（4）：529－535.

［16］ Zhang R，Castagna J.Seismic sparse－layer reflectivity inversion using basis pursuit decomposition.Geophysics，2011，76（6）：R147－R158.

［17］ Yuan S Y，Wang S X.Noise attenuation without spatial assumptions about seismic coherent events.／／80th Ann.Internat.Mtg.，Soc.Explor.Geophys.，Expanded Abstracts，2010，3524－3528.

［18］ Wiggins R.Entropy guided deconvolution.Geophysics，1985，50（12）：2720－2726.

［19］ 張廣智，劉洪，印興耀.井旁道地震子波精細(xì)提取方法.石油地球物理勘探，2005，40（2）：158－162.

Zhang G Z，Liu H，Yin X Y.Method for fine picking up seismic wavelet at up－h(huán)ole trace.OilGeophysicalProspecting（in Chinese），2005，40（2）：158－162.