張 華, 楊 會， 刁 塑， 劉 軍

(東華理工大學省部共建核資源與環(huán)境國家重點實驗室培育基地，江西 南昌 330013)

基于指數(shù)閾值迭代法的高精度地震數(shù)據(jù)重建

張 華, 楊 會， 刁 塑， 劉 軍

(東華理工大學省部共建核資源與環(huán)境國家重點實驗室培育基地，江西 南昌 330013)

地震勘探數(shù)據(jù)采集中，地震道缺失是不可能避免的現(xiàn)象。為了滿足后續(xù)處理和解釋的要求，缺失道數(shù)據(jù)重建是地震勘探數(shù)據(jù)處理中必不可少的預處理環(huán)節(jié)。為此，提出一種基于指數(shù)閾值迭代法的高精度重建方法進行疊前重建。引入能夠刻畫地震數(shù)據(jù)局部化特征的多尺度多方向二維曲波變換，采用閾值迭代法進行數(shù)據(jù)重建，并在迭代過程中采用軟閾值算子去除由欠采樣所引起的隨機噪聲干擾。同時在重建過程中針對傳統(tǒng)閾值參數(shù)收斂速度較慢的缺點，提出了一種新的指數(shù)閾值參數(shù)公式，降低計算迭代次數(shù)和提高重建精度。理論數(shù)據(jù)的模擬表明，該方法重建效果顯著，計算速度較快，應(yīng)用于實際地震勘探資料，獲得較好的重建效果。

曲波變換；軟閾值；數(shù)據(jù)重建；閾值迭代法

張華, 楊會，刁塑，等.2017. 基于指數(shù)閾值迭代法的高精度地震數(shù)據(jù)重建[J].東華理工大學學報：自然科學版，40(3)：253-260.

Zhang Hua, Yang Hui, Diao Su, et al.2017. High precision reconstruction of seismic data based on exponential threshold iteration method[J].Journal of East China University of Technology (Natural Science), 40(3):253-260.

在地震勘探中，受復雜地質(zhì)條件和采集環(huán)境的影響，例如障礙物、湖泊、禁采區(qū)以及人工干擾等，導致所采集到地震數(shù)據(jù)通常不完整、不規(guī)則，出現(xiàn)大量的地震缺失道，難以滿足后續(xù)處理以及資料解釋的要求，降低了地震勘探的分辨率(高建軍等，2011；陳曉等，2016；黃光南等，2016)。同時由于項目勘探成本的限制，不可能再到野外對缺失地震道進行補充采集和加密。但是在地下介質(zhì)為連續(xù)性和可預測性的基礎(chǔ)上，在室內(nèi)可以采用高效率高精度的數(shù)據(jù)重建方法重建出缺失的地震道，滿足高分辨率地震勘探的要求。另外一方面，這種高效率高精度的重建方法也可以主動指導復雜地區(qū)的地震數(shù)據(jù)采集，保證整個地震勘探的施工效率以及節(jié)省一定的野外成本。

目前，基于曲波變換的數(shù)據(jù)重建方法研究較多(徐衛(wèi)等，2016；張華等，2013；Zhang et al., 2015)，并且也發(fā)展到了三維數(shù)據(jù)重建，所采用的算法主要有閾值迭代法(Herrmann et al., 2008)、凸集投影算法(Zhang et al., 2013)及譜投影梯度法等(Berg et al., 2008)。一般來講，該類重建方法參數(shù)設(shè)置簡單且不需要地下速度結(jié)構(gòu)信息，計算效率較快，精度高。而對于閾值迭代算法，參數(shù)設(shè)置尤為簡單，并且重建效果也較為顯著，然而由于采用了曲波變換，運算速度仍然較慢，從而使得該方法得不到廣泛的應(yīng)用。實際上，該方法運算速度快慢的關(guān)鍵問題之一在于閾值參數(shù)的選取，然而在以往文章中，無論是采用傅里葉變換還是曲波變換作為稀疏基，都沒有詳細的研究閾值參數(shù)對重建結(jié)果的影響，盡管在這些文章中都保證了在每次迭代過程中閾值不斷減小，然而不同的下降閾值公式其收斂速度不一樣。為此，本文在前人研究的基礎(chǔ)上對不同的閾值參數(shù)公式進行研究，采用曲波變換作為稀疏表示基，提出了新的閾值參數(shù)公式，取得了滿意的重建效果。

