基于加速Bregman方法和閾值迭代法的聯合地震數據重建

龐 洋 張 華 郝亞炬 彭 清 梁 爽 韓紫璇

(東華理工大學核資源與環境國家重點實驗室，江西南昌 330013)

0 引言

在陸上地震勘探中，由于地震數據采集現場復雜的地震地質條件，加上采集設備、經濟成本及禁采區等因素的限制，往往會產生壞道、地震道缺失、道間距過大等情況，從而造成所采集地震數據的不完整[1]，包括有效信號缺失、信噪比低等。在海上，由于拖纜、洋流等因素也會造成空間假頻、繞射噪聲等問題，導致難以獲得較高質量地震數據。為了滿足后續處理對地震數據質量的較高要求，最直接有效的方法就是對缺失地震道進行補充采集和加密[2-3]。但受經濟成本制約，大多難以再回現場補采數據，因此只能在后期采用有效的地震數據重建方法重建野外缺失的地震道。

地震數據重建方法有以下四種主要類型。一是基于濾波器的方法[4]，該類方法通過插值濾波實現。因通常將非均勻網格采樣數據當作規則數據處理，故易造成較大誤差。二是波場延拓方法[5]，該方法充分利用地下信息。因地下結構等先驗信息通常是未知的，從而限制了該方法的應用。三是快速降秩方法[6-7]，該方法將插值問題看作圖像填充問題，計算速度快，參數設置簡單，但在非均勻網格采樣下的不規則缺失道重建及其反假頻能力方面還存在局限性。四是基于數學變換方法，這類方法不需預知地下結構等先驗信息，即能重建規則和不規則缺失地震道，且計算速度快，精度高。如Radon變換[8-9]、傅里葉變換[10-11]、小波變換[12]、Curvelet變換[13-15]、Seislet變換[16-17]等，是此類中的常用方法。本文研究的方法即屬于第四類，通過改進傳統算法，有效地重建出缺失地震數據。目前，基于傅里葉變換的地震數據重建方法已經產業化[18-19]，但其重建精度有待進一步提高。為此，眾多學者在傅里葉變換的基礎上對Curvelet進行研究。Hennenfent等[20]提出了基于Curvelet數據重建的稀疏促進反演方法，成功重建了含有缺失道的地震數據。張華等[15]提出基于曲波變換和新閾值參數的重建方法，且實現了頻率域中的三維地震數據重建。同時，人們還研究了基于曲波變換的重建方法[21-23]。但上述方法都是采用單一的迭代方法對地震數據進行重建，而單一方法難免存在局限性，如重建精度不夠、收斂速度較慢，或需迭代較多次數才能達到預期的重建精度，這些均限制了相應方法的工業應用。

于是，研究人員提出了許多加速迭代方法，其中一類為Bregman正則化方法。該方法最早由Osher等[24]、Yin等[25]引入圖像處理中，用于解決全變分(TV)的圖像恢復問題。隨后劉洋等[26]、郭萌等[27]、Cai等[28-29]應用線性Bregman方法進行地震數據重建，均取得較好效果。該方法參數設置簡單，通過只閾值化下一次的梯度項系數加速收斂，提高迭代速度； 但在每一次迭代過程中，由于部分未閾值化的系數參與重建，因此相比其他迭代法，線性Bregman方法重建效果有所下降。

由于迭代閾值收縮算法(Iterative Shrinkage Thresholding Algorithm，ISTA)參數設置簡單[30-32]，重建效果顯著，在反演領域得到了廣泛應用。但其運算速度較慢，需要迭代很多次才能達到滿意效果。為了綜合不同方法的優勢，白蘭淑等[33-34]將ISTA方法與線性Bregman方法結合，使聯合方法的前期迭代加快，后期迭代精度提高。由于只采用普通線性Bregman方法及傳統的線性閾值公式，因此重建精度和計算效率仍受到一定限制。

本文在此基礎上，采用Curvelet變換作為稀疏基，引入加速線性Bregman方法(Acceleration linea-rized Bregman Method，ALBM)[35]，并將其與ISTA聯合； 在此過程中，采用軟閾值函數和指數閾值公式，且以新型線性和指數加權因子調節加速Bregman與ISTA方法的比重，充分發揮這兩種算法的優點。理論和實際數據的重建結果表明，本文方法具有比傳統方法更高的重建精度和計算效率。

