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

張 華, 楊 會(huì)， 刁 塑， 劉 軍

(東華理工大學(xué)省部共建核資源與環(huán)境國(guó)家重點(diǎn)實(shí)驗(yàn)室培育基地，江西 南昌 330013)

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

張 華, 楊 會(huì)， 刁 塑， 劉 軍

(東華理工大學(xué)省部共建核資源與環(huán)境國(guó)家重點(diǎn)實(shí)驗(yàn)室培育基地，江西 南昌 330013)

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

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

張華, 楊會(huì)，刁塑，等.2017. 基于指數(shù)閾值迭代法的高精度地震數(shù)據(jù)重建[J].東華理工大學(xué)學(xué)報(bào)：自然科學(xué)版，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.

在地震勘探中，受復(fù)雜地質(zhì)條件和采集環(huán)境的影響，例如障礙物、湖泊、禁采區(qū)以及人工干擾等，導(dǎo)致所采集到地震數(shù)據(jù)通常不完整、不規(guī)則，出現(xiàn)大量的地震缺失道，難以滿足后續(xù)處理以及資料解釋的要求，降低了地震勘探的分辨率(高建軍等，2011；陳曉等，2016；黃光南等，2016)。同時(shí)由于項(xiàng)目勘探成本的限制，不可能再到野外對(duì)缺失地震道進(jìn)行補(bǔ)充采集和加密。但是在地下介質(zhì)為連續(xù)性和可預(yù)測(cè)性的基礎(chǔ)上，在室內(nèi)可以采用高效率高精度的數(shù)據(jù)重建方法重建出缺失的地震道，滿足高分辨率地震勘探的要求。另外一方面，這種高效率高精度的重建方法也可以主動(dòng)指導(dǎo)復(fù)雜地區(qū)的地震數(shù)據(jù)采集，保證整個(gè)地震勘探的施工效率以及節(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)。一般來(lái)講，該類重建方法參數(shù)設(shè)置簡(jiǎn)單且不需要地下速度結(jié)構(gòu)信息，計(jì)算效率較快，精度高。而對(duì)于閾值迭代算法，參數(shù)設(shè)置尤為簡(jiǎn)單，并且重建效果也較為顯著，然而由于采用了曲波變換，運(yùn)算速度仍然較慢，從而使得該方法得不到廣泛的應(yīng)用。實(shí)際上，該方法運(yùn)算速度快慢的關(guān)鍵問(wèn)題之一在于閾值參數(shù)的選取，然而在以往文章中，無(wú)論是采用傅里葉變換還是曲波變換作為稀疏基，都沒(méi)有詳細(xì)的研究閾值參數(shù)對(duì)重建結(jié)果的影響，盡管在這些文章中都保證了在每次迭代過(guò)程中閾值不斷減小，然而不同的下降閾值公式其收斂速度不一樣。為此，本文在前人研究的基礎(chǔ)上對(duì)不同的閾值參數(shù)公式進(jìn)行研究，采用曲波變換作為稀疏表示基，提出了新的閾值參數(shù)公式，取得了滿意的重建效果。

1 問(wèn)題陳述

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

yobs=Mf

(1)

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

yobs=MCHx

(2)

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

2 閾值迭代算法

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

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

(3)

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

(4)

對(duì)于方程(4)的解，閾值因子λ的選用特別關(guān)鍵，因?yàn)樗淖饔弥苯佑绊懙皆摲匠讨衛(wèi)1范數(shù)和l2范數(shù)兩項(xiàng)之間的權(quán)重。因此，在數(shù)據(jù)重建過(guò)程中，閾值參數(shù)λ需要不斷地促進(jìn)變化，直到最終解滿足一定的精度要求。為此，可以采用閾值迭代法進(jìn)行求解(Candes et al., 2005)，其迭代式如下：

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

(5)

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

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

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

(6)

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

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

3 理論模型

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

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

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

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

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

4 實(shí)際應(yīng)用

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

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

5 結(jié)論

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

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

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

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

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

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

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

張華, 陳小宏. 2013. 基于jitter采樣和曲波變換的三維地震數(shù)據(jù)重建[J].地球物理學(xué)報(bào),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

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

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

10.3969/j.issn.1674-3504.2017.03.006