1 問題陳述

地震數(shù)據(jù)重建主要通過采樣矩陣得到某未知信號在該矩陣下的線性測量值。為此，假設(shè)如下線性正演模型

yobs=Mf

(1)

其中yobs∈Rn表示不完整的地震數(shù)據(jù)；f∈RN，且N≥n，表示需要重建的完整地震數(shù)據(jù)；M∈Rn×N為采樣矩陣，希望能由不完整的測量數(shù)據(jù)yobs恢復出完整的地震數(shù)據(jù)f。為此，采用多尺度多方向曲波變換對地震數(shù)據(jù)進行稀疏，假設(shè)曲波變換用C表示，則地震數(shù)據(jù)f在曲波域C中的稀疏表示是系數(shù)x，方程(1)可以寫成

yobs=MCHx

(2)

其中H表示共軛轉(zhuǎn)置矩陣，因此，需要采用一定的算法從隨機缺失地震數(shù)據(jù)的稀疏系數(shù)x中重建出無假頻的數(shù)據(jù)f，并且保證重建精度滿足后續(xù)處理的要求。

2 閾值迭代算法

圖1 原始模型采樣及f-k頻譜Fig.1 Original model and the frequency spectruma.理論模型；b.50%隨機欠采樣；c.理論數(shù)據(jù)f-k頻譜；d.隨機采樣數(shù)據(jù)f-k頻譜

要從測量信號yobs中恢復出原始信號的全部信息，直接的辦法就是通過不斷優(yōu)化進行求解l0范數(shù)，而由于上式的求解是一個欠定的病態(tài)問題，不容易得到精確解。但在一定條件下，l1最小范數(shù)和l0最小范數(shù)可以得到同樣的近似解。那么方程(2)就可以變?yōu)閘1最小范數(shù)模型，即：

(3)

實際上，地震數(shù)據(jù)重建的關(guān)鍵技術(shù)之一，就是通過最小化策略求解上述欠定問題。為表述和求解方便，可將其寫為以下形式

(4)

對于方程(4)的解，閾值因子λ的選用特別關(guān)鍵，因為它的作用直接影響到該方程中l(wèi)1范數(shù)和l2范數(shù)兩項之間的權(quán)重。因此，在數(shù)據(jù)重建過程中，閾值參數(shù)λ需要不斷地促進變化，直到最終解滿足一定的精度要求。為此，可以采用閾值迭代法進行求解(Candes et al., 2005)，其迭代式如下：

xn=Sλ(xn-1+AT(y-Axn-1))

(5)

這里Sλ取軟閾值函數(shù)，x初始值可以設(shè)置為零。該式在每次迭代過程中都會更新閾值，最小化(4)式中的二次方項，繼而可通過閾值法投影到l1球上，不斷收斂到方程(4)的解，直到滿足精度要求，再對N次迭代后的系數(shù)xn做曲波反變換就可以恢復缺失道信息。

圖2 線性和指數(shù)閾值重建結(jié)果及f-k頻譜 Fig.2 Reuslts of different threshold parameter and frequency spectruma.線性閾值參數(shù)；b.指數(shù)閾值參數(shù)；c.線性閾值重建后的頻譜d.指數(shù)閾值重建后的頻譜

在解上述非線性問題過程中，首先需要保證閾值較大以強調(diào)稀疏促進的l1項，隨著迭代次數(shù)增加，需要使l2項在求解過程中的占較大影響，為此，使閾值因子λ慢慢遞減，通過不斷的迭代過程，逐漸逼近真實解。因此，閾值參數(shù)的選取工作尤為重要。一般來講，在閾值迭代過程中，閾值參數(shù)λ1滿足，‖Cy‖∞=λ1>λ2>…>ε，式中ε為接近0的小值。常用的線性閾值參數(shù)λi可由下式確定：

(6)

式中Max為|Cy|的最大值，即稀疏變換系數(shù)絕對值的最大值。

采用方程(6)所述的線性閾值公式最終可以重建出缺失道，然而在求解過程中，收斂速度較慢，從而導致該方法不能大規(guī)模工業(yè)生產(chǎn)應(yīng)用。為此，本文提出按指數(shù)λ1規(guī)律衰減的閾值參數(shù)，使得閾值參數(shù)逐漸下降的幅度更大，保證了求解的精度，并且收斂速度更快，節(jié)省計算工作量，該閾值參數(shù)公式為：