1 算法原理

1.1 地震數據重建基本原理

缺失道數據的重建過程，可假設為如下數學模型

yobs=Mf

(1)

式中：yobs∈Rm表示觀測到的缺失數據；f∈RM(M≥m)，表示需重建的完整地震數據；M∈Rm×M，表示采樣矩陣。數據重建即是由不完整的觀測數據yobs恢復完整的地震數據f的過程。由于地震數據在時間域不稀疏，本文擬采用具有多尺度多方向特性的曲波變換對地震數據進行稀疏處理，若用C表示曲波變換，x表示稀疏系數，則式(1)可寫成

yobs=MCHx

(2)

式中上標“H”表示共軛轉置矩陣。

由于恢復缺失數據的過程即是求解欠定方程，則可將該問題轉化為求L0范數稀疏約束最小化問題

(3)

由于求解零范數是一個NP-hard問題，則可將L0范數問題轉化為L1范數問題

(4)

由于野外采集數據都含有不同程度的噪聲，為了有效地求解，通常采用正則化思想，在壓縮感知理論框架下，式(4)的解可由以下泛函表達式得到

(5)

式中：A=MCH；λ為正則化參數，它決定L1項與L2項之間的權重。式(5)求解L1范數問題可看作基追蹤問題，要想從隨機缺失地震數據重建滿足后續處理要求的完整數據f，需采用有效反演算法。

1.2 迭代閾值收縮算法

迭代閾值收縮算法在每次迭代過程中都會更新閾值，最小化式(5)中的二次方項； 繼而可通過閾值法投影到L1球上，逐漸收斂到式(5)的解，直到滿足精度要求； 再對N次迭代后的系數x做曲波反變換就可恢復缺失道信息。

閾值迭代法的迭代公式如下

xk+1=T[xk+AH(yobs-Axk)]

k=1,2,…,N

(6)

式中T為閾值函數，分為硬閾值Th函數和軟閾值Ts函數，其公式分別為

(7)

(8)

式中τ表示閾值因子。

1.3 加速線性Bregman算法

線性Bregman方法通過將未閾值化的系數參與下一次迭代來加速收斂。該方法雖然在迭代初期收斂很快，但后期由于在恢復有效信號的同時也會將小于閾值的非有效信號吸收進來，導致最終的重建結果不佳。線性Bregman方法迭代公式為

(9)

為了進一步提高收斂速度，本文在傳統Bregman算法基礎上，引入加速線性Bregman方法(ALBM)[34]，其迭代公式為

(10)

在每次迭代過程中，通過加大未閾值化的曲波系數的比例，在迭代前期保留更大比例的有效信號，進一步提高了初期的收斂速度。

1.4 本文聯合算法

為了同時提高收斂速度和重建精度，將式(6)與式(10)聯合，得到如下迭代式

(11)

式中β為加權因子，它的選取對于提高重建精度和計算效率均具有重要影響。為此，本文引入如下線性加權函數和指數加權函數兩種公式進行對比

(12)

(13)

通過加權因子的作用，可使Bregman加速項zk前期占比較大，提高迭代速度，加速收斂，后期ISTA穩定項xk占比較大，提高迭代的精度，并且使兩種方法過渡得更平滑。

(14)

式中Max為|Cy|的最大值，即稀疏變換系數的最大值。本文閾值參數選擇為指數閾值參數[15]，該參數相比線性閾值參數在保證求解精度的同時收斂速度更快，節省計算工作量，該閾值參數公式為

(15)

2 理論模型

圖1 原始模型采樣及f-k頻譜

為了達到較好處理效果，本次重建模擬采用30次迭代，曲波變換的尺度為4，第二尺度上的方向值為32。圖2a、圖2b分別為閾值迭代法和加速線性Bregman法重建結果，其信噪比分別為14.99、14.87dB，圖2c為傳統線性Bregman法與閾值迭代法聯合(LBM+ISTA)重建結果，信噪比為15.44dB。圖2d為本文加速線性Bregman方法與閾值迭代法聯合(ALBM+ISTA)重建結果，信噪比為16.35dB?？梢娒糠N方法都能將缺失的地震道有效地恢復出來，但與其他方法相比，本文方法重建后的信噪比最高，重建后的有效波同相軸與理論數據一致。

