本發(fā)明屬于地球物理勘探領(lǐng)域，涉及一種地下介質(zhì)物性參數(shù)估計(jì)方法，特別涉及一種數(shù)據(jù)驅(qū)動(dòng)的地下介質(zhì)q場(chǎng)估計(jì)方法。

背景技術(shù)：

地震波在地下介質(zhì)中傳播時(shí)會(huì)隨著傳播時(shí)間衰減，主要表現(xiàn)為振幅衰減、相位畸變、主頻降低，且高頻部分比低頻部分衰減快，淺層比深層衰減快。引起衰減的因素包括幾何擴(kuò)散、散射、介質(zhì)的非完全彈性等。其中由介質(zhì)的非完全彈性引起的衰減常用品質(zhì)因子q來(lái)描述。如果知道地下介質(zhì)q值，可對(duì)地震記錄進(jìn)行有效的反q濾波補(bǔ)償，從而使地震記錄的高頻成分得到增強(qiáng)，地下淺、中、深層反射地震記錄的縱向分辨率都得到提高；另外，q值也與巖性、飽和度、孔隙率等有關(guān)，可用于儲(chǔ)層識(shí)別和烴類檢測(cè)。因此，估計(jì)地下介質(zhì)q值具有重要的實(shí)際意義。

迄今，國(guó)內(nèi)外學(xué)者提出了多種直接估計(jì)q值的方法。時(shí)間域方法有子波模擬法、解析信號(hào)法和脈沖上升時(shí)間法等，這些方法都需要振幅保真度高的數(shù)據(jù)，而實(shí)際記錄常受到幾何擴(kuò)散、散射等因素的影響，致使這些時(shí)間域方法精度降低。頻率域常用的方法有對(duì)數(shù)譜比值法(lsr,log-spectrumratio)、譜匹配法、中心頻率偏移法等，這些方法需要用時(shí)窗截取地震記錄，再對(duì)截出的地震記錄進(jìn)行分析。這類方法時(shí)窗的選擇對(duì)估計(jì)結(jié)果影響較大，如何最佳地選擇窗的類型及長(zhǎng)度是一個(gè)難題。利用脈沖瞬時(shí)屬性的變化來(lái)估計(jì)q值也是一條重要的途徑。mathneey和nowack提出了瞬時(shí)頻率匹配法，高靜懷和楊森林等人提出利用小波域地震子波包絡(luò)峰值處瞬時(shí)頻率估計(jì)q值的方法wepif(waveletenvelopepeakinstantaneousfrequency)。這類方法的優(yōu)點(diǎn)是時(shí)窗易于選取或者不需要加時(shí)窗，估計(jì)得到的q值分辨率較高，但由于只利用了脈沖包絡(luò)峰值處單個(gè)或相鄰幾個(gè)數(shù)據(jù)值，估計(jì)得到的q值受隨機(jī)噪聲的影響很大。tonn對(duì)近十種直接估計(jì)q值的方法進(jìn)行了比較，結(jié)果表明每種方法都有一定的適用條件。