3 理論模型

圖3 迭代次數(shù)與信噪比關(guān)系曲線圖Fig.3 The graph of SNR and iteration number after reconstructiona.迭代次數(shù)固定，每次迭代與信噪比圖；b.最大迭代次數(shù)不同時，迭代次數(shù)與信噪比圖

為了分析采用曲波變換后的指數(shù)閾值迭代法重建效果，定義信噪比公式SNR=20log10‖x0‖2/‖x-x0‖2，其中x0表示完整數(shù)據(jù)模型，x表示恢復結(jié)果，該值越高，表明重建精度越高。圖1a為256道合成地震記錄，該記錄總共有4層地震反射波，每一層反射波能量有所差異，采樣間隔為1 ms，道距為4 m，每道1 024個采樣點，然后對其進行隨機欠采樣以模擬野外數(shù)據(jù)缺失過程。圖1b為對模型數(shù)據(jù)進行50%隨機欠采樣后的缺失道地震數(shù)據(jù)剖面圖，圖1c為原始理論數(shù)據(jù)的f-k頻譜，圖1d為欠采樣后的f-k頻譜，從中可以看出隨機欠采樣會產(chǎn)生不相干的隨機噪聲，從而影響到了有效波的振幅能量，必須對其進行疊前重建，恢復缺失地震道能量，提高地震剖面的信噪比，滿足后續(xù)處理的要求。為此，采用曲波變換作為地震信號的稀疏基，引入閾值迭代算法進行求解計算，分別采用上述線性閾值和指數(shù)閾值進行重建。圖2a為采用線性閾值參數(shù)的重建結(jié)果，圖2b為指數(shù)閾值參數(shù)重建的結(jié)果，重建后的信噪比分別為13.22 dB和15.08 dB，其中迭代次數(shù)都為50次，曲波變換的尺度數(shù)為5，方向值為16，圖2c和圖2d分別為其f-k頻譜。從以上結(jié)果可以看出，線性和指數(shù)閾值參數(shù)都可以將缺失道信息有效地恢復出來，并且頻譜能量與原始數(shù)據(jù)非常接近，能量損失較少。但是從計算時間和重建的精度來看，本文所提出的指數(shù)閾值參數(shù)計算效率快，并且重建后的精度高于線性閾值近2 dB。盡管如此，由于隨機欠采樣沒有約束，采樣點完全隨機，造成局部連續(xù)缺失道較多，使得在地震同相軸彎曲較大的區(qū)域則重建精度不高，如圖2b中的第一條反射波同相軸近道部分能量沒有得到有效地收斂，局部誤差較大，這種情況下需要利用其它空間方向的信息進行重建，以彌補單一方向的不足。

為了進一步比較線性閾值參數(shù)和指數(shù)閾值參數(shù)重建后精度以及它們的計算效率，再對圖1b進行重建處理。圖3a表示采用50次迭代時信噪比與迭代次數(shù)關(guān)系曲線圖，它表示了每次迭代后兩種閾值參數(shù)重建結(jié)果的信噪比，從圖中可知指數(shù)閾值參數(shù)在每次迭代過程中重建后的信噪比都高，而線性閾值參數(shù)則信噪比相對較低。然后采用最大迭代次數(shù)5～100次進行分別重建，并且每次迭代次數(shù)都為5的倍數(shù)進行增加，從而計算這兩種閾值參數(shù)分別重建后的信噪比，圖3b為兩種閾值參數(shù)每次最大迭代次數(shù)與信噪比關(guān)系曲線圖，如果想要得到重建后信噪比為13 dB的計算結(jié)果，線性閾值參數(shù)需要經(jīng)過至少迭代50次，從而使得計算速度慢，不能處理海量的地震數(shù)據(jù)。而本文提出的閾值參數(shù)只需迭代18次左右，節(jié)省一大半計算時間，對于處理海量的地震數(shù)據(jù)具有明顯的優(yōu)勢。整體上來講，隨迭代次數(shù)的增大，這兩種閾值參數(shù)重建后的信噪比都會增加，但重建后在相同的信噪比基礎(chǔ)上，迭代次數(shù)太多顯然浪費計算時間，所以，不論從計算時間還是從重建后的信噪比來講，本文所提出的指數(shù)閾值參數(shù)公式具有較大優(yōu)勢。當?shù)螖?shù)超過一定值時，信噪比的增量變化不大，因此為了提高計算效率，本文指數(shù)閾值參數(shù)在隨后的重建處理中采用50次迭代。

