全局優化的時變稀疏反褶積方法及裝置的制造方法
【技術領域】
[0001] 本發明涉及地球物理勘探反演技術領域，尤其涉及一種全局優化的時變稀疏反褶 積方法及裝置。
【背景技術】
[0002] 在地球物理勘探領域中，采集的地震波形資料可看成是地面激發的子波在經過地 下傳播和界面反射后被接收時的波形疊加，是地震子波與反射系數褶積的結果，即符合后 續地震資料處理和反演中常用的褶積模型。其中，反射系數與地下介質的彈性特征有關，蘊 含地下巖性分界面的信息，反射系數相比于地震波形資料具有更寬頻帶和更高分辨率。因 此，為了獲得高分辨率的反射系數，常用反褶積、反演等方法進行反射系數的求取，以期達 到壓縮子波、拓寬頻帶和提高分辨率的效果。
[0003] 然而，常用的反褶積方法假設地震資料是穩態的，地下介質是完全彈性的，即子波 在時空上是沒有變化的，但這與實際情況(子波在地下傳播的過程中，吸收衰減等作用會改 造子波的波形、頻譜和相位)不符，直接影響了常規反褶積算法的穩定性和準確性。其次，常 規反褶積方法是在時間域對疊后地震資料進行處理，相當于運用了全頻帶的地震資料信 息，然而地震資料在低于約5Hz和高于約100Hz的頻帶內是相對不準確的，因此容易引入高、 低頻噪音，使結果的精度受到限制。再次，反射系數的位置求取是一個非線性的過程，運用 常規線性近似手段求取，其結果容易陷入局部解，且依賴于初始模型，因此精度難免受到限 制。
【發明內容】
[0004] 本發明提供了一種全局優化的時變稀疏反褶積方法及裝置，以至少解決目前地質 勘探過程中，使用常規反褶積方法獲得反射系數，結果精度不高的問題。
[0005] 根據本發明的一個方面，提供了一種全局優化的時變稀疏反褶積方法，包括:構建 非穩態地震資料的頻率域時變的部分譜反褶積模型;利用有限更新率原則和所述非穩態地 震資料的有效帶寬確定反射系數脈沖個數;利用非常快速模擬退火算法，結合所述部分譜 反褶積模型，對反射系數脈沖的時間位置進行搜索，針對搜索到的各個時間位置，利用最小 二乘法進行迭代，確定各個時間位置處的反射系數脈沖幅值;對于每次迭代，根據系統能量 確定可用的反射系數脈沖時間位置和幅值;利用部分譜觀測數據和部分譜合成數據的平方 誤差，判斷迭代是否收斂;如果收斂，則停止迭代，根據所述反射系數脈沖個數、所述可用的 反射系數脈沖時間位置和幅值，重構非零反射系數序列。
[0006] 在一個實施例中，構建非穩態地震資料的頻率域時變的部分譜反褶積模型，包括：
[0007] 從所述非穩態地震資料提取時變地震子波，所述時變地震子波的表示式如下：
實中，(69(6))表不時變地震子波的相位譜，Α( ω )表不 時變地震子波的振幅譜，ΗΤ表示希爾伯特變換，ω表示角頻率；
[0008] 結合所述時變地震子波構建所述部分譜反褶積模型，所述部分譜反褶積模型如 下：
[0009]
[0010]其中，S(fm)表示非穩態地震資料的傅里葉變換在頻率fm處的值;m = 1，2，…，Μ，Μ 2 Κ;Μ表示在采樣頻域(flOT，fhlgh)上采樣的頻點個數;Κ表示非零的反射系數脈沖個數;a k表示 第k個反射系數脈沖的幅值;表示第k個反射系數脈沖的時間位置;i表示虛數單位;N(fm) 表示噪音的傅里葉變換在頻率fm處的值;W(f m)表示時變地震子波的傅里葉變換在頻率匕處 的值。
[0011] 在一個實施例中，利用有限更新率原則和所述非穩態地震資料的有效帶寬確定反 射系數脈沖個數，包括:有效帶寬為(fi?，fh lgh)的非穩態地震資料能夠重構的非零反射系 數脈沖個數K為K=[|fhlgh-f lOT| XL/2]int，其中，L表示非零反射系數脈沖的周期，[· ]int表 示向下取整。
[0012] 在一個實施例中，利用最小二乘法進行迭代，確定各個時間位置處的反射系數脈 沖幅值，包括：
[0013] 如果無噪聲且稀疏層數先驗約束K小于第一預設值，采用以下公式計算反射系數 脈沖的幅值：
[0014]
[0015] 其中，?表示獲得的反射系數脈沖幅值向量；
=e〇s(2%/m)表示復指數函數exp(-i2JiTkf m)的實部，sin(2^/m)表示復指數 函數eXp(-i2Mkfm)的虛部;τ??表示第k個反射系數脈沖的時間位置表示第m個采樣頻率； d表示已知的部分譜觀測數據;T表示矩陣轉置；
[0016] 如果噪聲較強且稀疏層數先驗約束K大于第二預設值，采用以下公式計算反射系 數脈沖的幅值：
[0017]
[0018] 其中，λ表示阻尼系數;I表示單位矩陣。
[0019] 在一個實施例中，對于每次迭代，根據系統能量確定可用的反射系數脈沖時間位 置和幅值，包括：
[0020] 采用以下公式計算所述系統能量J:
[0021] J= I |Ga-d| |2+λ| |α| I2,
[0022] 其中
，6二=^〇8(2;^乂,)表示復指數函數 exp(-i23iTkfm)的實 部，C4 =-sin(2A./:J表示復指數函數exp(-i2Hkfm)的虛部; Tk表示第k個反射系數脈沖的 時間位置;fm表示第m個采樣頻率;a表示反射系數脈沖的幅值向量，a = [ai，a2，…，aK]T，ak表 示第k個反射系數脈沖的幅值;Ga表示部分譜合成數據;d表示已知的部分譜觀測數據；I卜 | 2表示求取數據的平方和，再開根號;λ表示阻尼系數；
[0023] 如果所述系統能量J達到預設范圍，則確定對應的反射系數脈沖時間位置和幅值 是可用的。
[0024] 在一個實施例中，利用部分譜觀測數據和部分譜合成數據的平方誤差，判斷迭代 是否收斂，包括：
[0025] 采用以下公式計算所述平方誤差Ja:
[0026] Ja= I |Ga-d| I2,
[0027] 其中，
良示復指數函數eXp(-i23iT kfm)的實 部，表示復指數函數eXp(-i2Mkfm)的虛部表示第k個反射系數脈沖的 時間位置;fm表示第m個采樣頻率;a表示反射系數脈沖的幅值向量，a = [ai，a2，…，aK]T，ak表 示第k個反射系數脈沖的幅值;Ga表示部分譜合成數據;d表示已知的部分譜觀測數據；I卜 | 2表示求取數據的平方和，再開根號；
[0028] 如果所述平方誤差Ja達到預設閾值范圍，則確定迭代收斂。
[0029]在一個實施例中，根據所述反射系數脈沖個數、所述可用的反射系數脈沖時間位 置和幅值，重構非零反射系數序列，包括：
[0030] 將所述反射系數脈沖個數、所述可用的反射系數脈沖時間位置和幅值帶入到非零 反射系數序列表達式中，所述非零反射系數序列r(t)的表達式如下：
[0031]
[0032]其中，K表示非零的反射系數脈沖個數;Tk表示第k個反射系數脈沖的時間位置;ak 表示第k個反射系數脈沖的幅值;δ (t)表示脈沖函數。
[0033] 在一個實施例中，所述方法還包括:在迭代過程中，如果檢測到迭代次數達到預設 的最大迭代次數，則停止迭代。
[0034] 根據本發明的另一個方面，提供了一種全局優化的時變稀疏反褶積裝置，包括:構 建單元，用于構建非穩態地震資料的頻率域時變的部分譜反褶積模型;第一確定單元，用于 利用有限更新率原則和所述非穩態地震資料的有效帶寬確定反射系數脈沖個數;迭代單 元，用于利用非常快速模擬退火算法，結合所述部分譜反褶積模型，對反射系數脈沖的時間 位置進行搜索，針對搜索到的各個時間位置，利用最小二乘法進行迭代，確定各個時間位置 處的反射系數脈沖幅值;第二確定單元，用于對于每次迭代，根據系統能量確定可用的反射 系數脈沖時間位置和幅值;判斷單元，用于利用部分譜觀測數據和部分譜合成數據的平方 誤差，判斷迭代是否收斂;重構單元，用于在迭代收斂的情況下，停止迭代，根據所述反射系 數脈沖個數、所述可用的反射系數脈沖時間位置和幅值，重構非零反射系數序列。
[0035]在一個實施例中，所述構建單元包括：
[0036] 提取模塊，用于從所述非穩態地震資料提取時變地震子波，所述時變地震子波的 表示式如下
其中，表示時變地震子波的相位 譜，Α( ω )表示時變地震子波的振幅譜，HT表示希爾伯特變換，ω表示角頻率；
[0037] 構建模塊，用于結合所述時變地震子波構建所述部分譜反褶積模型，所述部分譜 反褶積模型如下：
[0038]
[0039] 其中，S (fm)表示非穩態地震資料的傅里葉變換在頻率fm處的值;m = 1，2，…，Μ，Μ 2 Κ;Μ表示在采樣頻域(flOT，fhlgh)上采樣的頻點個數;Κ表示非零的反射系數脈沖個數;a k表示 第k個反射系數脈沖的幅值;表示第k個反射系數脈沖的時間位置;i表示虛數單位;N(fm)