由于垂直地震剖面(vsp)資料可對(duì)穿過(guò)介質(zhì)的波進(jìn)行直接測(cè)量，這為上述直接估計(jì)介質(zhì)q值的方法以及利用波形反演估計(jì)介質(zhì)q值及其它參數(shù)的方法(如：stewart提出的同時(shí)利用相鄰四道記錄的全波形在頻率域反演速度、衰減因子、上行波和下行波等的方法；amundsen和mittet只利用初至下行波和一次反射波在頻率域反演相速度和品質(zhì)因子q的方法)提供了途徑。vsp是在地震測(cè)井基礎(chǔ)上發(fā)展起來(lái)的，而測(cè)井?dāng)?shù)量有限，故由vsp資料估計(jì)得到的q值不能反映整個(gè)地下介質(zhì)的q場(chǎng)分布。為此，學(xué)者們也嘗試研究利用反射地震資料估計(jì)地下介質(zhì)q場(chǎng)的方法。目前常用的各種q值估計(jì)(或反演)方法在用于反射地震資料時(shí)往往會(huì)碰到許多問(wèn)題，如：復(fù)雜的地質(zhì)結(jié)構(gòu)使得反射波互相交疊，不能直接從地震記錄獲得時(shí)變子波的波形、振幅譜以及瞬時(shí)頻率等。另外，噪聲干擾影響瞬時(shí)屬性的提取，進(jìn)而影響q值估計(jì)的準(zhǔn)確性。也就是說(shuō)，復(fù)雜的地質(zhì)結(jié)構(gòu)和噪聲干擾常常使得估計(jì)出的地下介質(zhì)的q值不穩(wěn)定。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明的目的在于克服上述現(xiàn)有技術(shù)的缺點(diǎn)，提供一種數(shù)據(jù)驅(qū)動(dòng)的地下介質(zhì)q場(chǎng)估計(jì)方法，該方法為一種數(shù)據(jù)驅(qū)動(dòng)的基于自適應(yīng)分子分解利用非線性壓縮映射提取時(shí)變子波振幅譜并用質(zhì)心頻率偏移法估計(jì)地下介質(zhì)q場(chǎng)的方法，該方法采用地震包絡(luò)局部峰值作為地層結(jié)構(gòu)約束來(lái)構(gòu)造自適應(yīng)于地層結(jié)構(gòu)的分子時(shí)窗，然后采用非線性壓縮映射子波振幅譜估計(jì)方法從地震道分子分解時(shí)頻譜中估計(jì)正比于時(shí)變子波振幅譜的分量，接著計(jì)算質(zhì)心頻率，并對(duì)質(zhì)心頻率做篩選，減少子波干涉和噪聲的影響，最后用質(zhì)心頻率偏移法估計(jì)得到穩(wěn)健的地下介質(zhì)q場(chǎng)。

為實(shí)現(xiàn)上述目的，本發(fā)明的目的是通過(guò)以下技術(shù)方案來(lái)解決的：

一種數(shù)據(jù)驅(qū)動(dòng)的地下介質(zhì)q場(chǎng)估計(jì)方法，采用地震包絡(luò)局部峰值作為地層結(jié)構(gòu)約束來(lái)構(gòu)造自適應(yīng)于地層結(jié)構(gòu)的分子時(shí)窗，然后采用非線性壓縮映射子波振幅譜估計(jì)方法從地震道分子分解時(shí)頻譜中估計(jì)正比于時(shí)變子波振幅譜的分量，接著計(jì)算質(zhì)心頻率，并對(duì)質(zhì)心頻率做篩選，減少子波干涉和噪聲的影響，最后用質(zhì)心頻率偏移法估計(jì)得到穩(wěn)健的地下介質(zhì)q場(chǎng)。

本發(fā)明進(jìn)一步的改進(jìn)在于，具體包括以下步驟：

1)提取地震道包絡(luò)局部峰值：

2)生成滿足單位分割的原子窗集：

選擇基本原子窗函數(shù)g(t)，令gj(t)＝g(t-jδt)表示中心位于第j個(gè)采樣點(diǎn)的原子窗，對(duì)原子窗族{gj:1≤j≤n}按下式歸一化

得一組單位分割原子窗集{gj:1≤j≤n}，這里n為地震道采樣點(diǎn)的個(gè)數(shù)；

3)構(gòu)造自適應(yīng)分子時(shí)窗：

選擇每相鄰兩個(gè)包絡(luò)局部峰值點(diǎn)的中點(diǎn)作為各分子窗的邊界點(diǎn)，將邊界點(diǎn)間的滿足單位分割的原子窗疊加起來(lái)就形成了分子窗；設(shè)第k個(gè)分子窗對(duì)應(yīng)的包絡(luò)局部峰值點(diǎn)位于pk，前一個(gè)包絡(luò)局部峰值點(diǎn)位于pk-1(p0＝-p1)，后一個(gè)包絡(luò)局部峰值點(diǎn)位于pk+1，則此分子窗ψk(t)，對(duì)應(yīng)的第一個(gè)原子窗的中心位于mk-1+1＝(pk-1+pk)/2+1(m0＝0)，對(duì)應(yīng)的最后一個(gè)原子窗的中心位于mk＝(pk+pk+1)/2，分子窗由這中間的mk-mk-1個(gè)原子窗疊加得到，即分子窗ψk(t)由下式表示

