本發明涉及的是空間非均勻采樣下不規則缺失道的地震數據重建方法，具體是一種基于線性bregman算法的非均勻曲波三維地震數據重建方法。

技術背景

在野外數據采集過程中，為了得到相對完整而規則的地震數據，在數據采集之前必須進行野外觀測系統的設計，但由于采集設備、野外地形條件以及經濟成本限制等原因，地震炮點和檢波點通常會偏離原始設計位置，甚至有些炮點檢波點無法采集到有效的地震數據，從而導致地震數據沿空間方向常進行不規則欠采樣，出現空間假頻，影響到后續其它處理方法的效果，降低了地震勘探的分辨率(trad,2009；張華等，2013)。為了解決這種問題，最直接也是最有效的方法就是在野外重新進行數據采集，但從經濟角度出發，顯然不可能重新進行數據采集來完美的解決該問題。因此，必須在室內采取相應的不規則地震數據重建方法，使得缺失道數據得到有效地恢復。然而在目前已有的重建方法中，包括預測濾波方法(spitz,1991；naghizadehandsacchi,2007)，降秩方法(oropezaandsacchi,2011；ma,2013)，數學變換方法(白蘭淑等，2014；唐歡歡等,2014)等，很少能夠重建空間非均勻采樣下的地震缺失道。但由于野外復雜地形條件的限制或者海上電纜的水平偏移，很多情況下野外地震數據常進行空間非均勻采樣，如不加處理則會引起覆蓋次數的變化(地下不均勻照明)，在疊加成像時會形成扭曲的成像振幅(采集腳印)，加重了空間假頻現象的出現，影響后續成像處理。

為了解決空間非均勻采樣下地震道不規則重建問題，地球物理領域常規處理方法為共面元疊加，從而將非均勻采樣數據歸位到均勻采樣數據中來，滿足后續其他處理方法的要求。然而共面元疊加處理方法忽略了每個面元內各道共中心點的真實位置，改變了部分地震道的振幅和相位，從而導致部分地震道位置出現嚴重偏差，降低了地震勘探資料的分辨率。另外一種方法就是基于波動方程的重建方法，然而該方法需要地下結構的先驗信息，計算量非常巨大，對采樣率要求也較高，從而也不能較好的解決該問題(ronen,1987)。盡管如此，許多學者采用基于數學變換的重建方法對該問題進行處理，duijndam等(1999)提出基于傅立葉變換的二維非均勻數據重建方法，hindriks和duijndam(2000)將其擴展到三維，實現了三維非均勻采樣重建技術。但是，duijndam等人的傅立葉重建方法有其局限性，重建結果受最低速度和空間帶寬的影響很大，隨著采樣間隔的逐漸增大，重建結果會逐漸變差。xu等(2005)采用重新正交化的過程，提出了基于抗泄露傅里葉變換的重建方法，但該方法對具有嚴重假頻的地震數據則重建效果不好。zwartjes等(2007)也提出了反假頻的非均勻地震數據重建方法，該方法采用無假頻的低頻信息來重建有假頻的高頻信息，達到壓制假頻和恢復缺失道的目的，但當低頻也具有假頻時，則該方法則會失效。jin(2010)提出基于阻尼最小范數傅立葉反演下的五維地震數據重建，該方法引入非均勻傅里葉變換方法，能夠重建空間非均勻采樣下的不規則缺失地震數據，但該方法抗假頻能力不強。況且以上方法都是采用傅里葉變換作為稀疏基，盡管計算速度較快，但只適合處理近似線性同相軸或者平穩變化的地震信號，不能解決非線性同相軸或者非平穩地震數據的重建問題。

曲波變換能夠表征信號的局部細節特征，可以有效地重建非線性同相軸或者非平穩變化的地震數據，眾多的研究結果也證明，基于曲波變換的數據重建方法效果顯著(naghizadehandsacchi,2010；劉國昌等，2011；張華等,2015)。盡管如此，以往基于曲波變換的重建方法前提條件仍然是空間均勻采樣，而對于空間非均勻采樣信號，由于缺乏空間連續性，曲波變換則不能有效地探測出地震波前特征，從而導致以往曲波變換不能有效地重建空間非均勻采樣下的地震缺失道，從而限制了該方法的進一步應用。