圖3a～圖3d為圖2對應的頻譜圖，可看出有效波能量都得到了明顯收斂，消除了不相干的隨機噪聲，特別從圖3中的紅色矩形可見，相比其他方法，本文方法重建后頻譜與原始數據頻譜最吻合。從相應的誤差圖(圖4)也可看出，本文所提聯合算法重建誤差最小，進一步證明了本文方法的有效性。

圖2 不同算法重建結果

圖3 不同算法重建結果的頻譜

圖4 不同算法重建結果的誤差

為了從細節體現本文方法的有效性，特意抽取某一單道曲線與傳統聯合方法進行對比(圖5)。圖5a和圖5b分別表示第102和178缺失道重建前、后的單道曲線，時窗為0.2～1.8s。對比為原始數據(第1條)，LBM+ ISTA(第2條)和本文方法第3條重建后的單道曲線，以及其誤差單道曲線(第4和第5條)，可見本文方法重建結果與原始數據更吻合，重建后的精度高，誤差小。

圖5 重建結果與誤差的單道顯示

為了進一步比較各種算法在不同迭代次數下的重建信噪比，采用不同迭代次數進行重建，從而得到相應信噪比曲線。圖6a為最大迭代次數為30時的信噪比曲線，整體看初期線性ALBM收斂速度較快，但后期收斂速度變慢，且重建精度較低，而ISTA恰好相反。對于傳統聯合(LBM+ISTA) 方法，收斂速度和重建精度都較好，但整體來看，在相同迭代次數下，本文方法重建精度更高。圖6b為不同迭代次數下的信噪比曲線，可見ALBM在10次以內收斂速度較快，爾后收斂速度變慢，且重建精度也很難進一步提高。ISTA則隨著迭代次數的增加，重建后的信噪比較高，傳統聯合方法吸收了線性Bregman方法和ISTA的優勢，但相比本文所提方法，其收斂速度還是較慢，且可看出本文方法在20次以內就可達到較高信噪比，而后趨于穩定，顯然本文方法更具優勢。

圖6 重建后信噪比曲線

從以上結果還可看出聯合算法結合了加速ALBM和ISTA的優點。前期迭代速度快，這是由ALBM方法完成的，后期ISTA方法使得信噪比得以提高。本文方法在20次迭代后就已經完全收斂，信噪比達到16.29dB，所需時間僅為121.77s。而ISTA和LBM+ISTA均在60次迭代后才開始收斂，約需80次迭代才可完全收斂，所需時間為613s。表1為信噪比達到16dB的迭代次數和所需的時間，可見與ISTA和傳統聯合法相比，在達到滿足后續重建精度信噪比的情況下節省了約3/4的時間。

表1 信噪比達到16dB的迭代次數及時間

為了對比不同加權函數的重建效果，本文選取最大迭代次數5～60，然后記錄不同加權函數的迭代次數與重建信噪比之間的關系曲線。圖7a為本文選取的線性加權閾值和指數加權閾值，這兩個函數的加權值都是1～0，但指數閾值衰減更快，意味著ISTA項的比重快速上升。圖7b為兩種加權因子的不同迭代次數的信噪比曲線，由于兩種加權函數都是1～0，使得信噪比都是由低增至高，且隨著迭代次數的增加，最終信噪比幾乎相同。但由于指數加權函數衰減快，只需較少迭代次數就能使前期加速項衰減更快，在達到較高信噪比之后，穩定項的占比快速上升，最終獲得良好的重建結果。

圖7 不同加權函數重建信噪比曲線

原始地震數據大多含有噪聲，為了檢驗本文方法的抗噪重建能力，對圖1a加入隨機高斯噪聲(圖8a)，再對其進行50%隨機欠采樣(圖8b)，得到本文方法重建結果(圖8c)及其與原始數據的誤差(圖8d)。從重建結果看，該含噪缺失數據的地震波同相軸較連續，有效波信號恢復較好，損失的有效波信息較少，說明本文方法有較強抗噪能力，可滿足實際地震資料處理的要求。