令l為分子窗的個(gè)數(shù)，分子窗族{ψk(t):1≤k≤l}也構(gòu)成單位分割；

4)對(duì)步驟3)得到的自適應(yīng)分子時(shí)窗進(jìn)行能量歸一化：

令ek表示第k個(gè)分子窗的能量，即

能量歸一化以后的分子分析時(shí)窗為{ψk(t)/ek:1≤k≤l}；

對(duì)分子分析時(shí)窗進(jìn)行平移和調(diào)制后，得到一組分子標(biāo)架；此時(shí)，分子分解時(shí)頻變換定義為

其中f為頻率；

5)利用非線性壓縮映射提取時(shí)變子波振幅譜：

設(shè)x0∈(a,b)，取δ0＞0，使得對(duì)于任意給定的δ＞0，定義l²[a,b]空間中的子集

在[a,x0-δ0]上單增,且在[x0+δ0,b]上單減}

如果且q≤2，則有任取定義函數(shù)

此處0＜q≤1；顯然，0≤fq(u；x)≤1；設(shè)cq＞0，α,β＞1，對(duì)于定義上的算子p為

p是到的非線性壓縮映射算子；

對(duì)于第k個(gè)分子窗中的地震記錄片段，其振幅譜為設(shè)定迭代初值為：

其中0＜q≤1，

其中fm為振幅譜的截止頻率；建立迭代uk＝p(uk-1)，其中算子參數(shù)由最小二乘得到；由于算子是壓縮映射，得到不動(dòng)點(diǎn)u^*，則估計(jì)的子波振幅譜為記為lk(f)；

6)計(jì)算質(zhì)心頻率：

對(duì)于第k個(gè)分子窗中的地震記錄片段，質(zhì)心頻率fc,k為

式中fc為子波振幅譜的截止頻率；對(duì)計(jì)算得到的質(zhì)心頻率處理后得到最終的質(zhì)心頻率

7)估算地下介質(zhì)q場(chǎng)：

用質(zhì)心頻率偏移法估算q值，相應(yīng)的估算公式如下

式中：分別為t1時(shí)刻頻譜的質(zhì)心頻率和方差；是t2時(shí)刻頻譜的質(zhì)心頻率；時(shí)差δt＝t2-t1；為t1和t2之間地層介質(zhì)的品質(zhì)因子；進(jìn)一步推導(dǎo)，得到t時(shí)刻地下介質(zhì)的等效q值

式中fc,0和為初始時(shí)刻地震子波振幅譜的質(zhì)心頻率和方差，即

用步驟6)中計(jì)算得到的質(zhì)心頻率代替fc(t)，得到各個(gè)窗中心點(diǎn)處的q值

tk表示第k個(gè)分子窗的中心點(diǎn)所對(duì)應(yīng)的時(shí)間位置；對(duì)窗中心點(diǎn)處估得的q值插值后即得到穩(wěn)定的地下介質(zhì)q場(chǎng)。

本發(fā)明進(jìn)一步的改進(jìn)在于，步驟1)中提取地震道包絡(luò)局部峰值具體過(guò)程為：設(shè)s^*(t)為地震道s(t)的hilbert變換，則

a(t)＝[s(t)²+s^*(t)²]^1/2

為地震道s(t)的包絡(luò)；由上式計(jì)算地震道的包絡(luò)，并提取包絡(luò)局部峰值點(diǎn)。

本發(fā)明進(jìn)一步的改進(jìn)在于，步驟6)中對(duì)計(jì)算得到的質(zhì)心頻率做多次光滑濾波，并計(jì)算濾波前后的絕對(duì)誤差，去掉絕對(duì)誤差較大的一些質(zhì)心頻率后，再次擬合濾波，得到最終的質(zhì)心頻率

