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

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

0 引言

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

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

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

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

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

1 算法原理

1.1 地震數(shù)據(jù)重建基本原理

缺失道數(shù)據(jù)的重建過程，可假設(shè)為如下數(shù)學(xué)模型

yobs=Mf

(1)

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

yobs=MCHx

(2)

式中上標(biāo)“H”表示共軛轉(zhuǎn)置矩陣。

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

(3)

由于求解零范數(shù)是一個NP-hard問題，則可將L0范數(shù)問題轉(zhuǎn)化為L1范數(shù)問題

(4)

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

(5)

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

1.2 迭代閾值收縮算法

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

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

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

k=1,2,…,N

(6)

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

(7)

(8)

式中τ表示閾值因子。

1.3 加速線性Bregman算法

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

(9)

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

(10)

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

1.4 本文聯(lián)合算法

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

(11)

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

(12)

(13)

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

(14)

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

(15)

2 理論模型

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

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

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

圖2 不同算法重建結(jié)果

圖3 不同算法重建結(jié)果的頻譜

圖4 不同算法重建結(jié)果的誤差

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

圖5 重建結(jié)果與誤差的單道顯示

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

圖6 重建后信噪比曲線

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

表1 信噪比達(dá)到16dB的迭代次數(shù)及時間

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

圖7 不同加權(quán)函數(shù)重建信噪比曲線

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

圖8 含噪模型重建結(jié)果

3 實際應(yīng)用

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

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

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

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

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

圖12 實際數(shù)據(jù)不同算法迭代次數(shù)與信噪比關(guān)系

4 討論

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

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

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

5 結(jié)論

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