圖4 含噪地震數(shù)據(jù)及其重建結(jié)果(50%地震道缺失)Fig.4 Noisy seismic data and reconstruction result (50% missing)a.含噪理論模型；b.50%隨機欠采樣；c.由圖4b重建結(jié)果；d.圖4c與原始含噪模型誤差

由于野外地震數(shù)據(jù)噪聲一般比較發(fā)育，需要檢驗本文方法在噪聲較為發(fā)育下的重建效果。為此，在原始數(shù)據(jù)中加入一定比例的高斯隨機噪聲以模擬野外含噪地震記錄，其加噪結(jié)果如圖4a所示。然后對其進行50%一維隨機欠采樣，從而生產(chǎn)了大量的含噪缺失道，如圖4b所示，然后采用本文高精度指數(shù)閾值迭代法的對該含噪缺失道數(shù)據(jù)進行重建，檢驗在噪聲較為發(fā)育情況下本文方法的重建效果。圖4c為采用本文指數(shù)閾值參數(shù)迭代法的重建結(jié)果，重建后信噪比6.23 dB，含噪缺失地震道的全部信息都得到了較好的恢復，重建后的地震波同相軸非常連續(xù)，信噪比顯著提高。圖4d為重建后的地震剖面與原始加噪理論地震剖面的誤差剖面圖，從誤差剖面圖可以看出該剖面主要為噪聲能量，表明重建前后的有效信息能量損失較小，有效波信號恢復效果較好，說明基于本文所提出的高精度指數(shù)閾值迭代法在含噪地震數(shù)據(jù)重建中具有較強的抗噪能力，滿足處理實際資料的需要。

4 實際應(yīng)用

為了檢驗本文方法的應(yīng)用效果，特意對某區(qū)野外實際資料進行處理，圖5a為實際地震單炮記錄圖，該地震數(shù)據(jù)道距25 m，采樣率4 ms，180道接收。首先對原始單炮記錄進行50%隨機欠采樣，欠采樣結(jié)果如圖5b所示，然后采用本文方法對野外缺失道地震數(shù)據(jù)進行高精度重建，其中迭代次數(shù)為50次，曲波變換的尺度數(shù)為5，方向值為16。圖5c為本文方法重建后的結(jié)果，其信噪比為9.89 dB。在欠采樣50%的情況下，本文方法重建效果較好，缺失道信息得到了有效的恢復，并且弱能量有效波也得到了較好的保護。圖5d為重建后地震記錄與原始單炮記錄的誤差剖面圖，可以看出，除了在近道處誤差相對較大外，其他區(qū)域有效波能量損失較少。為了進一步從頻率域分析重建前后的能量損失情況，特意將其f-k頻譜圖進行顯示(圖6)。從圖6中也可以看出，本文方法重建后f-k譜有效波能量損傷最小，效果顯著。整體上來看，本文方法重建后地震波同相軸非常接近于原始實際地震記錄，提高了地震記錄的信噪比，完全可以滿足高分辨率地震勘探的要求。

圖5 野外數(shù)據(jù)重建過程(50%地震道缺失)Fig.5 Reconstruction result of the field data(50% missing) a.實際單炮記錄圖；b.隨機欠采樣記錄圖；c.本文方法重建結(jié)果；d.重建誤差剖面

5 結(jié)論

本文采用曲波變換作為稀疏表示基，引入閾值迭代算法進行重建，取得了較好的重建效果。通過研究可知，在得到同樣的重建精度時，傳統(tǒng)線性閾值迭代的重建方法效果有限，所需計算時間較多。而本文采用指數(shù)閾值迭代的進行數(shù)據(jù)重建，則可以得到了高精度的重建效果。通過理論對比分析表明本文方法在相同迭代次數(shù)下能夠得到更高的信噪比。而重建后達到相同的信噪比時，則所需要的迭代次數(shù)較少，節(jié)約了計算時間，從而使得本文方法具有參數(shù)設(shè)置簡單、用時少、精度高等特點，并且在實際資料的應(yīng)用中也得到了較好的重建效果。