與現(xiàn)有技術(shù)相比，本發(fā)明具有以下有益效果：

本發(fā)明利用地震包絡(luò)峰值作為地層結(jié)構(gòu)約束來(lái)構(gòu)造自適應(yīng)于地層結(jié)構(gòu)的分子時(shí)窗，由于地震數(shù)據(jù)的包絡(luò)峰值可以大致反映地層的層序結(jié)構(gòu)，因此構(gòu)造的時(shí)窗在橫向上與地層結(jié)構(gòu)有關(guān)，有利于保持估計(jì)的q場(chǎng)的橫向連續(xù)性，并且這種自適應(yīng)分子時(shí)窗構(gòu)造方法可以保證每個(gè)分子窗內(nèi)至少有一個(gè)反射子波，在減少時(shí)窗數(shù)量、提高計(jì)算效率的同時(shí)，可以有效減少窗端點(diǎn)對(duì)子波的截?cái)?；采用非線性壓縮映射子波振幅譜估計(jì)方法可以從較短數(shù)據(jù)中估計(jì)子波振幅譜，這為本發(fā)明從地震道分子分解時(shí)頻譜中估計(jì)正比于時(shí)變子波振幅譜的分量，從而得到穩(wěn)定的質(zhì)心頻率提供了保證；在計(jì)算q值之前對(duì)質(zhì)心頻率做了濾波篩選，可以有效減少子波干涉和隨機(jī)噪聲的影響，得到穩(wěn)定的地下介質(zhì)q場(chǎng)。

附圖說(shuō)明

圖1是合成地震道時(shí)反射系數(shù)示意圖；

圖2是合成地震道的包絡(luò)局部峰值示意圖；其中，實(shí)線表示非平穩(wěn)合成地震道，虛線表示地震道包絡(luò)，星號(hào)表示包絡(luò)局部峰值；

圖3是構(gòu)成單位分割的原子gabor時(shí)窗示例圖；

圖4是地震道包絡(luò)局部峰值約束下尋找分子gabor時(shí)窗邊界點(diǎn)示例圖；

圖5是將邊界點(diǎn)間的原子gabor時(shí)窗疊加形成分子gabor時(shí)窗示例圖；

圖6是能量歸一化后得到的分子gabor時(shí)窗示例圖；

圖7是一段實(shí)際地震資料剖面圖；

圖8是從地震資料剖面中提取出的質(zhì)心頻率圖；

圖9是對(duì)質(zhì)心頻率做濾波篩選后得到的質(zhì)心頻率圖；

圖10是基于濾波篩選后的質(zhì)心頻率計(jì)算得到的地下介質(zhì)q場(chǎng)圖。

具體實(shí)施方式

下面結(jié)合附圖對(duì)本發(fā)明做進(jìn)一步詳細(xì)描述：

本發(fā)明針對(duì)復(fù)雜的地質(zhì)結(jié)構(gòu)和噪聲干擾常常使得估計(jì)出的地下介質(zhì)的q值不穩(wěn)定的問(wèn)題，以地震道包絡(luò)局部峰值作為地層結(jié)構(gòu)，提出基于自適應(yīng)分子分解利用非線性壓縮映射提取時(shí)變子波振幅譜并用質(zhì)心頻率偏移法估計(jì)得到穩(wěn)健的地下介質(zhì)q場(chǎng)的方法。

本發(fā)明的物質(zhì)基礎(chǔ)是疊后反射地震數(shù)據(jù)，即在地表設(shè)置震源，在地表安置檢波器接收地震波，經(jīng)過(guò)一系列處理，最終疊加或偏移成像得到的地震數(shù)據(jù)。本發(fā)明的具體步驟為：

1)提取地震道包絡(luò)局部峰值；

根據(jù)復(fù)道分析原理，設(shè)s^*(t)為地震道s(t)的hilbert變換，則

a(t)＝[s(t)²+s^*(t)²]^1/2

