地震子波提取方法評(píng)述

1地震子波提取方法地震波提取是地震勘探中的一個(gè)非常重要的問(wèn)題。在正演問(wèn)題中,需要通過(guò)波動(dòng)方程或褶積模型結(jié)合地震子波來(lái)形成正演模擬地震數(shù)據(jù)。在反演和反褶積問(wèn)題中,也需要通過(guò)地震道提取一個(gè)子波,不同的子波往往對(duì)反演結(jié)果會(huì)有不同的影響。地震子波為具有確定的起始時(shí)間、能量有限且具有一定延續(xù)長(zhǎng)度的信號(hào),它是地震記錄中的基本單元。一個(gè)地震子波一般具有2～3個(gè)相位的延續(xù)長(zhǎng)度(約為90ms)。一般認(rèn)為,地震炸藥震源激發(fā)時(shí)所產(chǎn)生的地震波僅是一個(gè)延續(xù)時(shí)間極短的尖脈沖,隨著尖脈沖在黏彈性介質(zhì)中傳播,尖脈沖的高頻成分很快衰減,波形隨之加寬,便形成了地震子波。地震子波提取的基本框架是褶積模型,它包括子波、反射系數(shù)序列、含有噪聲的地震道。地震子波提取方法主要分為確定性子波提取方法、統(tǒng)計(jì)性子波提取方法兩大類。確定性子波提取方法首先利用測(cè)井資料計(jì)算出反射系數(shù)序列,然后結(jié)合井旁地震道由褶積模型求出地震子波。統(tǒng)計(jì)性子波提取方法通過(guò)地震道自身來(lái)估計(jì)子波,又可分為基于二階統(tǒng)計(jì)量和基于高階統(tǒng)計(jì)量?jī)煞N方法。基于二階統(tǒng)計(jì)量方法首先由Robinson(1975)提出,該法基于這樣的假設(shè):地震子波是時(shí)不變的且是單一相位的(零相位、最小相位或最大相位),地下反射系數(shù)是具有白噪譜的隨機(jī)序列。因此在上述假設(shè)下觀測(cè)到的地震道的自相關(guān)就給定了地震子波自相關(guān)的一個(gè)估計(jì)(即已知地震子波的振幅譜)。實(shí)際上地震子波是一種混合相位的,因此基于二階統(tǒng)計(jì)量的自相關(guān)統(tǒng)計(jì)方法提取的地震子波也是不準(zhǔn)確的。自20世紀(jì)80年代后期,許多學(xué)者開(kāi)始使用高階統(tǒng)計(jì)方法解決地震子波估計(jì)問(wèn)題,這些方法大多源于20世紀(jì)60年代發(fā)展起來(lái)的累積量和多譜理論,Matsuoka等人(1984)最早將它們用于混合相位地震子波的估計(jì)。高階統(tǒng)計(jì)量地震子波提取的新思想是由Lazear(1993)和Velis等人(1996)提出的,他們將非高斯信號(hào)處理中的四階累積量用于子波的估計(jì),為解決混合相位子波估計(jì)問(wèn)題提供了一條全新的思路。2密度測(cè)井資料確定性子波提取方法是利用聲波測(cè)井和密度測(cè)井資料,首先計(jì)算出反射系數(shù)序列,然后結(jié)合井旁地震道由褶積模型求出地震子波。其代表算法有以下幾種。2.1維納濾波法它是在時(shí)間域內(nèi)通過(guò)解線性方程組來(lái)估計(jì)子波,實(shí)際計(jì)算過(guò)程中只能截取一段地震記錄來(lái)估計(jì)子波,這就會(huì)產(chǎn)生截?cái)嗾`差。2.2除法律糾紛這種方法遇到的一個(gè)問(wèn)題是可能被零除,因?yàn)榉瓷湎禂?shù)譜中有零值或接近于零的值。習(xí)慣上是將反射系數(shù)序列譜白噪化,以避免這種現(xiàn)象的發(fā)生。2.3模式1:w選取目標(biāo)函數(shù)通過(guò)正演模型計(jì)算出合成地震記錄,然后與實(shí)際地震記錄比較,逐步修正子波,直到合成記錄與實(shí)際觀測(cè)記錄完全匹配為止。該法使用如下最小化目標(biāo)函數(shù)minf(w)=(r*w-d)(r*w-d)+Q(w)(1)其中:w為待求的地震子波;r為地層反射系數(shù)序列;d為地震數(shù)據(jù);Q(w)為對(duì)子波的先驗(yàn)約束。2.4的先驗(yàn)分布p該方法應(yīng)用貝葉斯公式將地震子波、反射系數(shù)、地震道噪聲的先驗(yàn)信息與似然函數(shù)相結(jié)合,從而得到各個(gè)參數(shù)的后驗(yàn)概率密度。假定要估計(jì)的地震子波為θ,為了進(jìn)行貝葉斯推斷,首先要來(lái)設(shè)定θ的先驗(yàn)分布p(θ),然后將得到的觀測(cè)信息融入先驗(yàn)信息來(lái)改進(jìn)先驗(yàn)分布。假定地震數(shù)據(jù)為x(與θ有關(guān)),用p(x|θ)表示其條件密度函數(shù)。有了先驗(yàn)分布和條件分布,就可以通過(guò)貝葉斯公式得到參數(shù)θ的后驗(yàn)分布,即p(θ|x)=h(x,θ)p(x)∝p(θ)p(x|θ)(2)通過(guò)最大后驗(yàn)估計(jì)或馬爾可夫鏈—蒙特卡羅模擬方法得到地震子波。2.5井旁地震子波的計(jì)算該方法考慮到地震子波在空間變化的特點(diǎn),首先由多道相關(guān)法提取初始地震子波的振幅譜,然后結(jié)合測(cè)井資料確定初始地震子波的相位譜,最后根據(jù)離散反演理論迭代求取精細(xì)的井旁地震子波。除上述五種確定性子波提取方法之外,還有很多其他的確定性地震子波提取方法,這里就不再一一贅述。3系統(tǒng)函數(shù)w都未知從盲信號(hào)處理的角度出發(fā),統(tǒng)計(jì)性子波提取應(yīng)該屬于盲系統(tǒng)辨識(shí)的范疇(圖1)。盲系統(tǒng)辨識(shí)是指在輸入信號(hào)R和系統(tǒng)函數(shù)W都未知的情況下,僅僅根據(jù)系統(tǒng)的輸出D來(lái)辨識(shí)系統(tǒng)函數(shù)W。統(tǒng)計(jì)性地震子波提取與盲系統(tǒng)辨識(shí)問(wèn)題類似,如果將地層反射系數(shù)看作輸入R,將系統(tǒng)地震子波看作是系統(tǒng)函數(shù)W,統(tǒng)計(jì)性地震子波提取就是在地層反射系數(shù)和地震子波都未知的情況下,僅僅根據(jù)觀測(cè)到的地震記錄D來(lái)估計(jì)地震子波。有關(guān)統(tǒng)計(jì)性地震子波主要提取方法分述如下。3.1相位最小波提取3.1.1子波時(shí)間域的反變?cè)诜瓷湎禂?shù)為白噪聲的假設(shè)條件下,地震記錄的自相關(guān)等于子波自相關(guān),通過(guò)地震記錄的自相關(guān)得到子波的振幅譜。當(dāng)子波是零相位時(shí),可以直接由地震記錄的振幅譜得到子波的頻譜,再由傅里葉反變換得到子波的時(shí)間域形式。當(dāng)子波為最小相位時(shí),可以由希爾伯特變換法得到相位譜,記子波的振幅譜為A(ω),待求的相位譜為φ(ω),則φ(ω)=-ΗΤ[A(ω)]=-lnA(ω)*1πω(3)最小相位譜也可以用Wold-Kolmogorov公式得到,然后與子波振幅譜組合再通過(guò)傅里葉反變換得到最小相位子波。3.1.2地震正演模型計(jì)算當(dāng)z為相位的子波在反射系數(shù)為白噪聲的假設(shè)條件下,有rbb(τ)=rxx(τ)(4)其中:rbb(τ)為子波的自相關(guān);rxx(τ)為地震道的自相關(guān)。對(duì)式(4)兩邊同時(shí)做Z變換Rbb(z)=B(z)B(z-1)=Rxx(z)=Μ∑τ=-Μrxx(τ)z-τ(5)對(duì)上式進(jìn)行因式分解,一般可以得到2M個(gè)互為倒數(shù)的零點(diǎn),設(shè)單位圓內(nèi)的零點(diǎn)為zm(m=1,?,k)?1zm(m=k+1,?,Μ),則B(Ζ)=Μk∏m=1(z-zm)Μ∏m=k+1(z-1zm)zΜ=b0+b1z-1+b2z-2+?+bΜz-Μ(6)由此可以得到最小相位地震子波b(t)=(b0,b1,?,bΜ)3.1.3相方位波的計(jì)算通過(guò)地震道的四階矩,應(yīng)用脈沖反褶積方法求得一個(gè)反子波,然后由反子波得到最小相位子波的振幅譜和相位譜,最終得到最小相位子波。該方法可以有效抑制地震道中的高斯白噪聲和高斯有色噪聲。3.1.4確定相位譜并求取地震子波在假設(shè)地震記錄所反映的波阻抗函數(shù)具有分形特性的條件下,反射系數(shù)的功率譜與頻率的關(guān)系為Ρr(f)=cf-β經(jīng)過(guò)整理,可以得到Ρw(f)=bΡx(f)fβ(7)對(duì)上式兩邊開(kāi)方并取絕對(duì)值得|W(f)|=|√bΡx(f)fβ|然后應(yīng)用同態(tài)濾波的方法結(jié)合希爾伯特變換得到相位譜,進(jìn)而求取地震子波。3.1.5arma給力量計(jì)算如果將地震數(shù)據(jù)表示為ARMA過(guò)程,即x(n)+p∑i=1aix(n-i)=e(n)+q∑j=1bje(n-j)(8)其中:x(n)為地震數(shù)據(jù);e(n)為地層反射系數(shù);ai、bj為系統(tǒng)參數(shù)?？梢缘玫絉x(τ)=σ2∞∑i=0h(i)h(i+τ)(9)線性系統(tǒng)的沖擊響應(yīng)就是沖擊信號(hào)激勵(lì)該系統(tǒng)時(shí)的輸出響應(yīng),故ARMA過(guò)程的定義式為p∑i=0aih(n-i)=q∑k=0bkδ(n-k)=bn(10)綜合式(9)、式(10)得到下面的方程Rx(l)+p∑i=1aiRx(l-i)=0?l＞q(11)式(11)就是著名的修正Yule-Walker方程,特別地,對(duì)于嚴(yán)格AR(p)過(guò)程,上式簡(jiǎn)化為Rx(l)+p∑i=1aiRx(l-i)=0?l＞0(12)通過(guò)求解以上方程組,就可以確定AR參數(shù)a1,…,ap,即{σ2(b20+b21+?+b2q)=c0σ2(b0b1+?+bq-1bq)=c0?σ2(b0bq)=c0(13)式(13)中共有q+2個(gè)未知參數(shù),為了保證解的惟一性,通常在q+1個(gè)方程組中約定b0=1,再通過(guò)Newton-Raphson算法求解上述方程組,就可以確定MA參數(shù)b1,…,bq,與AR參數(shù)相結(jié)合,從而最終確定ARMA最小相位系統(tǒng)。3.2混合相位波的提取3.2.1任意高斯統(tǒng)計(jì)的特性高階統(tǒng)計(jì)量通常包括高階累積量、高階矩以及高階累積量譜和高階矩譜,此外還包括倒高階累積量譜(即倒多譜)。相對(duì)于二階統(tǒng)計(jì)量及其功率譜而言,高階統(tǒng)計(jì)量除了包含完整的相位信息外,還具有任意高斯統(tǒng)計(jì)的優(yōu)良特性,給我們提供了十分豐富的信息,可以辨識(shí)非因果、非最小相位、非線性系統(tǒng),可以抑制高斯或非高斯的有色噪聲,可以抽取不同于高斯信號(hào)的多種信號(hào)特征等,從而給地震信號(hào)處理領(lǐng)域開(kāi)辟了一個(gè)嶄新而廣闊的前景。3.2.1.高階累積量的算法(1)GM算法三階累積量FIR系統(tǒng)的GM方程為q∑i=0b2(i)Rx(m-i)=∈3q∑i=0b(i)c3x(m-i,m-i)-q≤m≤2q(14)其中∈3=σ2e/γ3e。四階累積量FIR系統(tǒng)的GM方程表示為q∑i=0b3(i)Rx(m-i)=∈4q∑i=0b(i)c4x(m-i,m-i)-q≤m≤2q(15)其中∈4=σ2e/γ4e。GM算法是一種同時(shí)使用相關(guān)函數(shù)和高階累積量的方法,因此該方法只適用于無(wú)加性噪聲的特殊情況。(2)Tugnait算法為了使GM算法適用于加性白噪聲的情況,需要增加新的方程{q∑i=1b2(i)c3y(m-i,m-i)--q∑i=0[∈b2(i)]Ry(m,i)=-c3y(m,m)-q≤m≤1q+1≤m≤2qq∑i=1c3y(i-m,q)-[∈b(q)]Ry(m)=-c3y(-m,q)q+1≤m≤2q(16)求解上述方程組,可以得到問(wèn)題的解。(3)累積量矩陣方程法上述算法的主要缺點(diǎn)是只能適用于加性白噪聲(高斯或非高斯),為了在理論上完全抑制高斯有色噪聲,就必需避免自相關(guān)函數(shù)的使用,而只使用高階累計(jì)量?？紤]下面的方程q∑i=0w(i)ck-1kd(q,i-τ)=ckd(τ,0)ck-3kd(q,0)ckd(q,q)τ=-q,?,0,?,q(17)其中w(i)為待估計(jì)的子波。上式可以寫(xiě)成矩陣的形式[C1C2]w=[c1c2](18)通過(guò)最小二乘法可以求得地震子波w(t)。(4)MA階數(shù)的確定考慮擴(kuò)展階的累積量矩陣CΜA,e=[ckd(0,0)ckd(1,0)ckd(q,0)ckd(1,0)???ckd(qe,0)?0](19)式中qe>q,ckd是地震數(shù)據(jù)的高階量,使用奇異值分解確定CMA,e的有效秩(等于q+1),即可獲得MA階數(shù)q的估計(jì)。3.2.1.子波環(huán)保譜圖同態(tài)法子波估計(jì)時(shí)從復(fù)賽譜域分離子波,假設(shè)前提是子波比較光滑,而反射系數(shù)卻很不光滑,從而在復(fù)賽譜上可以將兩者分開(kāi)。高階復(fù)賽譜法非最小相位子波估計(jì)主要采用倒雙譜和倒三譜?？紤]下面的公式Bkd(z1,z2)=γkrH(z1)H(z2)H(z-11z-12)(20)其中:Bkd(z1,z2)為地震數(shù)據(jù)的雙譜;H(z)為待求子波的Z變換。由于子波是非最小相位,傳遞函數(shù)可以寫(xiě)為Η(z)=AzlΗmin(z)Ηmax(z-1)(21)其中:A為常數(shù);l表示線性相移;Hmin(z)、Hmax(z-1)分別為子波的最小相位部分和最大相位部分。將式(21)代入式(20),則雙譜變?yōu)锽kd(z1,z2)=βHmin(z1)Hmin(z2)··Hmax(z1-1)Hmax(z2-1)·(22)·Hmin(z1-1z2-1)Hmax(z1z2)其復(fù)賽譜為b^(m,n)=Ζ-1[lnBkd(z1,z2)](23)最終可以得到下式{b^(m,0)=ω^min(m)m＞0b^(m,0)=ω^max(m)m＜0(24)式中ω^(n)是子波的復(fù)賽譜。由式(24)可以看出,通過(guò)兩個(gè)半軸就可以定義子波的復(fù)賽譜。3.2.1.2223333343r.34444444444444444444444444444444423444444444444444444444343434343423434343423232323232222344444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444相位估計(jì)法是通過(guò)雙譜或三譜估計(jì)子波的相位譜,然后與子波的振幅譜相結(jié)合,通過(guò)傅里葉反變換得到地震子波。下面考慮雙譜的情況。設(shè)地震記錄和反射系數(shù)的雙譜分別為Bd(ω1,ω2)和Br(ω1,ω2),系統(tǒng)的頻率響應(yīng)記為H(ω)=|H(ω)|e-jφ(ω),則可以得到下面的方程{|Bd(ω1,ω2)|=γ3r|Η(ω1)|??|Η(ω2)|?|Η*(ω1+ω2)|(25)ψd(ω1,ω2)=φ(ω1)+φ(ω2)-φ(ω1+ω2)式中ψd(ω1,ω2)是地震記錄的雙譜相位。根據(jù)觀測(cè)數(shù)據(jù)可以估計(jì)Bd(ω1,ω2),即{|Bd(ω1,0)|=γ3r|Η(ω1)|?|Η(0)|?|Η*(ω1)|=γ3r?|Η(0)|?|Η(ω1)|2(26)|Η(ω1)|=c|Bd(ω1,0)|考慮相位公式的離散形式ψd(i,j)=φ(i)+φ(j)-φ(i+j)(27)式(27)可以寫(xiě)成矩陣的形式A?=ψ(28)最終相位譜的解可表示為?=(AΤA)-1AΤψ(29)實(shí)際計(jì)算表明,A的條件數(shù)比較低,因?yàn)檫@個(gè)方法在加性噪聲的干擾下相當(dāng)穩(wěn)健,該方法可以推廣到三譜,只是涉及的矩陣運(yùn)算量要比雙譜大得多。3.2.1.地震記錄vmax混合相位子波提取基于地震道高階累積量與子波高階矩匹配的思想,需建立如下的目標(biāo)函數(shù)min?=∑τ1,τ2,τ3[α(τ1,τ2,τ3)c4d(τ1,τ2,τ3)--γ4rm4w(τ1,τ2,τ3)]2(30)使得vdown≤w≤vup∑w=0其中:c4d(τ1,τ2,τ3)為地震記錄的四階累計(jì)量;m4w(τ1,τ2,τ3)為待求子波的四階矩;γ4r為反射系數(shù)的峰態(tài);w為子波向量;vdown和vup為子波的取值范圍,并滿足子波的和為零;α(τ1,τ2,τ3)為三維窗函數(shù)。由式(30)可以看出,目標(biāo)函數(shù)?為地震子波的函數(shù),該目標(biāo)函數(shù)具有很強(qiáng)的非線性,對(duì)其求解是一個(gè)非線性多參數(shù)多極值的優(yōu)化問(wèn)題。Lazear曾經(jīng)使用了基于梯度下降的線性優(yōu)化技術(shù)求解目標(biāo)函數(shù);Velis對(duì)地震記錄的四階累計(jì)量加了一個(gè)窗函數(shù),并采用模擬退火算法優(yōu)化該目標(biāo)函數(shù),取得了一定的效果。3.2.1.混合相位子波的制備該方法與累積量匹配法類似,將待求的混合相位子波分解為最小相位與一個(gè)全通子波的褶積,首先求出最小相位子波,然后用累積量匹配法求出全通子波,最后將最小相位子波與全通子波褶積就得到了混合相位的子波。3.2.2基于二次統(tǒng)計(jì)的方法(1)擴(kuò)展的yle-現(xiàn)實(shí)Porsani等人(1998)提出了一種混合子波估計(jì)的方法,利用離散時(shí)間信號(hào)的Z變換,并假設(shè)混合相位子波有α個(gè)零點(diǎn)在單位圓內(nèi),可以通過(guò)求解擴(kuò)展的Yule-Walker(EYW)方程來(lái)求其反子波,EYW方程系數(shù)矩陣的對(duì)角線上是自相關(guān)函數(shù)的α延遲。(2)合相原位波