鬼佬大哥大
  • / 23
  • 下載費用:30 金幣  

一種反褶積方法及裝置.pdf

摘要
申請專利號:

CN201410132417.3

申請日:

2014.04.03

公開號:

CN103954992A

公開日:

2014.07.30

當前法律狀態:

授權

有效性:

有權

法律詳情: 授權|||實質審查的生效IPC(主分類):G01V 1/28申請日:20140403|||公開
IPC分類號: G01V1/28 主分類號: G01V1/28
申請人: 中國石油天然氣股份有限公司
發明人: 王萬里; 魏新建; 何欣; 祿娟
地址: 100007 北京市東城區東直門北大街9號
優先權:
專利代理機構: 北京三友知識產權代理有限公司 11127 代理人: 黨曉林
PDF完整版下載: PDF下載
法律狀態
申請(專利)號:

CN201410132417.3

授權公告號:

||||||

法律狀態公告日:

2017.04.05|||2014.08.27|||2014.07.30

法律狀態類型:

授權|||實質審查的生效|||公開

摘要

本發明提供了一種反褶積方法及裝置。該方法包括以下處理步驟:將采集的地震數據處理成疊后地震數據記錄;構建包括反射系數序列在內的目標函數;在所述疊后地震數據記錄中選取一道所需處理的地震資料,從所述構建的目標函數中計算出所述地震資料中反射系數序列的振幅和位置;讀取下一道所需處理的地震資料,計算出所述下一道地震資料反射系數序列的振幅和位置;計算出所述疊后地震數據記錄中所有地震道反射系數序列的振幅和位置,輸出反射系數序列,完成所述疊后地震數據記錄的反褶積處理。本發明提供的反褶積方法能進一步提高地震資料的分辨率。

權利要求書

權利要求書
1.  一種反褶積方法,其特征在于,包括以下處理步驟:
S1:將所需處理的地震數據處理成疊后地震數據記錄;
S2:構建包括反射系數序列在內的目標函數;
S3:在所述疊后地震數據記錄中選取一道所需處理的地震資料,從所述構建的目標函數中計算出所述地震資料的反射系數序列的振幅和位置;
S4:讀取下一道所需處理的地震資料,計算出所述下一道地震資料反射系數序列的振幅和位置;
S5:計算出所述疊后地震數據記錄中所有地震道反射系數序列的振幅和位置,輸出反射系數序列,完成所述疊后地震數據記錄的反褶積處理。

2.  如權利要求1所述的一種反褶積方法,其特征在于,所述構建目標函數包括:
S201:建立地震數據褶積模型;
S202:將所述地震數據褶積模型轉換成矩陣形式;
S203:由所述矩陣形式構建目標函數;所述構建的目標函數為:
Jα=||Aα-S||2
上式中,Jα為關于反射系數序列α的范數,A為一個N×M的矩陣,其表達式為其中,ωt為地震子波,t=1,2,…N,N為所述疊后地震數據記錄的時間窗長度,M為預先設置的所述時間窗長度N內反射系數序列的個數,τM第M個反射系數的時延,s為t時刻所述疊后地震數據記錄的值。

3.  如權利要求2所述的一種反褶積方法,其特征在于,所述構建的目標函數為:
J=||Aα-S||2+μ||Cα-ε||2=Jα+μJc
上式中,μ為阻抗權系數,C為左下三角全為1的矩陣,ε為波阻抗。

4.  如權利要求3所述的一種反褶積方法,其特征在于,所述構建的目標函數為:
J=Jα+μJc+βJα
上式中,β為正則化的阻尼系數,其中,β=β0f0,β0為預設值,f0=max{(ATA+μCTCii)},i為矩陣(ATA+μCTC)的下標。

5.  如權利要求1-4任意一項所述的一種反褶積方法,其特征在于,從所述構建的目標函數中計算出所述地震資料的反射系數序列的振幅包括步驟:
S310:從所述疊后地震數據記錄中估計出地震子波;
S320:將所述估計出的地震子波代入所述目標函數中,計算出所述反射系數序列的振幅;
其中,所述估計地震子波的方法包括步驟:
S311:建立所述疊后地震數據記錄的褶積模型;
S312:將所述建立的褶積模型變換到復賽譜域;
S313:設計低通濾波器,在復賽譜域將所述建立的褶積模型的地震子波和反射系數序列進行分離,得到復賽譜域的地震子波;
S314:將所述復賽譜域的地震子波轉換到時間域,得到時間域的地震子波。

6.  如權利要求5所述的一種反褶積方法,其特征在于,所述地震子波的估計方法還包括步驟:
315:對所述估計出的地震子波相位進行校正。

7.  如權利要求5所述的一種反褶積方法,其特征在于,采用模擬退火算法計算反射系數序列中每一個反射系數的時間位置,其操作過程包括以下步驟:
S301:在所述時間窗長度N內隨機給定反射系數序列a(j)的初始時間位置,由所述初始時間位置計算所述目標函數Ja當前解,j=1,2…M;
S302:在所述時間窗長度N內對所述反射系數序列a(j)進行擾動,將所述擾動后的反射系數序列a(j)代入所述目標函數Ja中,計算出所述目標函數Ja的新解;
S303:將所述Jα的新解與當前解進行比較,根據Metropolis準則判斷是否接受所述目標函數Jα的新解;
如果判斷結果為接受所述新解,則接受反射系數α(j)當前的時間位置;此時,所述新解作為目標函數Ja當前解;
S304:重復S302、S303直到達到收斂條件后,停止迭代,得到所述反射系數α(j)的時間位置。

8.  如權利要求7所述的一種反褶積方法,其特征在于,所述模擬退火算法還包括在所述S303之后:
S3031:設置誤差項σ,將所述新解與當前解的差值與所述誤差項σ比較;如果所述新解與當前解的差值大于所述誤差項σ,則不接受所述目標函數Jα的新解。