為地震道s(t)的包絡(luò)。對(duì)比圖1中合成地震道所用的反射系數(shù)和圖2中合成地震道包絡(luò)局部峰值點(diǎn)可見，地震道的包絡(luò)局部峰值與反射界面有一定的對(duì)應(yīng)關(guān)系。

2)生成滿足單位分割的原子窗集；

恰當(dāng)?shù)剡x擇基本原子窗函數(shù)g(t)，令gj(t)＝g(t-jδt)表示中心位于第j個(gè)采樣點(diǎn)的原子窗，對(duì)原子窗族{gj:1≤j≤n}按下式歸一化

可得一組單位分割原子窗集{gj:1≤j≤n}，這里n為地震道采樣點(diǎn)的個(gè)數(shù)。如圖3中的一組窗就是滿足單位分割的原子窗集。

3)構(gòu)造自適應(yīng)分子時(shí)窗；

選擇每相鄰兩個(gè)包絡(luò)局部峰值點(diǎn)的中點(diǎn)作為各分子窗的邊界點(diǎn)，將邊界點(diǎn)間的滿足單位分割的原子窗疊加起來(lái)就形成了分子窗。如圖4所示，設(shè)第k個(gè)分子窗對(duì)應(yīng)的包絡(luò)局部峰值點(diǎn)位于pk，前一個(gè)包絡(luò)局部峰值點(diǎn)位于pk-1(p0＝-p1)，后一個(gè)包絡(luò)局部峰值點(diǎn)位于pk+1，則此分子窗ψk(t)，對(duì)應(yīng)的第一個(gè)原子窗的中心位于mk-1+1＝(pk-1+pk)/2+1(m0＝0)，對(duì)應(yīng)的最后一個(gè)原子窗的中心位于mk＝(pk+pk+1)/2，分子窗由這中間的mk-mk-1個(gè)原子窗疊加得到，即分子窗ψk(t)可由下式表示

令l為分子窗的個(gè)數(shù)，易知分子窗族{ψk(t):1≤k≤l}也構(gòu)成單位分割，如圖5所示。這樣，可以保證每個(gè)分子窗內(nèi)至少有一個(gè)反射子波，在減少時(shí)窗數(shù)量的同時(shí)，可以有效減少窗端點(diǎn)對(duì)子波的截?cái)嘈?yīng)。并且由于包絡(luò)峰值可以大致反映地層的層序結(jié)構(gòu)，構(gòu)成的分子窗在橫向上與地層結(jié)構(gòu)有關(guān)，有利于保持處理后資料的橫向連續(xù)性。

4)對(duì)步驟3)得到的自適應(yīng)分子時(shí)窗進(jìn)行能量歸一化；

分子分解時(shí)頻變換是保能量的變換，由圖5可以看出，在步驟3)得到的分子窗族中，不同的窗能量不同。如果用這些窗直接構(gòu)造分子gabor標(biāo)架，則一個(gè)地震道的時(shí)頻能量表示不僅與該地震道有關(guān)，還與分子窗有關(guān)。希望這一時(shí)頻能量表示只與地震道有關(guān)，為此，需要對(duì)每個(gè)分子窗進(jìn)行能量歸一化。令ek表示第k個(gè)分子窗的能量，即

能量歸一化以后的分子分析時(shí)窗為{ψk(t)/ek:1≤k≤l}，如圖6所示。

對(duì)分子分析時(shí)窗進(jìn)行平移和調(diào)制后，可得到一組分子標(biāo)架。此時(shí)，分子分解時(shí)頻變換可定義為

其中f為頻率。

5)利用非線性壓縮映射提取時(shí)變子波振幅譜；

設(shè)x0∈(a,b)，取δ0＞0足夠小，使得對(duì)于任意給定的δ＞0，定義l²[a,b]空間中的子集

在[a,x0-δ0]上單增,且在[x0+δ0,b]上單減}

如果且q≤2，則有任取定義函數(shù)

此處0＜q≤1。顯然，0≤fq(u；x)≤1。設(shè)cq＞0，α,β＞1，對(duì)于定義上的算子p為