技術實現要素：

本發明的目的是為了能夠利用二維非均勻地震數據空間信息，高精度重建野外非均勻采樣下不規則缺失道，并大幅度提高重建信號的保真度和信噪比，保護微弱的有效波信號，從而使反射波同相軸更加連續，而提出了一種基于線性bregman算法的非均勻曲波三維地震數據重建方法。

本發明提出了一種基于線性bregman算法的非均勻曲波三維地震數據重建方法，首先針對常規二維曲波變換方法難以對空間非均勻采樣下地震缺失道進行重建以及基于二維非均勻曲波變換的重建方法精度不高，應用不廣泛等問題。沿著兩個空間方向依次抽取非均勻三維地震數據的時間切片，在多尺度多方向二維曲波正變換的基礎上，引入二維空間非均勻快速傅里葉變換，建立均勻曲波系數與空間非均勻采樣下地震缺失數據之間的非均勻曲波反變換算子，然后采用線性bregman算法進行求解，通過選擇合適的閾值因子和動態步長，并采用軟閾值算子，精確地反演計算得到非均勻采樣下不規則地震數據的均勻曲波系數，最后再進行常規曲波反變換，從而形成了一種基于線性bregman算法的非均勻曲波三維地震數據重建方法。

進一步，地震數據的重建問題就是從不完整的數據中恢復出完整的地震數據，假設如下線性正演模型

y＝md

這里y∈rm代表采集的不完整地震數據；d∈rⁿ，且n＞＞m，表示待重建的完整數據；m∈r^n×m表示隨機采樣矩陣。假設數據x是d在曲波變換域c中的稀疏表示，則上述方程可以寫成：

y＝ax且

這里上標h代表共軛轉置矩陣。

從上述方程可知，以往重建方法采用的算子c^h為常規均勻曲波反變換算子。事實上，該算子也可以為非均勻曲波反變換算子。為此，本發明采用和來分別表示常規均勻曲波反變換算子和非均勻曲波反變換算子。

進一步，所述二維曲波正變換的定義為：

式中：ψj,l,k表示曲波函數，c(j,l,k)為曲波系數，j,l,k分別表示尺度，方

向和位置參數，d(x)表示地震數據，其反變換為：

進一步，常規的二維曲波正變換有四個步驟：(1)對地震數據應用二維傅立葉變換，得到頻率波數域系數；(2)在頻率波數域形成角度楔形；(3)將每一個楔形圍繞到原點進行重新裝配；(4)對每一個裝配好的楔形應用二維傅立葉反變換，得到離散曲波系數。為此，定義常規曲波正變換算子

在這個方程中，f表示二維傅里葉變換，它實現了離散曲波變換第1步。t表示曲波拼接算子，即將頻率波數域變換到曲波系數的過程，它實現了離散曲波變換第2～4步。

由于常規曲波正反變換算子滿足因此可以定義均勻曲波反變換算子：

其中f^h表示二維傅里葉反變換，t^h表示曲波平鋪算子。

在上述方程中，由于二維傅立葉變換參與了快速離散曲波變換之中，因此常規二維曲波變換不能處理非均勻采樣數據，然而可以沿著時間切片的兩個空間方向采用二維非均勻快速傅里葉反變換代替常規二維快速傅立葉反變換算子f^h。因此新的非均勻曲波反變換算子可以定義為：

該算子可以將離散曲波系數與非均勻采樣下不規則地震道建立相應的聯系。此時由于重建方程的稀疏解可以求解以下l1范數最優化問題得到：

subjecttoy＝ax

從而可以得出均勻曲波系數。在這個表達式中，x代表估計值，l1范數定義為x[i]是向量x中第i個元素。而對于該方程的求解，本發明選用線性bregman算法來求解此最小化問題；

在求解上式方程得到曲波系數后，重建后的地震波場d可以通過下式得到：