圖8 含噪模型重建結果

3 實際應用

為了檢驗本文方法的實際意義，選取一塊實際原始資料進行處理。圖9a是道距為25m、采樣率為4ms、180道接收的實際單炮記錄，對其進行50%隨機欠采樣和連續5道缺失采樣(圖9b)。設定迭代次數為20次，曲波變換的尺度為4，第二尺度上的方向值為16，信噪比為8.01dB，得到本文方法重建結果(圖10a)及其與原始記錄的誤差剖面(圖10b)?？梢娙笔У佬畔⒌玫接行е亟?，能量較弱的有效波也得到了保護，且在連續缺失5道的情況下，本文方法也能進行有效恢復。圖11是其對應的頻譜，易見本文方法對有效波能量損傷小，效果明顯。

圖9 實際單炮記錄(a)及其隨機欠采樣記錄(b)

圖10 本文方法對實際單炮隨機欠采樣記錄的重建結果(a)及誤差剖面(b)

圖11 實際(a)及其欠采樣(b)和重建(c)地震數據的f-k頻譜

同時，也采用其他算法進行重建并做對比。圖12為不同迭代次數時ALBM、ISTA、LBM+ISTA及本文方法重建后的信噪比與迭代次數關系曲線，可見本文方法在20次迭代后就已開始收斂，而其他方法需40次才開始收斂，迭代60次才可達到與本文方法同樣的重建效果。因此，在實際數據處理中，本文方法在迭代次數較少的情況下能有效地對缺失道進行了重建，提高了地震記錄的信噪比，滿足高精度地震勘探的要求。

圖12 實際數據不同算法迭代次數與信噪比關系

4 討論

前述加速因子的選取非常重要，這也是區別于傳統線性Bregman算法的關鍵所在。通常，該值隨迭代次數增加而增加，變化范圍是1～2。當然，不同的加速因子公式的計算效率不一樣，本文引入傳統加速因子公式，它隨迭代次數的增加而緩慢增大，使得在每次迭代過程中，都可加大未閾值化曲波系數的比例。盡管引入了部分噪聲，但保留了更大比例的有效信號，再通過與閾值迭代法進行聯合，進一步提高加速項前期的比例，可大幅度提高計算速度，而引入的部分噪聲可通過后期閾值迭代法濾除。

本文也嘗試過選擇指數加速因子公式做測試，但從最后的重建精度來看，其信噪比未能提高。這可能是加速因子呈指數快速增大后，會過多地加大未閾值化曲波系數的比例，盡管可提高部分計算精度，但后期閾值迭代法也不能有效濾除未閾值化的曲波系數帶來的大量噪聲，從而導致重建精度未能提高。當然還可選擇其他加速因子計算公式，在保證計算效率的同時，力爭或多或少提高重建精度。

另外，對于曲波變換的尺度和方向值的選擇，通常若尺度選擇較大，方向值選擇較多，則越能反映信號的細節部分，重建精度會有所提高，但計算時間越長。本文綜合考慮計算效率與重建精度，并通過多次測試對比，理論和實際數據尺度值都選擇為4，方向值分別選擇32和16能達到較好效果。

5 結論

本文主要在壓縮感知的基礎上，采用曲波基作為稀疏表示基，在綜合分析ALBM和ISTA兩種算法優點的基礎上，提出了加速聯合迭代算法進行快速高精度數據重建，并采用指數加權因子控制ALBM中的加速項與ISTA中的穩定項的比例，使該加速聯合迭代算法在前期計算中，ALBM起主要作用，加速算法的收斂，后期ISTA起主要作用，大幅度恢復前期所損失的有效波信號，從而提高重建精度，確保最后得到效率快、信噪比高的重建結果。與傳統閾值迭代法、LBM+ISTA方法相比，在達到相同最高信噪比的情況下，本文方法的重建時間節省約70%(表1)，充分突顯了本文方法用時少、精度高的特點。并且在加噪條件下和實際資料處理中也得到了較好重建效果，進一步體現了本文方法的廣泛有效性，可完全滿足海量地震數據處理的要求。