從本文的重建結(jié)果可以看出，由于曲波變換能夠反映地震信號的局部特征，更加稀疏的表示地震信號的波前特征，可以得到更好精度的重建結(jié)果，但是曲波變換冗余度高，與傅立葉變換相比，運算速度較慢，盡管本文提出了指數(shù)閾值參數(shù)公式，但在處理大量地震數(shù)據(jù)時還是會受到速度限制，因此也需要進一步優(yōu)化算法，以達到工業(yè)化生產(chǎn)的需要。

圖6 重建過程中的地震數(shù)據(jù)f-k頻譜圖Fig.6 The f-k spectrum of reconstruction processesa～c分別為圖5a～5c的f-k頻譜圖

陳曉, 周磊, 于鵬, 等. 2016. 基于模擬退火算法的大地電磁和地震數(shù)據(jù)貝葉斯聯(lián)合反演[J].東華理工大學學報：自然科學版,39(1):59-66.

高建軍, 陳小宏, 李景葉. 2011. 三維不規(guī)則地震數(shù)據(jù)重建方法[J]. 石油地球物理勘探,46(1):40-46.

黃光南，鄧居智,李紅星, 等. 2016. 三維VSP地震初至波走時層析成像研究[J]. 東華理工大學學報：自然科學版,39(3):266-272.

徐衛(wèi), 張華, 張落毅. 2016. 基于復值曲波變換的地震數(shù)據(jù)重建方法[J]. 物探與化探,40(4):750-755.

張華, 陳小宏. 2013. 基于jitter采樣和曲波變換的三維地震數(shù)據(jù)重建[J].地球物理學報,56(5):1637-1649.

Berg E D，F(xiàn)riedlander M P．2008. Probing the Pareto frontier for basis pursuit solutions[J]. SIAM J．on Scientific Computing，31(2):890-912.

Candes E J, Romberg J. 2005. Practical Signal Recovery from Random projections[J]. Proceedings of SPIE Computational Imaging III, 23(2):76-86.

Herrmann F J, Hennenfent G. 2008. 1Non-parametric seismic data recovery with curvelet frames[J]. Geophysical Journal International, 73(1):233-248.

Zhang H, Chen X H, Li H X. 2015. 3D seismic data reconstruction based on complex-valued curvelet transform in frequency domain[J]. Journal of Applied Geophysics, 113(6):64-73.

Zhang H, Chen X H, Wu X M. 2013. Seismic data reconstruction based on CS and Fourier theory[J]. Applied Geophysics, 10(3):170-180.

HighPrecisionReconstructionofSeismicDataBasedonExponentialThresholdIterationMethod

ZHANG Hua, YANG Hui, DIAO Su, LIU Jun

(State Key Laboratory Breeding Base of Nuclear Resources and Environment，East China University of Technology，Nanchang，JX 330013, China)

During the data acquisition of seismic exploration, the phenomenon of missing traces are unavoidable. In order to meet the requirements for subsequent processing and interpretation, missing data reconstruction is essential preprocessing steps in any seismic data processing chain. The data reconstruction method based on the exponential threshold iterative method has been introduced in the paper. Firstly, multi-scale and multi-directional curvelet transform to characterize the local features of seismic data has been introduced and the threshold iterative method is adopted. Meanwhile, a soft thresholding is introduced to remove the random noise arisen by the undersampling. Aiming at the disadvantage of slow convergence in traditional threshold parameter during the reconstruction process, an new exponential decreased threshold is also proposed and it can reduce iterations and improve reconstruction efficiency. The simulation results on synthetic seismic data showed that this method has a better effect and faster computation. At last, we apply this technology into real seismic data and obtain a good result.

curvelet transform; soft Threshold; data reconstruction; threshold iterative method

P631

A

1674-3504(2017)03-0253-08

2016-12-07

國家自然科學基金(41304097，41664006)；江西省自然科學基金(20151BAB203044, 20171BAB203031, 20171BAB202028)；江西省杰出青年人才資助計劃(20171BCB23068)

張 華( 1979—),男,博士,副教授,研究方向為數(shù)據(jù)重建及規(guī)則化反演。E-mail: zhhua1979@163.com

10.3969/j.issn.1674-3504.2017.03.006