9.  如權利要求7所述的一種反褶積方法,其特征在于,在達到收斂條件后,將所述模擬退火算法得出的反射系數序列α(j)的振幅值作為地震資料中反射系數序列的振幅。

10.  一種反褶積裝置,其特征在于,該裝置包括:數據預處理模塊、目標函數建立模塊、目標函數計算模塊、資料讀取模塊、反射系數序列輸出模塊,其中:
所述數據預處理模塊,用于將所需處理的地震數據處理成疊后地震數據記錄;
所述目標函數建立模塊,用于構建包括反射系數序列在內的目標函數;
所述目標函數計算模塊,用于在所述疊后地震數據記錄中選取一道所需處理的地震資料,從所述構建的目標函數中計算出所述地震資料的反射系數序列的振幅和位置;
所述資料讀取模塊,用于讀取下一道所需處理的地震資料,計算出所述下一道地震資料反射系數序列的振幅和位置;
所述反射系數序列輸出模塊,用于計算出所述疊后地震數據記錄中所有地震道反射系數序列的振幅和位置,輸出反射系數序列,完成所述疊后地震數據記錄的反褶積處理。

11.  如權利要求10所述的一種反褶積裝置,其特征在于,所述目標函數建立模塊構建的目標函數為:
Jα=||Aα-S||2
上式中,Jα為關于反射系數序列α的范數形式的函數,A為一個N×M的矩陣,其表達式為其中,ωt為地震子波,t=1,2,…N,N為所述疊后地震數據記錄的時間窗長度,M為預先設置的所述時間窗長度N內反射系數序列的個數,τM第M個反射系數的時延,S為t時刻所述疊后地震數據記錄的值。

12.  如權利要求11所述的一種反褶積裝置,其特征在于,所述目標函數建立模塊構建的目標函數為:
J=||Aα-S||2+μ||Cα-ε||2=Jα+μJc
上式中,μ為阻抗權系數,C為左下三角全為1的矩陣,ε為波阻抗。

13.  如權利要求12所述的一種反褶積裝置,其特征在于,所述目標函數建立模塊構建的目標函數為:
J=Jα+μJc+βJα
上式中,β為正則化的阻尼系數,其中,β=β0f0,β0為預設值,f0=max{(ATA+μCTCii)},i為矩陣(ATA+μCTC)的下標。

14.  如權利要求11所述的一種反褶積裝置,其特征在于,所述目標函數計算模塊包括振幅計算模塊、位置計算模塊,其中:
所述振幅計算模塊,用于從所述目標函數建立模塊中計算出所述地震資料中反射系數序列的振幅;
所述位置計算模塊,用于從所述目標函數建立模塊中計算出所述地震資料中反射系數序列的時間位置。

15.  如權利要求14所述的一種反褶積裝置,其特征在于,位置計算模塊包括初始模塊、擾動模塊、迭代模塊、收斂模塊,其中:
所述初始模塊1,用于在所述時間窗長度N內隨機給定反射系數序列a(j)的初始時間位置,由所述初始時間位置計算所述目標函數Ja當前解,j=1,2…M;
在所述擾動模塊,用于在所述時間窗長度N內對所述反射系數序列a(j)進行擾動,將所述擾動后的反射系數序列a(j)代入所述目標函數Ja中,計算出所述目標函數Ja的新解;
所述迭代模塊,用于將所述Jα的新解與當前解進行比較,根據Metropolis準則判斷是否接受所述目標函數Jα的新解;
如果判斷結果為接受所述新解,則接受反射系數α(j)當前的時間位置;此時,所述新解作為目標函數Ja當前解;
所述收斂模塊,用于繼續計算所述新解,并與所述當前解進行比較直到達到收斂條件后,停止迭代,得到所述反射系數α(j)的時間位置。

說明書

