截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演方法
【專利摘要】本發(fā)明涉及一種截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演方法。該方法首先在低頻段內(nèi)，用截?cái)鄷r(shí)窗反演策略進(jìn)行全波形反演，得到的反演結(jié)果作為初始模型用于中頻段反演。然后將中頻段的反演結(jié)果作為初始模型，用于常規(guī)時(shí)間域全波形反演并得到最終反演結(jié)果。截?cái)鄷r(shí)窗反演策略節(jié)約了大量正演時(shí)間。低通濾波多尺度反演策略，可將地震數(shù)據(jù)分為不同頻段進(jìn)行反演。將截?cái)鄷r(shí)窗與低通濾波組合使用，可在缺失低頻和強(qiáng)高斯噪聲干擾情況下，得到較好的反演結(jié)果。模型測(cè)試結(jié)果表明，該組合策略在提高計(jì)算效率的同時(shí)緩解了周跳現(xiàn)象，且截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演方法在反演精度、計(jì)算效率、低頻依賴、抗噪能力等方面優(yōu)于常規(guī)時(shí)間域全波形反演。
【專利說明】
截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演方法
技術(shù)領(lǐng)域
[0001] 本發(fā)明涉及一種地震勘探的數(shù)據(jù)處理方法，尤其是截?cái)鄷r(shí)窗和低通濾波多尺度反 演策略。整個(gè)計(jì)算過程利用超隨機(jī)震源編碼策略進(jìn)行全波形反演，目的是加快整個(gè)全波形 反演的進(jìn)程。
【背景技術(shù)】：
[0002] 隨著石油工業(yè)的發(fā)展和勘探開發(fā)的不斷深入，從構(gòu)造勘探階段逐漸走向巖性勘探 階段，其勘探難度逐漸加大。為了順應(yīng)這種要求，地震數(shù)據(jù)處理技術(shù)不斷完善，在這期間全 波形反演迅速發(fā)展起來，并成為當(dāng)今地球物理界的研究熱點(diǎn)。20世紀(jì)80年代，Tarantola提 出了時(shí)間域全波形反演理方法，并將目標(biāo)函數(shù)對(duì)模型參數(shù)求導(dǎo)通過殘差反傳波場(chǎng)和正傳波 場(chǎng)互相關(guān)運(yùn)算得到，從而避免了求取Jacobi矩陣。20世紀(jì)90年代，Pratt將全波形反演推廣 到了頻率域，提出只需要幾個(gè)離散的頻率就可以得到高精度的反演結(jié)果，而且低頻到高頻 的反演策略可以解決陷入局部極小值的問題。在頻率域全波形反演的基礎(chǔ)上Shin等人發(fā)展 了 Laplace-Fourier域全波形反演，該方法能夠在缺失低頻的情況下可以獲得一個(gè)高精度 初始模型。Wang等人(2011)提出在GPU上基于CUDA編程的加速全波形反演，該方法有效地加 快了全波形反演的計(jì)算效率。
[0003] Bunks(1995)提出基于低通濾波的多尺度全波形反演，該方法指出先對(duì)觀測(cè)數(shù)據(jù) 進(jìn)行低通濾波，再對(duì)估計(jì)的震源子波用以相同濾波器處理，可在一定程度上解決了全波形 反演的周跳問題。但是利用該多尺度方法進(jìn)行全波形反演依然存在著一些問題，如:地層是 一個(gè)粘彈性介質(zhì)，而且波動(dòng)方程正演是一個(gè)復(fù)雜的過程，當(dāng)?shù)卣鸩ㄔ诘叵聜鞑サ倪^程中會(huì) 出現(xiàn)頻移效應(yīng)或者其他復(fù)雜變化。利用低通濾波之后的震源子波進(jìn)行正演，會(huì)出現(xiàn)正演記 錄和觀測(cè)記錄不匹配的問題，關(guān)于這個(gè)問題Bunks等人并沒有給出很好的解決方案。Chen和 Wu等人利用時(shí)間衰減濾波器同時(shí)處理觀測(cè)數(shù)據(jù)和正演數(shù)據(jù)，即在地震數(shù)據(jù)上乘上一個(gè)隨著 時(shí)間增加而衰減的指數(shù)函數(shù)，這樣先演淺層信息，忽略深層構(gòu)造。隨著反演進(jìn)行逐漸減小衰 減函數(shù)的影響，使得反演由淺層到深層逐層反演。該方法極大的減小了每一步的反演參數(shù)， 一定程度上解決了全波形反演的的周跳問題。但是該方法在正演過程中浪費(fèi)了大量的時(shí) 間，因?yàn)檎莸玫降恼萦涗浽俪艘砸粋€(gè)衰減函數(shù)，深層信息基本衰減至零。該方法沒有利 用到深層信息，卻浪費(fèi)了大量時(shí)間來正演深部波場(chǎng)，很大程度上降低全波形反演的計(jì)算效 率。