其中，是標準的二維均勻曲波反變換算子，通過該式重建出空間均勻采樣下規則而完整的地震數據，實現一種基于線性bregman算法的非均勻曲波三維地震數據重建方法。

進一步，所述非均勻快速傅里葉變換主要實現過程如下：首先將非均勻地震數據與某高斯短濾波器進行褶積，并對其結果進行均勻網格下密集采樣，然后對密集采樣后的數據進行快速傅里葉變換到頻譜域，最后在頻譜域進行反褶積校正，得到非均勻地震數據的頻譜。

進一步，所述的線性bregman算法如下：

首先可以將l1范數下的最優化問題轉為求解下述bp規則化問題

其中，λ是一個閾值權衡因子，在平衡l1范數和l2范數起到重要作用，算法1給出了線性bregman算法求解非均勻曲波變換過程的偽代碼。

算法1.求解非均勻曲波變換的線性bregman算法

該函數用來處理數據中的噪聲，動態步長tk被定義為：

算法1中的閾值為

算法1中的軟閾值函數為

sλ＝sign(x)·max(|x|-λ,0)

本發明首先沿著兩個空間方向依次抽取非均勻三維地震數據的時間切片，在多尺度多方向二維曲波正變換的基礎上，引入二維空間非均勻傅里葉變換，建立均勻曲波系數與空間非均勻采樣下地震缺失數據之間的非均勻曲波反變換算子，然后采用線性bregman算法進行求解，通過選擇合適的閾值因子和動態步長，并采用軟閾值算子，精確地反演計算得到非均勻采樣下不規則地震數據的均勻曲波系數，最后再進行常規曲波反變換，從而形成了一種基于線性bregman算法的非均勻曲波三維地震數據重建方法。

本發明創造性主要體現在：

1、相對于共面元疊加和波動方程方法，本發明技術可以在不改變振幅和相位的前提下直接重建出缺失道地震數據，具有較好的保真度，并且不需要地下結構先驗信息，計算工作量較少。

2、相對于傅里葉變換，本發明技術更能夠表征具有曲線狀特征的地震信號局部細節特征，重建精度高，而且具有反假頻能力。

3、相對常規二維曲波變換方法，本發明技術不僅能夠處理空間均勻采樣下的地震道缺失重建，而且能夠處理空間非均勻采樣下的地震道缺失重建。大幅度提高了重建信號的信噪比，保護了微弱的有效波信號，從而使反射波同相軸更加連續、清晰，

4、相比于常規二維非均勻重建方法，本發明技術充分利用三維地震數據信息，從不同空間方向對缺失道進行重建，進一步提高重建精度，使得該發明在實際應用范圍更為廣泛、適用。

5、本發明技術對于其它非均勻信號處理領域也具有重要的借鑒意義。

綜述，本發明克服了常規二維曲波變換方法不能重建非均勻采樣下地震缺失道的缺點以及常規二維非均勻重建方法重建精度低，應用不廣泛的問題。并且該發明方法不僅可以重建非均勻帶假頻的缺失數據，而且也可以將非均勻網格數據歸位到任意指定的均勻采樣網格，大幅度提高了重建信號的分辨率和信噪比，保護了微弱的有效波信號，從而使反射波同相軸更加連續。這對于指導復雜地區非均勻地震數據采集、缺失道重建等方面具有重要的實用價值，同時對于其它非均勻信號處理領域也具有重要的借鑒意義。

附圖說明

圖1是本發明實施例中非均勻數據重建流程圖。

圖2是原始理論三維地震數據圖。

圖3是非均勻三維地震數據圖。

圖4是非均勻曲波變換三維地震數據規則化重建結果圖。

圖5是非均勻三維地震數據重建結果誤差圖。

圖6是50％非均勻采樣下的三維不規則缺失地震數據。

圖7是50％非均勻采樣下三維不規則缺失地震數據重建結果圖。

圖8是非均勻采樣下規則缺失三維地震數據圖。

圖9是非均勻采樣下規則缺失三維地震數據重建結果圖。

具體實施方式