說明書一種反褶積方法及裝置
技術領域
本方法涉及地震勘探中地震資料處理領域,特別涉及一種反褶積方法及裝置。
背景技術
在對地層蘊藏的石油天然氣資源進行開發之前,需要對地層進行勘測以測定地下的構造,為油氣資源的開采和地層識別提供依據。地震勘探是勘測地層結構的重要手段。
地震勘探是一種利用人工地震技術探測地下構造的勘探方法。它以人工方法按照一定的方式在地表附近激發地震波,產生稱之為地震子波的振動信號。來自不同反射深度界面的地震子波以不同的時間到達地表,通過在地表布設一種稱之為檢波器的接收裝置,接收來自不同深度地質界面反射的地震子波,其接收的信號的合集稱之為地震記錄。地震記錄可以認為是地震子波與一系列反射系數序列的褶積。當地層較薄時,薄層上下界面地震子波難以區分,來自不同界面的地震子波疊合在一起,使得從地震記錄上很難分辨這種很薄的地質界面。這種利用地震記錄分辨薄層的能力稱之為地震資料的分辨率。為了利用地震記錄分辨薄層,提高地震資料的分辨率,通常需要對地震記錄進行反褶積處理。
所述的反褶積處理,通常是指通過壓縮地震子波,將延續幾十甚至一百多毫秒的地震子波壓縮成原來的震源尖脈沖形式,使地震記錄變為反射系數序列的尖脈沖組合,提高分辨率的數據處理過程。現有技術中常用的反褶積方法通常是基于兩種基本假設(假設反射系數序列是白噪聲的,假設地震子波序列是最小相位的),其反褶積過程主要包括:假設地震記錄為
x(t)=S(t)+n(t)=Σt=0mb(t)ξ(t-τ)+n(t)---(1-1)]]>
其中S(t)為有效信號,n(t)為干擾波,b(τ)為地震子波,ξ(t)為反射系數序列,τ為反射系數序列的延時,m為反射系數個數。通常會假設不存在干擾波n(t),即:
x(t)=S(t)=b(t)*ξ(t)             (1-2)
對上式(1-2)兩邊求傅氏變換,則得到頻率域的地震記錄表示式:
X(ω)=B(ω)·ξ(ω)            (1-3)
上式(1-3)中,X(ω)、B(ω)和ξ(ω)分別為地震頻譜、子波頻譜和反射系數的頻譜。
上式(1-3)經過變化可以得到:
ξ(ω)1B(ω)·X(ω)---(1-4)]]>
如果令:
A(ω)=1B(ω)---(1-5)]]>
則有:
ξ(ω)=A(ω)·X(ω)            (1-6)
再對上式(1-6)式做反傅氏變換,將其變換至時間域,可得到:
ξ(t)=a(t)*x(t)=a(t)*b(t)*ξ(t)          (1-7)
上式(1-7)中,a(t)為A(ω)的時間函數。根據(1-7)式可知:
a(t)*b(t)=δ(t)           (1-8)
因為b(t)為地震子波,而a(t)和b(t)之間又存在著頻譜互為倒數的關系(即A(ω)=1/B(ω)),所以把通常a(t)稱為反子波,又叫做反褶積因子。
由此可知,若已知地震子波,利用數學方法求出a(t),再利用(1-7)式讓反子波a(t)與地震記錄x(t)做褶積,就可以求出反射系數序列ξ(t),即
ξ(t)=Στα(τ)x(t-τ)---(1-9)]]>
估計出地震子波的方法很多,常用的有同態濾波方法或統計性子波估算方法。估計出地震子波后,就可以求出反射系數的序列,達到地震子波壓縮成尖脈沖、提高地震記錄縱向分辨能力的目的。
上述反褶積方法基于兩種基本假設估計出地震子波,求出反射系數序列。但是所述的兩種基本假設常常與實際地層結構存在較大的差異,尤其是地震子波對弱信號的壓制,使得該方法在地震數據處理時降低了地震數據記錄的分辨率。
發明內容
本發明目的在于提供一種新的反褶積方法,該方法擺脫了兩種基本假設,減少了地震子波對地震數據中弱信號的干涉作用,能提高地震數據記錄的分辨率。
本發明的實現方法包括以下處理步驟:
S1:將所需處理的地震數據處理成疊后地震數據記錄;
S2:構建包括反射系數序列在內的目標函數;
S3:在所述疊后地震數據記錄中選取一道所需處理的地震資料,從所述構建的目標函數中計算出所述地震資料的反射系數序列的振幅和位置;
S4:讀取下一道所需處理的地震資料,計算出所述下一道地震資料反射系數序列的振幅和位置;
S5:計算出所述疊后地震數據記錄中所有地震道反射系數序列的振幅和位置,輸出反射系數序列,完成所述疊后地震數據記錄的反褶積處理。
上述所述的一種反褶積方法,其優選方案為,所述構建目標函數包括:
S201:建立地震數據褶積模型;
S202:將所述地震數據褶積模型轉換成矩陣形式;
S203:由所述矩陣形式構建目標函數;所述構建的目標函數為:
Jα=||Aα-S||2
上式中,Jα為關于反射系數序列α的范數,A為一個N×M的矩陣,其表達式為其中,ωt為地震子波,t=1,2,…N,N為所述疊后地震數據記錄的時間窗長度,M為預先設置的所述時間窗長度N內反射系數序列的個數,τM第M個反射系數的時延,S為t時刻所述疊后地震數據記錄的值。
上述所述的一種反褶積方法,其優選方案為,所述構建的目標函數為:
J=||Aα-S||2+μ||Cα-ε||2=Jα+μJc
上式中,μ為阻抗權系數,C為左下三角全為1的矩陣,ε為波阻抗。
上述所述的一種反褶積方法,其優選方案為,所述構建的目標函數為:
J=Jα+μJc+βJα
上式中,β為正則化的阻尼系數,其中,β=β0f0,β0為預設值,f0=max{(ATA+μCTCii)},i為矩陣(ATA+μCTC)的下標。
上述所述的一種反褶積方法,其優選方案為,從所述構建的目標函數中計算出所述地震資料中反射系數序列的振幅包括步驟:
S310:從所述疊后地震數據記錄中估計出地震子波;
S320:將所述估計出的地震子波代入所述目標函數中,計算出所述反射系數序列的振幅;
其中,所述估計地震子波的方法包括步驟:
S311:建立所述疊后地震數據記錄的褶積模型;
S312:將所述建立的褶積模型變換到復賽譜域;
S313:設計低通濾波器,在復賽譜域將所述建立的褶積模型的地震子波和反射系數序列 進行分離,得到復賽譜域的地震子波;
S314:將所述復賽譜域的地震子波轉換到時間域,得到時間域的地震子波。
上述所述的一種反褶積方法,其優選方案為,所述地震子波的估計方法還包括步驟:
315:對所述估計出的地震子波相位進行校正。
上述所述的一種反褶積方法,其優選方案為,采用模擬退火算法計算反射系數序列中每一個反射系數的時間位置,其操作過程包括以下步驟:
S301:在所述時間窗長度N內隨機給定反射系數序列a(j)的初始時間位置,由所述初始時間位置計算所述目標函數Ja當前解,j=1,2…M;
S302:在所述時間窗長度N內對所述反射系數序列a(j)進行擾動,將所述擾動后的反射系數序列a(j)代入所述目標函數Ja中,計算出所述目標函數Ja的新解;
S303:將所述Jα的新解與當前解進行比較,根據Metropolis準則判斷是否接受所述目標函數Jα的新解;
如果判斷結果為接受所述新解,則接受反射系數α(j)當前的時間位置;此時,所述新解作為目標函數Ja當前解;
S304:重復S302、S303直到達到收斂條件后,停止迭代,得到所述反射系數α(j)的時間位置。
上述所述的一種反褶積方法,其優選方案為,所述模擬退火算法還包括在所述S303之后:
S3031:設置誤差項σ,將所述新解與當前解的差值與所述誤差項σ比較;如果所述新解與當前解的差值大于所述誤差項σ,則不接受所述目標函數Jα的新解。
上述所述的一種反褶積方法,其優選方案為,在達到收斂條件后,將所述模擬退火算法得出的反射系數序列α(j)的振幅值作為地震資料中反射系數序列的振幅。
本發明還提供一種反褶積裝置,該裝置包括:數據預處理模塊、目標函數建立模塊、目標函數計算模塊、資料讀取模塊、反射系數序列輸出模塊,其中:
所述數據預處理模塊,用于將所需處理的地震數據處理成疊后地震數據記錄;
所述目標函數建立模塊,用于構建包括反射系數序列在內的目標函數;
所述目標函數計算模塊,用于在所述疊后地震數據記錄中選取一道所需處理的地震資料,從所述構建的目標函數中計算出所述地震資料的反射系數序列的振幅和位置;
所述資料讀取模塊,用于讀取下一道所需處理的地震資料,計算出所述下一道地震資料反射系數序列的振幅和位置;
所述反射系數序列輸出模塊,用于計算出所述疊后地震數據記錄中所有地震道反射系數 序列的振幅和位置,輸出反射系數序列,完成所述疊后地震數據記錄的反褶積處理。
上述所述的一種反褶積裝置,其優選方案為,所述目標函數建立模塊構建的目標函數為:
Jα=||Aα-S||2
上式中,Jα為關于反射系數序列α的范數形式的函數,A為一個N×M的矩陣,其表達式為其中,ωt為地震子波,t=1,2,…N,N為所述疊后地震數據記錄的時間窗長度,M為預先設置的所述時間窗長度N內反射系數序列的個數,τM第M個反射系數的時延,S為t時刻所述疊后地震數據記錄的值。
上述所述的一種反褶積裝置,其優選方案為,所述目標函數建立模塊構建的目標函數為:
J=||Aα-S||2+μ||Cα-ε||2=Jα+μJc
上式中,μ為阻抗權系數,C為左下三角全為1的矩陣,ε為波阻抗。
上述所述的一種反褶積裝置,其優選方案為,所述目標函數建立模塊構建的目標函數為:
J=Jα+μJc+βJα
上式中,β為正則化的阻尼系數,其中,β=β0f0,β0為預設值,f0=max{(ATA+μCTCii)},i為矩陣(ATA+μCTC)的下標。
上述所述的一種反褶積裝置,其優選方案為,所述目標函數計算模塊包括振幅計算模塊、位置計算模塊,其中:
所述振幅計算模塊,用于從所述目標函數建立模塊中計算出所述地震資料中反射系數序列的振幅;
所述位置計算模塊,用于從所述目標函數建立模塊中計算出所述地震資料中反射系數序列的時間位置。
上述所述的一種反褶積裝置,其優選方案為,位置計算模塊包括初始模塊、擾動模塊、迭代模塊、收斂模塊,其中:
所述初始模塊1,用于在所述時間窗長度N內隨機給定反射系數序列a(j)的初始時間位置,由所述初始時間位置計算所述目標函數Ja當前解,j=1,2…M;
在所述擾動模塊,用于在所述時間窗長度N內對所述反射系數序列a(j)進行擾動,將所述擾動后的反射系數序列a(j)代入所述目標函數Ja中,計算出所述目標函數Ja的新解;
所述迭代模塊,用于將所述Jα的新解與當前解進行比較,根據Metropolis準則判斷是否接受所述目標函數Jα的新解;
如果判斷結果為接受所述新解,則接受反射系數α(j)當前的時間位置;此時,所述新解作為目標函數Ja當前解;
所述收斂模塊,用于繼續計算所述新解,并與所述當前解進行比較直到達到收斂條件后,停止迭代,得到所述反射系數α(j)的時間位置。
本發明提供的反褶積方法是一種稀疏脈沖非線性反褶積方法,該方法將傳統的以子波壓縮為核心理念的反褶積方法轉移到求解反射系數位置和振幅上來,減小了地震子波對地震數據處理的影響,較大幅度地提高地震記錄的分辨率。同時通過對反射系數統計特征的有效約束,減小了反褶積結果的多解性,提高了該反褶積方法數據處理的結果穩定性。
附圖說明
圖1是本發明提供的一種反褶積方法實施例1的方法流程圖;
圖2是本發明實施例1中模擬退火算法的方法流程圖;
圖3是本發明實施例1中模擬退火算法的另一種實施方式的方法流程圖;
圖4是本發明提供的一種反褶積裝置模塊結構示意圖;
圖5是本發明反褶積裝置中目標函數計算模塊的模塊結構示意圖;
圖6是本發明所述目標函數中位置計算模塊的模塊結構示意圖;
圖7是采用基于兩種基本假設對數據進行反褶積處理后的效果圖;
圖8是采用本發明反褶積方法對同一地震數據反褶積處理后的效果圖。
具體實施方式
為了使本技術領域的人員更好地理解本申請中的技術方案,下面將結合本申請實施例中的附圖,對本申請實施例中的技術方案進行清楚、完整地描述,顯然,所描述的實施例僅僅是本申請一部分實施例,而不是全部的實施例。基于本申請中的實施例,本領域普通技術人員在沒有做出創造性勞動前提下所獲得的所有其它實施例,都應當屬于本發明保護的范圍。
圖1是本發明提供的一種反褶積方法實施例1的方法流程圖。如圖1所示,該方法包括以下處理步驟:
S1:采集地震數據,將所需處理的地震數據處理成疊后地震數據記錄。
地震數據的采集,其主要是指在油氣勘探的區域,布置二維或三維測線,使用炸藥震源或可控震源激發地震波,在測線上等間距布置多個檢波器來接收地震波信號,以等時間間隔離散采樣地震數據,并以數字形式收集記錄。在地震數據采集記錄中,每個檢波點上記錄的地震波數據稱為地震道。一個檢波點上的地震道為一個單道或一個道。可以為一個地震道設置一個地震道序號。采集后的地震數據通常會包括壞道、空道、隨機干擾及其他誤差,因此還需要進行預處理、校正等消除或減小地震數據記錄的誤差。
所述的預處理通常包括對采集的原始數據進行解碼、抽取地震道、剔除空道或壞道等一系列的前期處理過程。在校正過程中,可以包括靜校正和動校正處理。所述的靜校正可以是按靜校正量所確定的移動量,把地震道的振幅離散值進行整體移動,消除地面山谷、地表風化層和低速帶等地表異常的影響。所述的動校正可以是通過動校正公式計算出動校正量來消除地震波到達各不同檢波點的正常時差。可以在經動、靜校正處理后的地震數據記錄上地震道序號相同的地震道中,取其相對應采樣點的算術平均值作為所述相對應采樣點疊加后的新值,這樣由多個地震數據疊加組成的地震道的集合稱為疊后地震數據記錄。地震數據疊加的方法有很多,本發明方法地震數據疊加的目的是壓制多次波和隨機干擾,因此不限于本實施例中所述的地震數據疊加方法。
在地震數據處理中,對地震數據記錄的反褶積可以是對整個地震數據記錄反褶積處理,也可以僅對某一時間窗內的一段地震數據記錄反褶積處理。本申請中所述的所需處理的地震資料,可以為上述所述的兩種情況中任意一種。
采集地震數據后,將所需處理的地震數據處理成疊后地震數據記錄。
S2:構建包括反射系數序列在內的目標函數。
所述構建的目標函數可以包括反射系數序列,也可以同時包括地震子波和反射系數序列。所構建的目標函數可以為多次試驗或理論推導得出,并且能夠求解出至少一組可接受的反射系數序列的解。本實施例1中所述的構建目標函數可以包括以下步驟:
S201:建立地震數據褶積模型。
將地震數據記錄的褶積模型表示為st=wt*rt+nt,t=1,2,…N,N為所述疊后地震數據記錄的時間窗長度。上式中,st為地震數據記錄,ωt為地震子波,rt為稀疏脈沖(即反射系數序列),nt為隨機噪音。其中,所述稀疏脈沖rt的表達式可以為如下的形式:
rt=Σj=1Mαjδt-τj]]>
上式中,M是非零脈沖的個數,δt是單位脈沖函數,τj是脈沖的時延,αj是脈沖的振幅,j=1,2…M。
將所述系數脈沖rt表達式代入所述褶積模型中,可以得到:
st=Σj=1Mαjωt-τj+1+nt,t=1,2...N---(2-1)]]>
S202:將所述地震數據褶積模型轉換成矩陣形式。
上述公式(2-1)寫成矩陣形式進行展開,可以表示為:
w1-τ1+1w1-τ2+1···w1-τM+1w2-τ1+1w2-τ2+1···w2-τM+1············wN-τ1+1wN-τ2+1···wN-τm+1α1α2···αM+n1n2···nN=S1S2···Sn]]>
上述矩陣表達式中,至可以表達為一個N×M矩陣A,因此上述矩陣表達式也可以表示為:
Aα+n=s            (2-2)
上式(2-2)中,A是一個N×M的矩陣,其中
S203:由所述矩陣形式構建目標函數。
上述公式(2-2)中的A是一個N維的非線性方程組,表示為在地震數據記錄中t時刻地震子波的值。褶積模型合成的地震數據記錄是由地震子波和反射系數序列褶積而成,若已知實際地震數據記錄和地震子波,可以通過所述疊后地震記錄與合成的地震記錄的最小平方差來求解反射系數序列的振幅。基于上述思想,本實施例1中可以將所述目標函數定義如下:
Jα=||Aα-S||2        (2-3)
上式(2-3)定義的目標函數Jα是一個關于反射系數序列α的范數。A為一個N×M的矩陣,其表達式為其中,ωt為地震子波,t=1,2,…N,N為所述疊后地震數據記錄的時間窗長度,M為預先設置的所述時間窗長度N內反射系數序列的個數,τM第M個反射系數的時延,S為t時刻所述疊后地震數據記錄的值。
S3:在所述疊后地震數據記錄中選取一道所需處理的地震資料,從所述構建的目標函數中計算出所述地震資料的反射系數序列的振幅和位置。
從所構建的目標函數中求解出反射系數序列的振幅和位置,可以根據所構建的目標函數不同而有不同計算方法,具體的可以根據所構建的目標函數的特性進行計算。在本實施例1中,所構建的目標函數Jα=||Aα-S||2中,其計算反射系數序列振幅的過程可以包括以下步驟:
S310:從所述疊后地震數據記錄中估計出地震子波。
在本實施例1構建的目標函數中,涉及到未知參數量地震子波。因此,在求解目標函數時可以先由上述S1所述的疊后地震數據記錄估計出地震子波ωt。地震數據子波的估計方法有很多,包括同態濾波方法和統計性子波估計方法。不同的估計方法得出地震子波會有偏差。 在常規的反褶積處理中,所述地震子波ωt的偏差對最后反褶積處理的結果影響很大。而在本發明的反褶積方法中,即使所述估計出的地震子波ωt存在一定的誤差,也可以得出較為準確反射系數序列。
本實施例1中所述地震子波的估計方法可以包括以下步驟:
S311:建立所述疊后地震數據記錄的褶積模型。
所述疊后地震數據記錄的褶積模型可以表示為:
st=ωt*rt+nt
其中,st是已知的疊后地震數據記錄,ωt為地震子波,rt為反射系數序列,nt為噪聲項。
S312:將所述建立的褶積模型變換到復賽譜域。
不考慮噪聲項nt。對上述褶積模型等號兩邊同時做傅里葉變換,得到所述疊后地震數據記錄褶積模型在頻率域的表達式:
x(ω)=w(ω)ξ(ω)          (3-1)
對上述公式(2-1)x(ω)=w(ω)ξ(ω)等號兩邊取對數,轉化為線性系統的公式:
lnx(ω)=lnw(ω)+lnξ(ω)              (3-2)
對上述公式(2-2)lnx(ω)=lnw(ω)+lnξ(ω)作反傅里葉變換,得到:
x~(t)=w~(t)+ζ~(t)---(3-3)]]>
其中,分別稱為x(t),w(t),ξ(t)的復賽譜序列。
S313:設計低通濾波器,在復賽譜域將所述建立的褶積模型的地震子波和反射系數序列進行分離,得到復賽譜域的地震子波。
將地震記錄由時間域轉換為復賽譜域后,地震記錄可以表示為地震子波與反射系數序列之和。其中地震子波能量集中在原點附近。此時,可以通過設定一個“分割點”,對復賽譜域的地震記錄進行低通濾波,可以得到復賽譜域的地震子波。
S314:將所述復賽譜域的地震子波轉換到時間域,得到時間域的地震子波ωt。
在估計出地震子波后,為了得到最佳的地震子波相位,還可以包括以下處理步驟:
S315:對所述估計出的地震子波相位進行校正。
對所述地震子波相位進行校正的方法有很多,常用的相位校正方法有常相位校正方法。所述的常相位校正方法通常是指基于地震子波的剩余相位不依賴頻率變化的假設,對記錄的相位譜進行一個不依賴于頻率的常數的調整。本實施例1中所述的對地震子波相位進行校正不限于所述的常相位校正方法。
估計出一個地震子波后,將所述地震子波ωt作為已知量求解所述目標函數。
S320:將所述估計出的地震子波代入所述目標函數中,計算出所述反射系數序列的振幅。
將上述所述建立的目標函數(2-3)展開為:
Jα=(Aα-S)T(Aα-S)
=(αTAT-ST)(Aα-S)           (3-4)
=αTATAα-αTATS-STAα+STS
對上述公式(3-4)進行求導,可以得到:
∂Jα∂αT=AT-ATS=0---(3-5)]]>
由上述公式(3-5)可以得到求解所述反射系數序列α的表達式:
α=(ATA)-1ATS           (3-6)
上式(3-6)中,A為一個N×M的矩陣,AT為A的轉置,其表達式為其中,ωt為地震子波,t=1,2,…N,N為所述疊后地震數據記錄的時間窗長度,M為預先設置的所述時間窗長度N內反射系數序列的個數,S為t時刻所述疊后地震數據記錄的值。
所述M值是一個經驗值,該值可以由作業人員根據采集的地震數據記錄中單位時間內可能含有的反射系數序列個數自行設定。多次試驗證明,所述M值的設置可以略大于正常的經驗值,盡量減少遺漏的反射系數序列。例如,一個地震道的一段時間窗N內設置的M為20,實際該時間窗長度N內有15個反射系數序列,那么在本實施例1求解過程中可以求解出所述15個反射系數序列,其余的5個反射系數序列為0或近似為0。相反,如果所述時間窗N內實際有20個反射系數序列,而在實施本實施例1過程中設置為15,那么在求解過程中將會僅得出15個反射系數序列,其余5個反射系數序列將會被遺漏。
上述公式(3-6)計算得出的是一個地震道中所選取的地震資料中一組M個反射系數序列的振幅,因此還需要計算出所述M個反射系數序列中每個反射系數的時間位置。計算所述反射系數序列位置的方法可以根據所構建的目標函數選擇合適的方法。針對本實施例1所述構建的目標函數,本實施例1中可以采用模擬退火算法計算所述反射系數序列中每一個反射系數的時間位置。模擬退火(Simulated Annealing,簡稱SA)是一種通用概率算法,用來在一個大的搜尋空間(解空間)內采用迭代算法找尋命題的最優解。
圖2是本發明實施例1中模擬退火算法的方法流程圖。如圖2所示,所述采用模擬退火算法計算所述每一個反射系數的時間位置的過程可以包括以下步驟:
S301:在所述時間窗長度N內隨機給定反射系數序列a(j)的初始時間位置,由初始時間位置計算所述目標函數Ja當前解,j=1…M。
M個反射系數序列中的每個反射系數都在所述時間窗長度N內隨機給定初始時間位置。
S302:在所述時間窗長度N內對所述反射系數序列a(j)進行擾動,將所述擾動后的反射系數序列a(j)代入所述目標函數Ja中,計算出所述目標函數Ja的新解。
每次擾動后更新所述反射系數序列a(j)的時間位置。
所述的對反射系數序列a(j)進行擾動,在模擬退火算法迭代中通常是指對反射系數序列a(j)的時間位置給定一個擾動量,使其時間位置發生變化,進而使與所述反射系數序列a(j)有關的目標函數Jα的值發生變化。所述的擾動量可以是由其他函數隨機產生設定的,可以是一個增加的量,也可以是一個減小的量。當然,所述的擾動量也可以是自行設定的具有一定變化規律的擾動量。本實施例1中所述的擾動量可以使所述該模擬退火算法在所述時間窗長度N內對反射系數序列的時間位置搜索一遍即可,本發明方法對具體的擾動量的定義不做限定。
所述目標函數Jα的新解是相對于上一次擾動或初始值的解。所述目標函數Jα的上一次擾動的解或初始值的解可以統稱為當前解。
S303:將所述Jα的新解與當前解進行比較,根據Metropolis準則判斷是否接受所述目標函數Jα的新解。
若由Metropolis準則判斷接受所述目標函數Jα的新解,則將所述目標函數Jα的新解作為當前解繼續迭代計算,同時也接受所述反射系數α(j)當前的時間位置。此時,將所述新解作為目標函數Ja當前解。
S304:重復S302、S303直到達到收斂條件后,停止迭代,得到所述反射系數α(j)的時間位置。
所述的收斂條件,可以是所述模擬退火中定義的“退火計劃”,所述的“退火計劃”通常是指模擬退火中穩定的衰減函數,該函數是一個收斂函數。所述的收斂條件也可以是自定義的迭代次數或搜索解結束的停止條件。本實施例1中所述的收斂條件可以為所述模擬退火算法對當前所求的反射系數在所述時間窗長度N內的完成一遍時間位置搜索。
停止迭代后,可以得到的所述反射系數α(j)的最佳時間位置。
圖3是所述模擬退火算法的另一種實施方式。如圖3所示,上述所述的模擬退火中,為消除異常壞值,所述S303中可以包括以下步驟:
S3031:設置誤差項σ,將所述新解與當前解的差值與所述誤差項σ比較;如果所述新解與當前解的差值大于所述誤差項σ,則不接受所述目標函數Jα的新解。
所述誤差項σ為一個經驗值,其目的是防止求解結果與原所述疊后地震數據記錄誤差太大。該值可以根據計算需求自行設定。
需要說明的是,在上述所述采用模擬退火計算得出所述反射系數序列最佳時間位置的同 時,也計算得出了所述目標Jα的新解,即計算得出了所述反射系數序列α(j)的振幅。試驗證明,該模擬退火算法計算得出的所述反射系數序列α(j)最佳時間位置時其所述反射系數序列α(j)的振幅也為最優的解。因此,本發明方法實施例1中可以在達到收斂條件后,將所述模擬退火算法得出的所述反射系數序列α(j)的振幅值作為所述反射系數序列的振幅。
S4:讀取下一道所需處理的地震資料,計算出所述下一道地震資料反射系數序列的振幅和位置;
計算完所述選取的地震資料反射系數序列的振幅和位置后,繼續選取下一道地震資料計算其反射系數序列的振幅和位置。所述下一道地震資料的選取,可以按照當前所選取的地震道序列的順序依次選取,也可以按照其他可以選取完所有所需處理的地震道的方法。
S5:計算出所述疊后地震數據記錄中所有地震道反射系數序列的振幅和位置,輸出反射系數序列,完成所述疊后地震數據記錄的反褶積處理。
計算完所需處理的地震數據中所有地震道的反射系數序列的振幅和位置,即完成了所述疊后地震數據記錄的反褶積處理。
本實施例1所述的反褶積方法將傳統的以子波壓縮或兩種基本假設的方法轉移到求解反射系數序列的位置和振幅上,提供了一種新的反褶積思路。在本發明方法中可以采用模擬退火算法進行多次迭代計算出最佳的所述反射系數序列的振幅和位置,在此過程中可以將一些由人為制造的假象所破壞的地震數據橫向連續性降到最低,很大程度上提高反褶積處理結果的準確性。
本發明還提供所述反褶積方法的第二種實施方式。本發明實施例2提供的反褶積方法允許在最優化過程中可引入波阻抗約束。所述波阻抗具有阻力的含義,在反褶積過程中將波阻抗的因素考慮在內,可以提高地震數據記錄反褶積的準確性和橫向連續性。地震波在介質中傳播時,作用于某個面積上的壓力與單位時間內垂直通過此面積的質點流量之比,稱為波阻抗。
在地震數據記錄處理中,反射系數序列αj與波阻抗εi有如下的關系:
Σj=1MαjuTi-τj=ϵi---(4-1)]]>
將公式(4-1)經整理和變換后,寫成矩陣形式為:
10···011···0············11···1r1r2···rM=ϵ0ϵ1···ϵk-1]]>
上述矩陣形式可以表示為:Cα=ε,其中,C為左下三角全為1的矩陣。
取上述反射系數序列與波阻抗差的范數形式引入所述構建的目標函數中,則所述目標函數可以表示為:
J=||Aα-S||2+μ||Cα-ε||2=Jα+μJc          (4-2)
公式(4-2)中μ為阻抗權系數,是一個可以預先設置經驗值。將公式(4-2)展開為:
J=(Aα-S)T(Aα-S)+μ(Cα-ε)T(Cα-ε)
=(αTAT-ST)(Aα-S)+μ(αTCT-εT)(Cα-ε)
=αTATAα-αTATS-STAα+STS+μ(αTCTCα-αTCTε-εTCα+εTε)
對公式(4-2)求導可得:即ATAα-ATS+μCTCα-μCTε=0,
由上可以得出所述反射系數序列α的表達式可以為:
α=(ATA+μCTC)-1(ATS+μCTε)         (4-3)
上式(4-3)中,A為一個N×M的矩陣,其表達式為AT為A的轉置,ωt為地震子波,t=1,2,…N,N為所述疊后地震數據記錄的時間窗長度,M為預先設置的所述時間窗長度N內反射系數序列的個數,S為t時刻所述疊后地震數據記錄的值,μ為阻抗權系數,C為左下三角全為1的矩陣,ε為波阻抗。
上述實施例2提供的反褶積方法中引入波阻抗因素,可以提高該反褶積方法處理的準確性。
在上述所述求解目標函數Jα=||Aα-S||2的過程中,還可以引入正則化操作。所述的正則化是為了防止矩陣求解過程中出現“病態”,而加入的輔助性因子,可以增加方程求解的穩定性。在本實施例中1,在所述目標函數的正則化可以為加入因子βJα,則,所述構建的目標函數可以表達為:
J=Jα+μJc+βJα         (4-4)
上式(3-7)中,β為正則化的阻尼系數,其中β=β0f0。其中β0可以根據經驗值進行預 設定,通常設置在0.001左右。本實施例1中所述的β0可以設置為0.001。可以f0可以定義為:
f0=max{(ATA+μCTC)ii}           (4-5)
上式中,i為矩陣(ATA+μCTC)的下標。上式f0的含義可以表示為取矩陣(ATA+μCTC)對角線中的最大值。
所述目標函數加入正則化處理,很大程度上減少了解的非唯一性的可能性,保障了解的穩定性,提高了本申請發明方法對地震數據記錄進行反褶積的處理結果的準確性和穩定性。
本發明根據所述反褶積方法提供一種反褶積裝置。圖4是本發明提供的一種反褶積裝置模塊結構示意圖。如圖4所述,該裝置包括:數據預處理模塊1、目標函數建立模塊2、目標函數計算模塊3、資料讀取模塊4、反射系數序列輸出模塊5,其中:
所述數據預處理模塊1,可以用于將采集的所需處理的地震數據處理成疊后地震數據記錄。
所述目標函數建立模塊2,可以用于構建包括反射系數序列在內的目標函數。所述構建的目標函數可以為:
Jα=||Aα-S||2
上式中,Jα為關于反射系數序列α的范數形式的函數,A為一個N×M的矩陣,其表達式為其中,ωt為地震子波,t=1,2,…N,N為所述疊后地震數據記錄的時間窗長度,M為預先設置的所述時間窗長度N內反射系數序列的個數,τM第M個反射系數的時延,S為t時刻所述疊后地震數據記錄的值。
所述目標函數建立模塊2構建的目標函數也可以為:
J=||Aα-S||2+μ||Cα-ε||2=Jα+μJc
上式中,μ為阻抗權系數,C為左下三角全為1的矩陣,ε為波阻抗。
所述目標函數建立模塊2構建的目標函數還可以為:
J=Jα+μJc+βJα
上式中,β為正則化的阻尼系數,其中,β=β0f0,β0為預設值,f0=max{(ATA+μCTCii)},i為矩陣(ATA+μCTC)的下標。
所述目標函數計算模塊3,可以用于在所述疊后地震數據記錄中選取一道所需處理的地震資料,從所述構建的目標函數中計算出所述地震資料的反射系數序列的振幅和位置。
所述資料讀取模塊4,可以用于度取下一道所需處理的地震資料,計算出所述下一道地 震資料反射系數序列的振幅和位置。
所述反射系數序列輸出模塊5,可以用于計算出所述疊后地震數據記錄中所有地震道反射系數序列的振幅和位置,輸出反射系數序列,完成所述疊后地震數據記錄的反褶積處理。
圖5是所述目標函數計算模塊3的模塊結構示意圖。如圖5所示,所述目標函數計算模塊3,可以包括振幅計算模塊310、位置計算模塊300,其中:
所述振幅計算模塊310,可以用于從所述目標函數建立模塊2中計算出所述地震資料中反射系數序列的振幅。
所述位置計算模塊300,可以用于從所述目標函數建立模塊2中計算出所述地震資料中反射系數序列的時間位置。
圖6是所述位置計算模塊300的模塊結構示意圖。如圖6所示,所述位置計算模塊300可以包括初始模塊301、擾動模塊302、迭代模塊303、收斂模塊304,
其中,所述初始模塊301,可以用于在所述時間窗長度N內隨機給定反射系數序列a(j)的初始時間位置,由所述初始時間位置計算所述目標函數Ja當前解,j=1,2…M。
所述擾動模塊302,可以用于在所述時間窗長度N內對所述反射系數序列a(j)進行擾動,將所述擾動后的反射系數序列a(j)代入所述目標函數Ja中,計算出所述目標函數Ja的新解。
所述迭代模塊303,可以用于將所述Jα的新解與當前解進行比較,根據Metropolis準則判斷是否接受所述目標函數Jα的新解。如果判斷結果為接受所述新解,則接受反射系數α(j)當前的時間位置。此時,所述新解作為目標函數Ja當前解。
所述收斂模塊304,可以用于繼續計算所述新解,并與所述當前解進行比較直到達到收斂條件后,停止迭代,得到所述反射系數α(j)的時間位置。
本發明提供的一種反褶積裝置,采用了最小二乘法和正則化方法能準確預測地震數據記錄中反射系數序列的振幅,并且用模擬退火算法計算反射系數序列振幅的位置,提高地震數據處理后的分辨率。
試驗證明,本發明方法和裝置相比常規的基于兩種基本假設的方法,將傳統的以子波壓縮為核心理念的反褶積方法轉移到求解反射系數位置和振幅上來,可以較大幅度地提高地震記錄的分辨率。圖7采用基于兩種基本假設對數據進行反褶積處理后的效果圖,圖8是采用本發明反褶積方法對同一地震數據反褶積處理后的效果。由圖7和圖8對比可以看出,采用本發明反褶積方法處理之后的地震數據分辨率得到明顯改善,沉積內幕和接觸關系清晰自然。尤為難得的是,由于消除了反褶積之后剩余相位對弱反射信號的干涉影響,1700-1900ms之間的波阻抗差異較小的弱反射信號得到了有效恢復,為精細地質解釋和儲層描述提供更加 可靠的基礎數據。

關 鍵 詞:
一種 反褶積 方法 裝置
  專利查詢網所有資源均是用戶自行上傳分享,僅供網友學習交流,未經上傳用戶書面授權,請勿作他用。
關于本文
本文標題:一種反褶積方法及裝置.pdf
鏈接地址:http://www.wwszu.club/p-6140788.html
關于我們 - 網站聲明 - 網站地圖 - 資源地圖 - 友情鏈接 - 網站客服客服 - 聯系我們

[email protected] 2017-2018 zhuanlichaxun.net網站版權所有
經營許可證編號:粵ICP備17046363號-1 
 


收起
展開
鬼佬大哥大