顯然，p[u；x]在區(qū)間[a,b]上是連續(xù)的，并且δ≤p[u；x]≤cq+δ，由此可得p[u；x]∈l²[a,b]。因此，p是到l²[a,b]的非線性算子。根據(jù)banach空間算子不動(dòng)點(diǎn)定理，可以證明：p是到的壓縮算子，并且在中存在唯一不動(dòng)點(diǎn)。

對(duì)于第k個(gè)分子窗中的地震記錄片段，其振幅譜為設(shè)定迭代初值為：

其中0＜q≤1，

其中fm為振幅譜的截止頻率。建立迭代uk＝p(uk-1)，其中算子參數(shù)由最小二乘得到。由于算子是壓縮映射，可以得到不動(dòng)點(diǎn)u^*，則估計(jì)的子波振幅譜為記為lk(f)。

6)計(jì)算質(zhì)心頻率；

對(duì)于第k個(gè)分子窗中的地震記錄片段，質(zhì)心頻率為

式中fc為子波振幅譜的截止頻率。理論上，地震子波在傳播的過(guò)程中質(zhì)心頻率是隨著傳播時(shí)間衰減的，然而實(shí)際中(如圖7所示)，當(dāng)?shù)貙虞^薄時(shí)，地震波之間的互相干涉會(huì)使得質(zhì)心頻率突然增大或者突然減小(如圖8所示)。為此，本發(fā)明在計(jì)算q值之前對(duì)計(jì)算得到的質(zhì)心頻率做多次光滑濾波，并計(jì)算濾波前后的絕對(duì)誤差，去掉絕對(duì)誤差較大的一些質(zhì)心頻率后，再次擬合濾波，得到最終的質(zhì)心頻率以減少地震波干涉的影響(如圖9所示)。

7)估算地下介質(zhì)q場(chǎng)；

用質(zhì)心頻率偏移法估算q值，相應(yīng)的估算公式如下

式中：分別為t1時(shí)刻頻譜的質(zhì)心頻率和方差；是t2時(shí)刻頻譜的質(zhì)心頻率；時(shí)差δt＝t2-t1；為t1和t2之間地層介質(zhì)的品質(zhì)因子。進(jìn)一步推導(dǎo)，可得t時(shí)刻地下介質(zhì)的等效q值

式中fc,0和為初始時(shí)刻地震子波振幅譜的質(zhì)心頻率和方差，即

用步驟6)中計(jì)算得到的質(zhì)心頻率代替fc(t)，可得到各個(gè)窗中心點(diǎn)處的q值

tk表示第k個(gè)分子窗的中心點(diǎn)所在的時(shí)間位置。對(duì)窗中心點(diǎn)處估得的q值插值后，就可以得到穩(wěn)定的地下介質(zhì)q場(chǎng)(如圖10所示)。

該方法用地震道包絡(luò)局部峰值作為地層結(jié)構(gòu)約束來(lái)構(gòu)造自適應(yīng)分子時(shí)窗，并用這些時(shí)窗生成分子標(biāo)架對(duì)地震記錄做分子分解；接著在分子分解時(shí)頻域利用非線性壓縮映射提取時(shí)變子波振幅譜，并計(jì)算質(zhì)心頻率；然后對(duì)質(zhì)心頻率做濾波篩選，以減少子波干涉和噪聲影響；最后用質(zhì)心頻率偏移法計(jì)算每個(gè)時(shí)窗對(duì)應(yīng)的等效衰減因子q，并插值得到穩(wěn)定的地下介質(zhì)q場(chǎng)。所述分子時(shí)窗的構(gòu)造方法，可以保證每個(gè)時(shí)窗內(nèi)至少有一個(gè)反射子波，在減少時(shí)窗數(shù)量、提高計(jì)算效率的同時(shí)，使得這些時(shí)窗在橫向上與地層結(jié)構(gòu)有關(guān)，有利于保持所估計(jì)q場(chǎng)的橫向連續(xù)性。本發(fā)明能夠有效克服復(fù)雜地質(zhì)結(jié)構(gòu)變化和噪聲干擾的影響，得到穩(wěn)定的地下介質(zhì)q場(chǎng)。