以下實施案例用于說明本發明，但不用來限制本發明的范圍。

實施例1

實現該方法的步驟主要包括，地震數據重建模型建立，常規二維曲波正反變換，非均勻快速傅里葉變換，建立非均勻曲波正反變換算子，地震波場重建，線性bregman算法求解等。具體步驟如下：

步驟1：數據重建模型建立。地震數據的重建問題就是從不完整的數據中恢復出完整的地震數據，假設如下線性正演模型

y＝md

這里y∈r^m代表采集的不完整地震數據；d∈rⁿ，且n＞＞m，表示待重建的完整數據；m∈r^n×m表示隨機采樣矩陣。假設數據x是d在曲波變換域c中的稀疏表示，則上述方程可以寫成：

y＝ax且

這里上標h代表共軛轉置矩陣。

從上述方程可知，以往重建方法采用的算子c^h為常規均勻曲波反變換算子，事實上，該算子也可以為非均勻曲波反變換算子。為此，本發明采用和來分別表示常規均勻曲波反變換算子和非均勻曲波反變換算子。

步驟2：二維曲波正反變換。為了得到高精度地震數據重建結果，需要在二維曲波變換的基礎上進行數據重建，二維曲波正變換的定義為：

式中：ψj,l,k表示曲波函數，c(j,l,k)為曲波系數，j,l,k分別表示尺度，方向和位置參數，d(x)表示地震數據，其反變換為：

步驟3：非均勻傅里葉變換。首先將非均勻地震數據與某高斯短濾波器進行褶積，并對其結果進行均勻網格下密集采樣，然后對密集采樣后的數據進行快速傅里葉變換到頻譜域，最后在頻譜域進行反褶積校正，得到非均勻地震數據的頻譜。

步驟4：建立非均勻曲波正反變換算子。常規的二維曲波變換主要進行了四個步驟，(1)對地震數據應用二維傅立葉變換，得到頻率波數域系數；(2)在頻率波數域形成角度楔形；(3)將每一個楔形圍繞到原點進行重新裝配；(4)對每一個裝配好的楔形應用二維傅立葉反變換，得到離散曲波系數。為此，定義常規曲波正變換算子

在這個方程中，f表示二維傅里葉變換，它實現了離散曲波變換第1步。t表示曲波拼接算子，即將頻率波數域變換到曲波系數的過程，它實現了離散曲波變換第2～4步。

由于常規曲波正反變換算子滿足因此可以定義均勻曲波反變換算子：

其中f^h表示二維傅里葉反變換，t^h表示曲波平鋪算子。

在上述方程中，由于二維傅立葉變換參與了快速離散曲波變換之中，因此常規二維曲波變換不能處理非均勻采樣數據，然而可以沿著時間切片的兩個空間方向采用二維非均勻快速傅里葉反變換代替常規二維快速傅立葉反變換算子f^h。因此新的非均勻曲波反變換算子可以定義為：

該算子可以將離散曲波系數與非均勻采樣下不規則地震道建立相應的聯系。此時由于重建方程的稀疏解可以求解以下l1范數最優化問題得到：

subjecttoy＝ax

從而可以得出均勻曲波系數。在這個表達式中，x代表估計值，l1范數定義為x[i]是向量x中第i個元素。而對于該方程的求解，本發明選用線性bregman算法來求解此最小化問題；

步驟5：地震波場重建。在求解上式方程得到曲波系數后，重建后的地震波場d可以通過下式得到：

其中，是標準的二維均勻曲波反變換算子，通過該式重建出空間均勻采樣下規則而完整的地震數據，實現一種基于線性bregman算法的非均勻曲波三維地震數據重建方法。

步驟6：線性bregman算法。首先可以將l1范數下的最優化問題轉為求解下述bp規則化問題

其中，λ是一個閾值權衡因子，在平衡l1范數和l2范數起到重要作用，算法1給出了線性bregman算法求解非均勻曲波變換過程的偽代碼。

算法1.求解非均勻曲波變換的線性bregman算法

該函數用來處理數據中的噪聲，動態步長tk被定義為：

