基于全變差最小化約束的地震反射系數反演方法
【技術領域】
[0001] 本發明屬于地震油氣勘探技術領域，具體為一種基于全變差最小化約束共輒梯度 算法的地震反射系數反演方法。
【背景技術】
[0002] 作為地震油氣勘探中的重要步驟和環節，地震儲層預測一直在油氣勘探中扮演著 不可或缺的角色。隨著巖性油氣藏勘探開發的不斷深入，對小于地震波長四分之一的薄儲 層(薄層)預測問題日益突出。模型反演，稀疏脈沖反演等常規地震反演方法對薄儲層識別 能力有限，直接影響到后續的油氣藏評價和井位設計。
[0003] 地震反射系數反演是在頻率域實施的，可獲得高分辨率的時間域反射系數信息， 實現對小于調諧厚度的薄儲層進行有效識別，以達到提高地震數據分辨率和提高儲層預測 精度的目的。
[0004] 地震反射系數反演是用有限的頻域信息得到全域頻域信息，是一個欠定反演問 題，傳統共輒梯度方法求解欠定問題效率低且易受異常值影響，會導致反演結果背離實際 地層特性，與地震剖面匹配較差等不利后果，不能正確的反應出真實的地層信息。
[0005] 目前應用于地震反射系數估計的反演算法包括：Charles I .Puryear和John P.Castagna采用的最小二乘共輒梯度法。其基本思想是把共輒性與最小二乘方法相結合， 利用已知點處的梯度構造一組共輒方向，并沿這組方向進行搜素，求出目標函數的極小點。
[0006] Rui Zhang和John P.Castagna采用的基追蹤算法。它采用反射系數的范數作為稀 疏性的度量，通過最小化L1范數將反射系數稀疏表示問題定義為一類有約束的極值問題， 進而轉化為線性規劃問題進行求解。
[0007] 袁三一和王尚旭等采用的粒子群算法和列文伯格-馬夸爾特法的聯合算法。其首 先采用粒子群算法進行全局尋優，來反演出地層反射系數位置，然后采用列文伯格-馬夸爾 特法算法來精準地求取反射系數的數值大小。粒子群算法和列文伯格-馬夸爾特法算法相 聯合，反演結果更精確，收斂性更快。
[0008] 秦德文等采用的基于隨機搜索的蒙特卡羅算法。它是以概率和統計理論方法為基 礎的一種計算方法。將所求解的問題同一定的概率模型相聯系，用計算機實現統計模擬或 抽樣，以獲得問題的近似解。蒙特卡羅方法有很強的適應性，該方法的收斂性是指概率意義 下的收斂，因此問題維數的增加不會影響它的收斂速度。
[0009] Kelyn PaolaCastafio.和German Ojeda采用的模擬退火算法和遺傳算法。模擬退 火算法從某一較高初溫出發，伴隨溫度參數的不斷下降，結合概率突跳特性在解空間中隨 機尋找目標函數的全局最優解。而遺傳算法是模擬達爾文生物進化論的自然選擇和遺傳學 機理的生物進化過程的計算模型，是一種通過模擬自然進化過程搜索最優解的方法。應用 于反射系數估計問題時，遺傳算法效果比模擬退火算法稍好。
[0010] 柴新濤，李振春等采用的最小二乘QR分解算法。它是利用Lanezos方法求解最小二 乘問題的一種投影法。由于在求解過程中用到QR因子分解，在對數據誤差傳遞的壓制和求 解收斂效率上具有明顯的優越性。總體來講，這些算法都可以不同程度的解決反射系數的 反演問題。
【發明內容】
[0011] 針對現有技術中存在的不足，本發明的目的之一在于解決上述現有技術中存在的 一個或多個問題。例如，本發明的目的之一在于解決現有技術中存在的反射系數反演效率 低和易受異常值影響的問題，提供一種基于全變差最小化約束的共輒梯度地震反射系數反 演方法。該方法利用地層連續性約束，在頻率域提升共輒梯度反演的準確性和效率。通過全 變差最小化約束增加對反射系數的識別能力，抑制共輒梯度方法易受到數據變化擾動以及 反演結果地層跳變的問題，獲得高分辨率的時間域反射系數，提升地震數據的分辨率，提高 對小于調諧厚度的薄儲層反射系數識別能力。
[0012] 為了實現上述目的，本發明提供了一種基于全變差最小化約束的地震反射系數反 演方法。所述反演方法包括以下步驟：
[0013] A、從疊后地震數據S中提取地震子波w并進行傅立葉變換得到地震子波w的頻域表 示W(f)，其中：
[0014] ff(f)=FFT(w) (1)
[0015] B、從所述疊后地震數據S中取一個地震道的一個時窗內的地震數據s進行傅立葉 變換，得到地震數據s的頻域表示S(f)，其中：
[0016] S(f)=FFT(s) (2)
[0017] C、根據地震數據形成原理由所述地震數據s的頻域表示S(f)和所述地震子波w的 頻域表示w(f)得到反射系數的頻域表示R(f)，其中：
[0018] R(f)=S(f)/ff(f) (3)
[0019] D、以所述時窗中心為分析點，對所述時窗內的地震數據s的時域反射系數進行傅 立葉變換，得到反射系數的頻域表達R(fT，其中：
[0020]
⑷
[0021 ]在等式(4)中，N為由點數表示的分析數據長度，Ti代表反射系數對之間的時間間 隔，j為虛數單位。
[0022] E、在所述步驟C得到的反射系數的頻域表示R(f)和所述步驟D得到的反射系數的 頻域表達R(fV相等的基礎上，根據反射系數奇偶分解原理構建所述時窗內反射系數反演 的目標函數方程。
[0023] F、利用共輒梯度算法求解所述目標函數方程以得到所述時窗內反射系數的奇分 量和偶分量，并重構出所述時窗內的反射系數，其中，在每次迭代過程中對共輒梯度算法所 求得的解進行全變差最小化約束，同時將經全變差最小約束之后的解作為下一次迭代計算 的初始解。
[0024] G、在所述地震道上按步長step滑動所述時窗，得到下一個時窗的地震數據，重復 步驟B至F，直到所述時窗遍歷所述地震道的地震數據，完成所述地震道的反演，得到所述地 震道的反射系數。
[0025] H、取下一地震道的地震數據并按照步驟B到G進行處理，直到得到每一個地震道的 反射系數。
[0026] 根據本發明的基于全變差最小化約束的地震反射系數反演方法的一個實施例，所 述步驟E包括:用歐拉公式將所述步驟D得到反射系數的頻域表達R(fV中的指數項展開，得 到等式(5):
[0027]
(5)
[0028]根據奇偶分解原理，等式(5)中實部是反射系數對^與^1+1的偶分量，用re(i，N_i +1)來表示，虛部是反射系數對ri與rN-i+i的奇分量，用r〇(i，N-i+l)來表示，得到等式(6):
[0029]
(6)
[0030]將所述等式⑶分析點移至時窗中點，并結合等式(6)建立目標函數方程，得到等式(7):
[0031]
(7 )
[0032] 在等式(7)中，flOT為低頻截止頻率，fhigh為高頻截止頻率，~為偶分量權重，α。為奇 分量權重，Re代表著實部，Im代表著虛部，h為反射系數的偶分量，:r。為反射系數的奇分量， △ t為相對于分析點的位移量，tw為時窗win的半窗長。
[0033]將等式(7)在頻域和時域離散化并寫為奇分量和偶分量的形式，得到下面的等式 (8)和等式(9):
[0034]
[0035]
[0036] 在等式(8)和等式(9)中，是起始頻率，取值為所述flOT; fm是截止頻率，取值為所 high 〇
[0037] 將所述等式(8)和等式(9)分別寫為矩陣形式，得到下面的等式(10)和等式(11):
[0038] be = AeXre (10)
[0039] b〇 = A〇Xr〇 (11)
[0040] 建立目標函數方程：
[0041 ] (12)
[0042] 在等式（10)、（ 11)和（12)中，為偶分量權重，α。為奇分量權重;是頻域反射系數 的實部;b。是頻域反射系數的虛部;I為偶分量變換矩陣，BP:
[0043]
(13)
[0044] A。為奇分量變換矩陣，即：
[0045]
(.1.4).
[0046] 將所述目標函數方程(12)寫為更一般的優化問題，得到目標函數方程(15):
[0047] b=Ax (15)
[0048] 在等式(15)中，X為待求的反射系數。
[0049] 根據本發明的基于全變差最小化約束的地震反射系數反演方法的一個實施例，所 述步驟F包括：
[0050] 將所述目標函數方程（15)的求解問題轉化為| |Ax=b| |2優化問題，進一步寫成二 次函數形式：
[0051 ]
(16)
[0052] 在等式（16)中，Q=ATA，y = ATb，T表示轉置運算。
[0053] 初始化:叉。=0，1'。= 13-厶叉。，〇=||(>)||2，1^=1，1^為迭代次數。
[0054] 循環步驟：
[0055] (1)更新梯度方向：
[0056] gk = Qxk-i-y (17)
[0057] (2)更新共輒方向：
[0058]若k = l，則pk = r。；若k關 1，則
[0060] pk = _gk+akPk-1 (19)[0061] (3)更新第k次迭代的解：
[0059] (18) _2]
(2〇) [0063] (4)對xk進行全變差最小化，
[0064] ,~、 (21)
[0065]其中，D表示整體剖面集，采用對稱邊界條件，xu表示第i行第j個元素。
[0066]