【發(fā)明內(nèi)容】
：
[0004] 本發(fā)明的目的是針對(duì)上述現(xiàn)有低通濾波之后數(shù)據(jù)不匹配問題，提供一種截?cái)鄷r(shí)窗 的低通濾波多尺度全波形反演方法。
[0005] 本發(fā)明的解決方案是:利用巴特沃斯最平滑低通濾波器同時(shí)處理觀測(cè)記錄和正演 記錄，這樣有效避免了因數(shù)據(jù)不匹配而產(chǎn)生錯(cuò)誤的反演結(jié)果，該方法更適合實(shí)際地震數(shù)據(jù) 處理流程，可以得到精確的反演結(jié)果。針對(duì)時(shí)間域全波形反演計(jì)算效率低的問題，本發(fā)明提 出截?cái)鄷r(shí)窗全波形反演策略，即在剛開始反演的時(shí)候，只截取實(shí)際數(shù)據(jù)前部分時(shí)間進(jìn)行反 演，所以也只需要正演對(duì)應(yīng)的時(shí)窗數(shù)據(jù)進(jìn)行數(shù)據(jù)匹配，這樣通過正演截?cái)鄷r(shí)窗內(nèi)的波場(chǎng)來 代替正演全部時(shí)間波場(chǎng)，極大地縮短了正演所需要的時(shí)間。隨著反演的進(jìn)行，逐漸增加反演 時(shí)窗長(zhǎng)度，該截?cái)鄷r(shí)窗反演策略能夠減少每一步反演參數(shù)的個(gè)數(shù)，同時(shí)增加了程序的穩(wěn)定 性和反演結(jié)果的可靠性。
[0006] 本發(fā)明的目的是通過以下技術(shù)方案實(shí)現(xiàn)的：
[0007] 截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演方法的核心是在MATLAB 2013a平臺(tái)上完 成的。利用巴特沃斯濾波器對(duì)觀測(cè)數(shù)據(jù)和正演數(shù)據(jù)進(jìn)行低通濾波，通過設(shè)定巴特沃斯最平 滑濾波器不同的截?cái)囝l率，來實(shí)現(xiàn)時(shí)間域多尺度反演策略。低通濾波器能夠有效濾除地震 數(shù)據(jù)中的高頻成分，保留地震數(shù)據(jù)中的低頻成分，由于低頻成分含有地下構(gòu)造宏觀信息，有 利于初始模型的建立，同時(shí)避免周跳現(xiàn)象的發(fā)生。截?cái)鄷r(shí)窗反演策略，即選擇不同長(zhǎng)度的時(shí) 間窗口來提取觀測(cè)數(shù)據(jù)用于全波形反演，同時(shí)按照該時(shí)窗對(duì)應(yīng)的時(shí)間長(zhǎng)度進(jìn)行正演，這樣 可以避免每一步都正演所有時(shí)間數(shù)據(jù)，極大程度上節(jié)約了計(jì)算時(shí)間。
[0008] 截?cái)鄷r(shí)窗反演策略和低通濾波多尺度反演策略本身都具有壓制周跳和提高反演 精度的作用。為了得到更穩(wěn)定的反演能流程、更精確的反演結(jié)果和更高效的反演方法，本發(fā) 明是進(jìn)一步把二者結(jié)合起來使用(如圖1所示）。
[0009] 截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演方法，包括以下步驟：
[0010] a、安裝MATLAB2013a軟件，并安裝MATLAB語(yǔ)言編寫的地震數(shù)據(jù)處理軟件Crews工具 包；
[0011] b、對(duì)采集的地震數(shù)據(jù)進(jìn)行預(yù)處理，并將觀測(cè)記錄和觀測(cè)系統(tǒng)輸入計(jì)算機(jī)進(jìn)行時(shí)間 域全波形反演；
[0012] c、從觀測(cè)記錄中估計(jì)震源子波，利用估計(jì)得到的震源子波進(jìn)行正演(本發(fā)明利用 雷克子波代替實(shí)際數(shù)據(jù)中的震源子波）；
[0013] d、建立線性遞增初始模型，該速度模型作為反演的起始值，通過計(jì)算模型的更新 方向，實(shí)現(xiàn)由初始模型逐漸向真實(shí)模型靠近的過程；
[0014] e、根據(jù)最小二乘原理構(gòu)造目標(biāo)函數(shù)：
，其中V表示地下P波速 度，dobs為采集的觀測(cè)記錄，(1。31在速度模型通過正演得到的正演記錄。在求目標(biāo)函數(shù)的梯度 過程中對(duì)目標(biāo)函數(shù)兩端關(guān)于P波速度V求導(dǎo)，得到梯度表達(dá)式：
[0016] 其中Pb表示由伴隨震源fs向地下傳播得到的波場(chǎng)，Pf表示正演波場(chǎng)，s表示震源的 數(shù)目，r表示檢波器的數(shù)目。
[0017] f、選取截?cái)鄷r(shí)窗的時(shí)間長(zhǎng)度，截取觀測(cè)記錄前部分地震信號(hào)，并記錄截?cái)鄷r(shí)窗長(zhǎng) 度，該截?cái)鄷r(shí)窗的時(shí)間長(zhǎng)度作為正演的時(shí)間長(zhǎng)度；
[0018] g、按照聲波方程的動(dòng)力學(xué)規(guī)律，利用輸入的觀測(cè)系統(tǒng)、估計(jì)得到的震源子波和初 始速度模型進(jìn)行地震波場(chǎng)正演。該步驟只需要正演截?cái)鄷r(shí)窗內(nèi)的地震數(shù)據(jù)，不需要將所有 時(shí)間地震數(shù)據(jù)全部正演出來。存儲(chǔ)正傳波場(chǎng)數(shù)據(jù)和正演記錄。
[0019] 根據(jù)要求設(shè)定截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演方法相關(guān)參數(shù)，包括模型大 小nzXnx，網(wǎng)格距dx，dz，正演時(shí)窗長(zhǎng)度It，時(shí)間采樣間隔dt，雷克子波主頻fo,巴特沃斯濾 波器的截?cái)囝l率f?t，每個(gè)超隨機(jī)震源的最大迭代次數(shù)itermax，更新梯度最小值gtol，目標(biāo) 函數(shù)精度tol，反演結(jié)果的最大值vmax與最小值vmin;
[0020] h、對(duì)截?cái)鄷r(shí)窗內(nèi)的觀測(cè)記錄和正演記錄同時(shí)進(jìn)行低通濾波，確保用相同的低通濾 波器處理數(shù)據(jù)。濾波之后，得到處理觀測(cè)記錄和處理正演記錄；
[0021] i、將處理觀測(cè)記錄和處理正演記錄做差，得到殘差數(shù)據(jù)；
[0022] j、將殘差數(shù)據(jù)作為伴隨震源，用反傳算子作用于伴隨震源上，得到模型空間的殘 差反傳波場(chǎng)；
[0023] k、正傳波場(chǎng)與殘差反傳波場(chǎng)做互相關(guān)，得到模型更新梯度；
[0024] 1、用超記憶梯度優(yōu)化算法計(jì)算模型更新方向，并通過Wolfe收斂準(zhǔn)則尋找合適的 步長(zhǎng)；
[0025] m、判斷是夠滿足終止條件，若滿足則輸出截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演 方法反演結(jié)果。若不滿足終止條件，將反演結(jié)果作為下一個(gè)循環(huán)的初始模型，返回d步驟；
[0026] n、將上一步的輸出結(jié)果作為常規(guī)時(shí)間域全波形反演的初始速度模型，然后進(jìn)行常 規(guī)時(shí)間域全波形反演。迭代200次之后輸出最終反演結(jié)果。
[0027] 有益效果:通過對(duì)原有理論方法的改進(jìn)和相關(guān)地震全波形反演技術(shù)的整合，本發(fā) 明成功地將截?cái)鄷r(shí)窗和低通濾波等技術(shù)應(yīng)用到了全波形反演中。
[0028] 提出的截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演方法，該方法不僅能夠緩解因低頻 缺失帶來的周跳問題，而且極大程度的節(jié)約了正演的時(shí)間。利用截?cái)鄷r(shí)窗反演策略顯著提 高了全波形反演的計(jì)算效率和反演精度。利用巴特沃斯濾波器對(duì)觀測(cè)數(shù)據(jù)和正演數(shù)據(jù)進(jìn)行 低通濾波，濾除高頻成分，保留低頻成分，通過設(shè)定巴特沃斯最平滑濾波器不同的截?cái)囝l 率，來實(shí)現(xiàn)頻率由低到高的時(shí)間域多尺度反演策略，有效地壓制了全波形反演的周跳現(xiàn)象。 該反演策略一定程度上提高了全波形反演的精度，且更適合實(shí)際地震數(shù)據(jù)處理，同時(shí)本發(fā) 明提出的低通濾波反演策略無需對(duì)地震子波進(jìn)行濾波，有效避免因數(shù)據(jù)不匹配而造成的反 演錯(cuò)誤。在整個(gè)計(jì)算過程中震源都采用超隨機(jī)震源編碼策略，目的是在加快反演速率的同 時(shí)壓制超級(jí)炮內(nèi)部的串?dāng)_噪聲。與常規(guī)時(shí)間域全波形反演相比，本發(fā)明反演地下速度的過 程中用了截?cái)鄷r(shí)窗和低通濾波反演策略，解決了以下問題：
[0029] 1、首先將截?cái)鄷r(shí)窗和低通濾波策略組合起來，應(yīng)用到全波形反演中，可快速得到 一個(gè)高精度的初始速度模型，該初始速度模型基本恢復(fù)出了Marmousi模型的宏觀信息，利 用該初始速度模型進(jìn)行常規(guī)全波形反演，可以緩解周跳問題。
[0030] 2、將截?cái)鄷r(shí)窗函數(shù)引入到全波形反演中，使得每次正演不需要計(jì)算所有時(shí)刻的波 場(chǎng)，只計(jì)算截?cái)鄷r(shí)窗內(nèi)的波場(chǎng)即可，極大地縮短了正演所需要的計(jì)算時(shí)間。
[0031] 3、截?cái)鄷r(shí)窗低通濾波多尺度反演策略能有效減少每一次反演的未知數(shù)個(gè)數(shù)，根據(jù) 互相關(guān)求得的梯度和模型更新方向更加準(zhǔn)確。由于在模型更新迭代的過程中，需要通過 Wolfe線性搜索尋找合適的迭代步長(zhǎng)，當(dāng)梯度方向或者下降方向不準(zhǔn)確的情況下，Wolfe線 性搜索會(huì)浪費(fèi)大量的時(shí)間，而本發(fā)明提出的截?cái)鄷r(shí)窗低通濾波反演策略能較準(zhǔn)確的得到模 型更新的梯度方向，這樣在線性搜索的過程中節(jié)約大量的計(jì)算時(shí)間。
[0032] 4、常規(guī)全波形反演需要存儲(chǔ)正傳波場(chǎng)和殘差反傳波場(chǎng)，存儲(chǔ)這兩個(gè)波場(chǎng)需要大量 的計(jì)算機(jī)內(nèi)存。本發(fā)明提出的截?cái)鄷r(shí)窗反演策略，只需要存儲(chǔ)時(shí)窗內(nèi)的正傳波場(chǎng)和殘差反 傳波場(chǎng)即可，一定程度上節(jié)約了所需計(jì)算機(jī)的內(nèi)存空間，尤其在實(shí)際數(shù)據(jù)處理中。
[0033] 5、超隨機(jī)震源編碼策略能有效壓制超級(jí)炮內(nèi)部的串?dāng)_噪聲，該方法在保證反演精 度的同時(shí)一定程度上減少了正演的次數(shù)，有效地縮短了全波形反演的計(jì)算時(shí)間。
[0034] 模型測(cè)試結(jié)果表明，即使在缺失低頻數(shù)據(jù)和很強(qiáng)的高斯噪聲干擾的情況下，利用 截?cái)鄷r(shí)窗低通濾波策略的全波形反演也可以得到較好的反演結(jié)果。該反演策略既能實(shí)現(xiàn)在 低頻段由淺層到深層的逐層反演方法，又能在計(jì)算效率方面得到顯著的改善，同時(shí)緩解了 因缺失低頻數(shù)據(jù)而導(dǎo)致的周跳現(xiàn)象。可以看出本發(fā)明提出的截?cái)鄷r(shí)窗低通濾波反演策略可 為常規(guī)全波形反演提供一個(gè)高精度的初始速度模型。
【附圖說明】
[0035] 圖1截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演方法流程圖。
[0036]圖2雷克子波
[0037](左)雷克子波波形圖，（右)雷克子波頻譜圖。
[0038]圖3巴特沃斯最平滑濾波器圖。
[0039]圖4速度模型
[0040] (左)線性遞增初始速度模型圖，（右)真實(shí)速度模型圖。
[0041] 圖5濾波波形圖和頻譜圖
[0042] (a) 15Hz低通濾波波形圖和真實(shí)地震記錄波形圖，
[0043] (b)25Hz低通濾波波形圖和真實(shí)地震記錄波形圖，
[0044] (c)15Hz低通濾波，25Hz低通濾波和真實(shí)地震記錄頻譜圖，
[0045] (d)lOHz高通濾波波形圖，
[0046] (e)lOHz高通濾波頻譜圖。
[0047]圖6超隨機(jī)震源觀測(cè)記錄 [0048](左)無噪數(shù)據(jù)，（右)含噪數(shù)據(jù)。
[0049] 圖7截?cái)鄷r(shí)窗示意圖。
[0050] 圖8觀測(cè)記錄低通濾波 [0051](左）15Hz低通濾波觀測(cè)記錄，
[0052](右)25Hz低通濾波觀測(cè)記錄。
[0053]圖9常規(guī)時(shí)間域全波形反演結(jié)果圖。
[0054]圖10低通濾波多尺度全波形反演結(jié)果圖 [0055] (a) 15Hz低通濾波全波形反演結(jié)果圖，
[0056] (b)25Hz低通濾波全波形反演結(jié)果圖，
[0057] (c)低通濾波最終反演結(jié)果圖。
[0058]圖11截?cái)鄷r(shí)窗全波形反演結(jié)果圖。
[0059] 圖12截?cái)鄷r(shí)窗低通濾波多尺度全波形反演結(jié)果圖。
[0060] 圖13常規(guī)時(shí)間域全波形反演結(jié)果圖（缺失低頻測(cè)試）。
[0061 ]圖14低通濾波多尺度全波形反演結(jié)果圖（缺失低頻測(cè)試）。
[0062] (a) 15Hz低通濾波全波形反演結(jié)果圖，
[0063] (b)25Hz低通濾波全波形反演結(jié)果圖，
[0064] (c)低通濾波最終反演結(jié)果圖。
[0065] 圖15截?cái)鄷r(shí)窗全波形反演結(jié)果圖（缺失低頻測(cè)試）。
[0066] 圖16截?cái)鄷r(shí)窗低通濾波多尺度全波形反演結(jié)果圖（缺失低頻測(cè)試）。
[0067] (a)截?cái)鄷r(shí)窗+15Hz低通濾波全波形反演結(jié)果圖，
[0068] (b)25Hz低通濾波全波形反演結(jié)果圖，
[0069] (c)截?cái)鄷r(shí)窗低通濾波最終反演結(jié)果圖。
[0070]圖17常規(guī)時(shí)間域全波形反演結(jié)果圖（含高斯噪聲測(cè)試）。
[0071 ]圖18截?cái)鄷r(shí)窗全波形反演結(jié)果圖（含高斯噪聲測(cè)試）。
[0072] 圖19低通濾波多尺度全波形反演結(jié)果圖（含高斯噪聲測(cè)試）。
[0073] 圖20截?cái)鄷r(shí)窗低通濾波多尺度全波形反演結(jié)果圖（含高斯噪聲測(cè)試）。
【具體實(shí)施方式】
[0074]下面結(jié)合附圖和實(shí)例對(duì)本發(fā)明進(jìn)一步的詳細(xì)說明
[0075]截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演方法，包括以下步驟：
[0076] a、程序是在MATLAB2013a軟件框架下編寫完成，首先需要安裝MATLAB2013a軟件， 并安裝MATLAB語(yǔ)言編寫的地震數(shù)據(jù)處理軟件Crews工具包；
[0077] b、采集的地震數(shù)據(jù)首先需進(jìn)行預(yù)處理，并將觀測(cè)記錄和觀測(cè)系統(tǒng)輸入計(jì)算機(jī)進(jìn)行 時(shí)間域全波形反演。其中預(yù)處理包括：
[0078] B1、多次波衰減、消除交混回響和壓制虛反射等不能正演的地震波。
[0079] B2、低頻保護(hù)去噪、缺失地震道補(bǔ)償、保幅處理等。
[0080] B3、定義觀測(cè)系統(tǒng)。
[0081] c、從觀測(cè)記錄中估計(jì)震源子波，該震源子波用于正演(本發(fā)明利用雷克子波代替 實(shí)際數(shù)據(jù)中的震源子波(圖2));
[0082] d、建立線性遞增初始模型（圖4左），該速度模型作為反演的起始值，通過計(jì)算模型 的更新方向，來實(shí)現(xiàn)由初始模型逐漸向真實(shí)模型靠近的過程；
[0083] e、根據(jù)最小二乘原理構(gòu)造目標(biāo)函數(shù)：
，其中v表示地下P波速 度，cUs為采集的觀測(cè)記錄，(1。31在速度模型通過正演得到的正演記錄。在求目標(biāo)函數(shù)的梯度 過程中對(duì)目標(biāo)函數(shù)兩端關(guān)于P波速度v求導(dǎo)，得到梯度表達(dá)式：
[0085]其中Pb表示由伴隨震源fs向地下傳播得到的波場(chǎng)，pf表示正演波場(chǎng)，s表示震源的 數(shù)目，r表示檢波器的數(shù)目。
[0086] f、選取截?cái)鄷r(shí)窗（圖7)，截取觀測(cè)記錄前部分地震信號(hào)，并記錄截?cái)鄷r(shí)窗長(zhǎng)度，該 截?cái)鄷r(shí)窗的時(shí)間長(zhǎng)度作為正演的時(shí)間長(zhǎng)度；
[0087] g、按照聲波方程的動(dòng)力學(xué)規(guī)律，利用輸入的觀測(cè)系統(tǒng)、估計(jì)得到的震源子波和初 始速度模型進(jìn)行地震波場(chǎng)正演。該步驟只需要正演截?cái)鄷r(shí)窗內(nèi)的地震數(shù)據(jù)，不需要將所有 時(shí)間地震數(shù)據(jù)全部正演出來。存儲(chǔ)正傳波場(chǎng)數(shù)據(jù)和正演記錄。
[0088] 根據(jù)要求設(shè)定截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演方法相關(guān)參數(shù)，包括模型大 小nzXnx，網(wǎng)格距dx，dz，正演時(shí)窗長(zhǎng)度It，時(shí)間采樣間隔dt，雷克子波主頻fo,巴特沃斯濾 波器的截?cái)囝l率(f?t)，每個(gè)超隨機(jī)震源的最大迭代次數(shù)itermax，更新梯度最小值gtol，目 標(biāo)函數(shù)精度tol，反演結(jié)果的最大值vmax和最小值vmin;
[0089] 超隨機(jī)震源編碼:時(shí)間域全波形反演在單炮正演的過程中耗費(fèi)大量的時(shí)間，有人 提出利用同時(shí)激發(fā)多個(gè)震源(超級(jí)炮)來代替單炮激發(fā)震源，這樣能縮短正演所需要的時(shí) 間，但是利用同時(shí)激發(fā)震源進(jìn)行全波形反演又給反演結(jié)果帶來了串?dāng)_噪聲。為了提高全波 形反演的計(jì)算效率同時(shí)壓制串?dāng)_噪聲本軟件利用超隨機(jī)震源編碼策略。
[0090] 隨機(jī)相位編碼地表超級(jí)炮激發(fā)出的雷克子波，要求相位正負(fù)隨機(jī)分配，同時(shí)每個(gè) 超級(jí)炮中的單炮激發(fā)延時(shí)不同，延時(shí)最大不超過400ms。設(shè)隨機(jī)相位編碼信息為：
[0092]其中s'US={1，2…Ns}，s中元素的個(gè)數(shù)固定，其表示一個(gè)超級(jí)炮中所含單炮數(shù)目 個(gè)數(shù)。超級(jí)炮中單炮數(shù)量可根據(jù)實(shí)際需求調(diào)整混合比例。s中元素大小隨機(jī)分配。iGN+表示 隨機(jī)的正整數(shù)。超隨機(jī)震源編碼可以表示為：
[0093] spf = Dy ? A ? r ? T ? L ? F ? f (4)
[0094]其中Dy表示動(dòng)態(tài)震源編碼函數(shù)(動(dòng)態(tài)編碼的意思就是在超記憶梯度法每迭代10次 之后重新更換隨機(jī)相位編碼信息，最終可實(shí)現(xiàn)地表全覆蓋。該方法能夠有效地增加地表炮 的覆蓋密度，同時(shí)減弱串?dāng)_噪聲對(duì)全波形反演的影響，最終實(shí)現(xiàn)計(jì)算效率上的巨大飛躍。）， A表示隨機(jī)振幅編碼函數(shù)，r表示隨機(jī)相位編碼函數(shù)，T表示隨機(jī)延時(shí)編碼函數(shù)，L表示隨機(jī) 位置編碼函數(shù)，F(xiàn)表示隨機(jī)雷克子波主頻編碼函數(shù)，f震源函數(shù)。
[0095] h、對(duì)截?cái)鄷r(shí)窗內(nèi)的觀測(cè)記錄和正演記錄同時(shí)進(jìn)行低通濾波，確保用相同的低通濾 波器(圖3)處理數(shù)據(jù)。濾波之后，得到處理觀測(cè)記錄和處理正演記錄；
[0096]巴特沃斯低通濾波器：以巴特沃斯函數(shù)來近似濾波器的系統(tǒng)函數(shù)。巴特沃斯濾波 器是根據(jù)幅頻特性在通頻帶內(nèi)具有最平滑特性定義的濾波器。
[0097]巴特沃思濾波器的低通模平方函數(shù)表示
[0099]巴特沃斯濾波器的主要特征如下：
[0102]把、|凡(^)|2是如勺單調(diào)下降函數(shù)；
[0103] H4、|Ha( j Q ) |2隨著階次N的增大而更接近于理想低通濾波器；
[0104] 濾波器的幅頻特性隨著濾波器階次N的增加而變得越來越好，在截止頻率n。處的 函數(shù)值始終為1/2的情況下，通帶內(nèi)有更多的頻帶區(qū)的值接近于1;在阻帶內(nèi)更迅速的趨近 于零。
[0105] i、將處理觀測(cè)記錄和處理正演記錄做差，得到殘差數(shù)據(jù)；
[0106] j、將殘差數(shù)據(jù)作為伴隨震源，用反傳算子作用于伴隨震源上，得到模型空間的殘 差反傳波場(chǎng)；
[0107] k、正傳波場(chǎng)與殘差反傳波場(chǎng)做互相關(guān)得到模型更新梯度；
[0108] 全波形反演的目標(biāo)函數(shù)一般定義為正演數(shù)據(jù)與觀測(cè)數(shù)據(jù)殘差的L2泛函，形式如 下：
[0110] 其中0。31便是正演記錄，Dcal表示實(shí)際野外的觀測(cè)記錄。將目標(biāo)函數(shù)對(duì)速度v求導(dǎo)得 到：
[0112]本發(fā)明為了更清晰的闡述理論方法，假設(shè)密度為常數(shù)，則只需要對(duì)速度進(jìn)行反演。 常密度聲波方程為：
[0114]其中u為聲波波場(chǎng)，f為震源函數(shù)。v表示地下介質(zhì)的聲速度。
[0115] 將聲波方程簡(jiǎn)化寫成算子的形式：
[0116] Lu = s (9)
[0118] 方程(9)兩邊同時(shí)對(duì)速度v求導(dǎo)得到(Pratt et al.，1998):
[0120]得到波場(chǎng)對(duì)地下速度v的導(dǎo)數(shù)為：
[0122]將方程(11)代入方程(7)中得到：
[0124] 其中I/1表示伴隨算子或稱為反傳算子，表示殘差反傳波場(chǎng)。fs = Dcal-D〇bs稱為伴隨震源。
[0125] 1、用超記憶梯度優(yōu)化算法計(jì)算模型更新方向，并通過Wolfe收斂準(zhǔn)則尋找合適的 步長(zhǎng)。
[0126] 超記憶梯度法：用當(dāng)前梯度及之前多步梯度的信息來對(duì)當(dāng)前梯度的信息進(jìn)行校 正，以加快目標(biāo)函數(shù)的收斂速度。其迭代形式為：
[0127] vk+i = Vk+akdk (13)
[0128] 其中Vk代表迭代第k部速度模型，ak表示步長(zhǎng)，dk表示下降方向。
[0129] 超記憶梯度法下降方向可以通過如下形式來表示：
[0133] 其中，0<P<l，mem>0是一個(gè)正整數(shù)，mem我們稱之為記憶度。不同記憶度目標(biāo)函 數(shù)收斂速度及模型反演精度不同，記憶梯度數(shù)量過多則會(huì)導(dǎo)致梯度方向校正過度，目標(biāo)函 數(shù)下降變慢，記憶梯度數(shù)量少則退化為共輒梯度法。經(jīng)過實(shí)驗(yàn)驗(yàn)證得到，在地震全波形反演 中最好的記憶度大致在mem=4，5。超記憶梯度法有著一定的優(yōu)勢(shì)，但當(dāng)梯度記憶過多，計(jì)算 速度變慢，內(nèi)存需求增大。所以選擇合適的記憶度才能發(fā)揮超記憶梯度法最大的優(yōu)勢(shì)。
[0134] m、判斷是夠滿足終止條件，若滿足則輸出截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演 方法反演結(jié)果。若不滿足終止條件，將反演結(jié)果作為下一個(gè)循環(huán)的初始模型，返回d步驟。
[0135] n、將上一步的輸出結(jié)果作為常規(guī)時(shí)間域全波形反演的初始速度模型，然后進(jìn)行常 規(guī)時(shí)間域全波形反演。迭代200次之后輸出最終反演結(jié)果。
[0136] 實(shí)施例1
[0137] 根據(jù)勘探要求，將MATLAB2013a在Windows 7旗艦版系統(tǒng)下進(jìn)行安裝，并安裝 MATLAB語(yǔ)言編寫的地震數(shù)據(jù)處理軟件Crews工具包。
[0138] 本發(fā)明利用Marmousi P波速度模型進(jìn)行測(cè)試，Marmousi速度模型具有復(fù)雜的構(gòu) 造、細(xì)微的層理和深部?jī)?chǔ)油藏構(gòu)造，很適合進(jìn)行理論方法測(cè)試。原始Marmousi速度模型較 大，考慮到硬件內(nèi)存以及計(jì)算時(shí)間因素，本發(fā)明對(duì)原始模型進(jìn)行抽稀處理，對(duì)抽稀之后的 Marmousi速度模型進(jìn)行截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演方法測(cè)試。真實(shí)模型（圖4 右)和線性遞增初始模型（圖4左）。
[0139] 模型參數(shù)如下：
[0140] 表1截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演方法測(cè)試參數(shù)(時(shí)間域）
[0142]截?cái)鄷r(shí)窗反演策略參數(shù)如下：
[0143] 對(duì)觀測(cè)記錄利用截?cái)鄷r(shí)窗進(jìn)行處理，設(shè)定截?cái)鄷r(shí)窗長(zhǎng)度It = 0.2:0.1: 2. Is，即選 擇起始截?cái)鄷r(shí)窗〇-〇. 2s，然后每迭代10次截?cái)鄷r(shí)窗的時(shí)間長(zhǎng)度增加0.1 s，且起始0時(shí)刻不變 (圖7)，直到截?cái)鄷r(shí)窗的總長(zhǎng)度為2. Is。該過程共20個(gè)截?cái)鄷r(shí)窗，全波形反演共迭代200次。 正演只需要正演對(duì)應(yīng)時(shí)窗的時(shí)間長(zhǎng)度，這樣能很大程度上節(jié)約正演所需要的計(jì)算時(shí)間。截 斷時(shí)窗反演策略由淺層到深層，逐層進(jìn)行反演，極大地減小了每一次反演的參數(shù)，有效避免 了周跳現(xiàn)象的發(fā)生，同時(shí)增強(qiáng)了程序的穩(wěn)定性。
[0144] 低通濾波多尺度反演策略參數(shù)如下：
[0145] 為了克服全波形反演過程中陷入局部極小和周跳等問題，同時(shí)也為了提高成像質(zhì) 量，用低通濾波多尺度方法進(jìn)行全波形反演，設(shè)定15Hz和25Hz兩個(gè)截?cái)囝l率。兩個(gè)截?cái)囝l率 將整個(gè)全波形反演過程分成三個(gè)頻段：
[0146] I、0-15Hz:低頻段。在低頻段進(jìn)行反演可以得到地下模型的宏觀信息，在低頻段的 反演結(jié)果作為下一個(gè)頻率段反演的初始模型，這樣可以逐漸把速度模型精確地反演出來。 但是實(shí)際數(shù)據(jù)中往往缺失低頻信息，陸地地震數(shù)據(jù)缺失10Hz以下的數(shù)據(jù)，低通濾波只是在 存在低頻數(shù)據(jù)的情況下才能發(fā)揮作用，當(dāng)缺失低頻數(shù)據(jù)時(shí)低通濾波多尺度全波形反演的效 果并不是很好(圖13)。所以本發(fā)明在低頻段內(nèi)結(jié)合截?cái)鄷r(shí)窗反演策略來消除因缺失低頻而 導(dǎo)致的周跳現(xiàn)象。
[0147] n、0-25Hz:中頻段。該頻段主要起到一個(gè)中間過度的作用，低頻段的反演結(jié)果作 為中頻段的初始速度模型。
[0148] m、〇-end:全頻段。中頻帶的反演結(jié)果全頻帶的初始速度模型，該步驟即可以得到 時(shí)間域全波形的反演結(jié)果。
[0149] 基于以上參數(shù)，在表2所示的測(cè)試環(huán)境下進(jìn)行全波形反演。
[0150] 表2計(jì)算機(jī)性能測(cè)試環(huán)境
[0152] 表3全波形反演方法對(duì)比
[0153]
[0154] 從表3可以得出，常規(guī)時(shí)間域全波形反演與本發(fā)明提出的截?cái)鄷r(shí)窗的低通濾波多 尺度全波形反演具有明顯的差別(模型誤差，計(jì)算時(shí)間方面）。常規(guī)時(shí)間域全波形反演的計(jì) 算時(shí)間為101.3分鐘，而截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演只需43.6分鐘就可以完成 整個(gè)計(jì)算過程。從模型誤差的角度結(jié)合反演結(jié)果圖（圖9，10，11，12)可以看出截?cái)鄷r(shí)窗的低 通濾波多尺度全波形反演反演結(jié)果精度更高，反演得到的速度結(jié)果更加準(zhǔn)確，而常規(guī)時(shí)間 域全波形反演的結(jié)果甚至出現(xiàn)了錯(cuò)誤，與真實(shí)模型相差很遠(yuǎn)。
[0155] 基于以上參數(shù)，在表2所示的測(cè)試環(huán)境下進(jìn)行多種影響因素測(cè)試：
[0156] 表4不同影響因素下截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演方法測(cè)試結(jié)果
[0158] 注:表3和表4中模型反演歸一化誤差計(jì)算公式為
，其中vtrue代 表真實(shí)速度模型，vk代表反演結(jié)果。
[0159] 從表4和圖13,14,15中可以看出，缺失低頻對(duì)全波形反演有著嚴(yán)重的影響。當(dāng)?shù)卣?數(shù)據(jù)中缺失低頻的時(shí)候，常規(guī)時(shí)間域全波形反演結(jié)果（圖13)，低通濾波多尺度全波形反演 結(jié)果（圖14)和截?cái)鄷r(shí)窗全波形反演結(jié)果（圖15)都出現(xiàn)了明顯的錯(cuò)誤，模型反演誤差明顯偏 高。同時(shí)，由于缺失低頻的因素導(dǎo)致模型更新梯度計(jì)算結(jié)果不正確，使得在進(jìn)行Wolf e線性 搜索步長(zhǎng)的過程中浪費(fèi)了大量的時(shí)間，這就解釋了為何相同的迭代次數(shù)，當(dāng)缺失低頻的時(shí) 候全波形反演的計(jì)算時(shí)間要明顯長(zhǎng)于完整頻帶全波形反演的計(jì)算時(shí)間。為了更好地解決因 缺失低頻成分帶來的反演錯(cuò)誤問題，本發(fā)明將截?cái)鄷r(shí)窗反演策略和低通濾波多尺度策略結(jié) 合使用，從圖16可以看出，該結(jié)合策略使得全波形反演的結(jié)果得到了明顯的改善。模型誤差 也說明了，該策略得到的反演結(jié)果與真實(shí)模型更加接近。
[0160]從表4和圖17,18,19,20中可以看出，當(dāng)?shù)卣饠?shù)據(jù)中含有高斯噪聲的時(shí)候，反演結(jié) 果變得不清晰了。常規(guī)時(shí)間域全波形反演（圖17)和截?cái)鄷r(shí)窗全波形反演（圖18)都出現(xiàn)了明 顯的錯(cuò)誤，說明截?cái)鄷r(shí)窗反演策略雖然在一定程度上加快了反演進(jìn)程，但是其抗噪能力較 弱。低通濾波多尺度全波形反演（圖19)較前兩者具有明顯的優(yōu)勢(shì)，證明了低通濾波多尺度 策略具有較強(qiáng)的抗噪能力。將截?cái)鄷r(shí)窗反演策略和低通濾波多尺度反演策略結(jié)合起來使用 (圖20)，不僅加快了反演的進(jìn)程，而且彌補(bǔ)了截?cái)鄷r(shí)窗抗噪能力弱的缺陷。
[0161]圖1是整個(gè)反演過程的流程圖，從流程圖可以看出首先在低頻段利用截?cái)鄷r(shí)窗反 演策略進(jìn)行全波形反演，然后再把低頻段截?cái)鄷r(shí)窗反演結(jié)果作為初始模型用于下一步中頻 段反演，最后將中頻段的反演結(jié)果最為常規(guī)時(shí)間域全波形反演的初始模型。該截?cái)鄷r(shí)窗低 通濾波多尺度全波形反演策略，在反演精度、計(jì)算效率、低頻依賴、抗噪能力等方面優(yōu)于常 規(guī)時(shí)間域全波形反演。
【主權(quán)項(xiàng)】
1. 一種截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演方法，其特征在于，包括以下步驟： a、 安裝MATLAB2013a軟件，并安裝MATLAB語(yǔ)言編寫的地震數(shù)據(jù)處理軟件Crews工具包； b、 對(duì)采集的地震數(shù)據(jù)進(jìn)行預(yù)處理，并將觀測(cè)記錄和觀測(cè)系統(tǒng)輸入計(jì)算機(jī)，進(jìn)行時(shí)間域 全波形反演； c、 從觀測(cè)記錄中估計(jì)震源子波，該估計(jì)得到的震源子波用于正演模擬； d、 建立線性遞增初始模型，該速度模型作為反演的起始值，通過計(jì)算模型的更新方向， 實(shí)現(xiàn)由初始模型逐漸向真實(shí)模型逼近的過程； e、 根據(jù)最小二乘法原理構(gòu)建目標(biāo)函數(shù)其中v表示P波速度，(Us為采集的觀測(cè)記錄，dw在速度模型通過正演得到的正演記錄， 在求目標(biāo)函數(shù)的梯度過程中，對(duì)目標(biāo)函數(shù)兩端的P波速度v求導(dǎo)，得到梯度表達(dá)式：其中pb表示由伴隨震源fs向地下傳播得到的波場(chǎng)，pf表示正演模擬波場(chǎng)，s表示震源的 數(shù)目，r表示檢波器的數(shù)目； f、 選取截?cái)鄷r(shí)窗的時(shí)間長(zhǎng)度，截取觀測(cè)記錄前部分地震信號(hào)，并記錄截?cái)鄷r(shí)窗長(zhǎng)度，該 截?cái)鄷r(shí)窗的時(shí)間長(zhǎng)度作為正演模擬的時(shí)間長(zhǎng)度； g、 按照聲波方程的動(dòng)力學(xué)規(guī)律，利用輸入的觀測(cè)系統(tǒng)、估計(jì)得到的震源子波和初始速 度模型進(jìn)行地震波場(chǎng)正演，存儲(chǔ)正傳波場(chǎng)數(shù)據(jù)和正演記錄； 根據(jù)要求設(shè)定截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演方法相關(guān)參數(shù)，包括模型大小nz Xnx，網(wǎng)格距dx，dz，正演模擬時(shí)窗長(zhǎng)度It，時(shí)間采樣間隔dt，雷克子波主頻fo,巴特沃斯濾 波器的截?cái)囝l率f?t，每個(gè)超隨機(jī)震源的最大迭代次數(shù)iter max，更新梯度最小值gtol，目 標(biāo)函數(shù)精度tol，反演結(jié)果的最大值v max和最小值v min; h、 對(duì)截?cái)鄷r(shí)窗內(nèi)的觀測(cè)記錄和正演記錄同時(shí)進(jìn)行低通濾波，確保用相同的低通濾波器 處理數(shù)據(jù)，濾波之后，得到處理觀測(cè)記錄和處理正演記錄； i、 將處理觀測(cè)記錄和處理正演記錄做差，得到殘差數(shù)據(jù)； j、 將殘差數(shù)據(jù)作為伴隨震源，用反傳算子作用于伴隨震源上，得到模型空間的殘差反 傳波場(chǎng)； k、 正傳波場(chǎng)與反傳殘差波場(chǎng)做互相關(guān)，得到模型更新梯度； l、 用超記憶梯度優(yōu)化算法計(jì)算模型更新方向，并通過Wolfe收斂準(zhǔn)則尋找合適的步長(zhǎng)； m、 判斷是否滿足終止條件，若滿足則輸出截?cái)鄷r(shí)窗的低通濾波多尺度全波形反演方法 的反演結(jié)果;若不滿足終止條件，將反演結(jié)果作為下一個(gè)循環(huán)的初始模型，返回d步驟； n、 將m步驟輸出結(jié)果作為常規(guī)時(shí)間域全波形反演的初始速度模型，然后進(jìn)行常規(guī)時(shí)間 域全波形反演，迭代200次之后輸出最終反演結(jié)果。
【文檔編號(hào)】G01V1/28GK106054244SQ201610429648
【公開日】2016年10月26日
【申請(qǐng)日】2016年6月16日
【發(fā)明人】胡勇, 韓立國(guó), 張盼, 鞏向博, 張鳳蛟
【申請(qǐng)人】吉林大學(xué)