算法1中的閾值為

算法1中的軟閾值函數為

sλ＝sign(x)·max(|x|-λ,0)

實現該方法具體操作為：

為了衡量數據重建方法的效果，定義信噪比公式snr＝20log10||d0||2/||d-d0||2來進行對比，單位為db，其中d0表示原始模型數據，d表示重建結果，信噪比越高，代表重建結果與模型數據越接近，處理效果越理想。

首先建立一個二維不均勻速度模型，設置256炮，每炮256道接收，炮距和道距都為12米，采樣率4毫秒，采用聲波有限差分算法得到該模型所有的二維地震記錄，然后對所得到的正演地震數據沿檢波器，炮點以及時間坐標排列成三維數據體。圖2為原始理論三維地震數據圖，由于不能全部顯示理論的三維數據模型，只能從三個不同的方向進行顯示，其中時間切片為0.44s，炮點和檢波點距離都為1524m。首先為了驗證本發明所提方法的重建效果，對理論三維地震數據的時間切片進行二維均勻傅里葉變換，然后再進行空間二維非均勻快速傅里葉反變換，得到新的空間非均勻采樣下的256×256道三維地震數據，如圖3所示(圖3表示非均勻三維地震數據圖)，此時信噪比為6.47db，盡管名義上的采樣網格(道距和炮距)還是12×12米，但每道地震數據道距不均勻，其道距范圍為0m～24m，顯然，如果將非均勻采樣下地震數據在均勻采樣網格上進行顯示，連續的地震波場則會被破壞，而常規重建方法的前提條件是均勻道距下不規則缺失的地震記錄。為此采用本發明方法進行規則化重建，規則化重建后的采樣網格為12m×12m，結果如圖4所示(圖4表示非均勻曲波變換三維地震數據規則化重建結果圖)，規則化重建后的信噪比為42.23db，地震波場連續性顯著提高，重建后的地震記錄與原始記錄非常接近。圖5表示非均勻三維地震數據重建結果誤差圖，可以看出誤差幾乎忽略不計，重建后信噪比非常高，與原始地震記錄相比較幾乎沒有視覺上的差異，從而說明本發明方法能夠反映出地震波場的局部細節特征，因此重建方法精度高，保真度較好。

為了檢驗本發明方法在非均勻采樣下不規則缺失道重建效果，同樣采用空間非均勻傅里葉變換，并且對其進行50％非均勻采樣，得到50％非均勻采樣下的三維不規則地震數據，如圖6所示(圖6表示50％非均勻采樣下的三維不規則缺失地震數據)，此時三維地震記錄炮距或道距范圍為0m～84m，可以看出地震數據道距極不均勻，同相軸不連續，不能直接應用于后續資料的處理，必須采用非均勻三維地震數據重建方法進行處理。為此利用本發明方法進行重建，重建后的采樣網格為12m×12m，重建結果如圖7所示(圖7表示50％非均勻采樣下三維不規則缺失地震數據重建結果圖),重建后的信噪比為21.04db，可以看出盡管不均勻采樣下地震數據缺失50％地震道，但是重建效果精度仍然較高，重建后同相軸更連續，能量損失較少。

為了進一步檢驗本發明方法的反假頻能力，對原始理論非均勻三維地震數據進行規則欠采樣，如圖7所示(圖7表示非均勻采樣下規則缺失三維地震數據圖)，規則欠采樣帶來較為嚴重的假頻成分，與信號的真實頻譜存在部分重疊，然后采用本發明方法進行重建，重建后道距為12m×12m，重建結果如圖9所示(圖9表示非均勻采樣下規則缺失三維地震數據重建結果圖)，重建后信噪比為18.55db，重建效果較好，缺失的地震道得到了有效的恢復，表明本發明方法具有較強的反假頻能力，能夠進行復雜地區非均勻采樣下的不規則和規則缺失地震數據重建。

對于本領域技術人員而言，顯然本發明不限于上述示范性實施例的細節而且在不背離本發明的精神或基本特征的情況下，能夠以其他的具體形式來實現本發明的實用功能。