基于L_1/2正則化的拋物線(xiàn)Radon變換多次波壓制方法

吳秋瑩 胡斌 劉財(cái) 高銳

摘要：在地震數(shù)據(jù)處理中，多次波的存在會(huì)對(duì)地震數(shù)據(jù)成像和地震資料解釋帶來(lái)影響，如何有效地壓制多次波干擾是地震勘探中的重要問(wèn)題。拋物線(xiàn)Radon變換因其高效的特點(diǎn)被廣泛應(yīng)用于多次波壓制中，但在野外地震數(shù)據(jù)采集時(shí)，炮檢距的有限性會(huì)導(dǎo)致變換域中的能量擴(kuò)散，產(chǎn)生假象，使多次波壓制達(dá)不到理想的效果。針對(duì)此問(wèn)題，提出一種基于L_1/2正則化的稀疏反演高分辨拋物線(xiàn)Radon變換，并應(yīng)用廣義迭代收縮算法（generalized iterated shrinkage algorithm， GISA）進(jìn)行求解。研究結(jié)果表明，L_1/2正則化有很強(qiáng)的稀疏約束能力，能提高解的稀疏度，改進(jìn)信噪分離的效果。與最小二乘反演和基于L₁正則化的稀疏反演相比，基于L_1/2正則化的稀疏反演高分辨拋物線(xiàn)Radon變換能更有效地壓制多次波，并確保了重構(gòu)數(shù)據(jù)與原始數(shù)據(jù)的一致性。

關(guān)鍵詞：多次波壓制；高分辨率拋物線(xiàn)Radon變換；L_1/2正則化

doi：10.13278/j.cnki.jjuese.20220307

中圖分類(lèi)號(hào)：P631.4

文獻(xiàn)標(biāo)志碼：A

收稿日期：2022-11-08

作者簡(jiǎn)介：吳秋瑩（1993—），女，博士研究生，主要從事地震數(shù)據(jù)處理工作，E-mail： qywu20@mails.jlu.edu.cn

通信作者：劉財(cái)（1963—），男，教授，博士生導(dǎo)師，主要從事地震波場(chǎng)正反演理論、綜合地球物理等方面研究，E-mail： liucai@jlu.edu.cn

基金項(xiàng)目：國(guó)家自然科學(xué)基金項(xiàng)目（41874125）

Supported by the National Natural Science Foundation of China （41874125）

Multiple Suppression Method of Parabolic Radon Transform

Based on L_1/2RegularizationWu Qiuying ¹， Hu Bin¹， Liu Cai¹， Gao Rui^2，3

1. College of GeoExploration Science and Technology， Jilin University， Changchun 130026， China

2. School of Earth Sciences and Engineering， Sun Yat-Sen University， Guangzhou 510275， China

3. Key Laboratory of Deep-Earth Dynamics of Ministry of Natural Resources， Institute of Geology， Chinese Academy of

Geological Sciences， Beijing 100037， China

Abstract： In the context of seismic data processing， the presence of multiples poses inherent challenges to the imaging and interpretation of seismic data. The effective suppression of these multiples stands as a key issue in seismic exploration. Leveraging its high efficiency， the parabolic Radon transform emerges as a widely used technique for multiple suppression. However， in field seismic data acqisition， due to the limited offset， energy diffusion and illusions reduce the effect of multiple suppression in the Radon domain. In response to this challenge， we propose a L_1/2-regularized high-resolution parabolic Radon transform with sparse inversion， where the inverse problem is solved by? generalized iterated shrinkage algorithm （GISA）. The L_1/2 regularization chosen for its robust sparse constraint capabilities plays an important role in enhancing the solution sparsity and improving the signal-noise separation. Compared with the least square inversion and the sparse inversion method based on L₁regularization， the L_1/2-regularized sparse inversion of using the high-resolution parabolic Radon transform can suppress multiples effectively and ensure the consistency between the reconstructed data and the original data.

Key words： multiple suppression; high-resolution parabolic Radon transform; L_1/2regularization

0 引言

多次波在反射波地震勘探中通常被視為干擾波，較強(qiáng)的多次波會(huì)影響深層一次波的成像效果，在地震資料解釋中造成誤解。因此，針對(duì)壓制多次波的問(wèn)題，國(guó)內(nèi)外學(xué)者對(duì)此做了大量的研究。壓制多次波的方法一般可分為兩個(gè)類(lèi)別^[1-2]：一是濾波類(lèi)方法，這類(lèi)方法基于一次波和多次波的周期性、空間特性等差異來(lái)壓制多次波，常用的方法有預(yù)測(cè)反褶積^[3]、Radon域?yàn)V波^[4]、f-k（頻率-波數(shù)）域?yàn)V波^[5]、聚束濾波^[6]等；二是預(yù)測(cè)相減類(lèi)方法，這類(lèi)方法是對(duì)多次波進(jìn)行預(yù)測(cè)并從地震記錄中減去，常用的方法有逆散射級(jí)數(shù)法^[7]、波場(chǎng)延拓法^[8]、自由表面多次波壓制法^[9-10]等。由于一次反射波與具有相同法線(xiàn)時(shí)間的多次反射波相比，傳播速度更高，兩者同相軸之間存在明顯差異，當(dāng)一次波和多次波間的速度差異較大時(shí)，濾波方法能有效且高效地壓制多次波。由于一次反射和多次反射在τ-q（截距-曲率）域中的分離效果較f-k、τ-p（斜率）域中的分離效果更好，且更高效、成本更低，故拋物線(xiàn)Radon變換被廣泛地應(yīng)用于壓制多次波。高分辨率拋物線(xiàn)Radon變換能使一次波和多次波在速度差異較小時(shí)也得以分辨。因此，本文研究了拋物線(xiàn)Radon變換及其高分辨率方法。

Radon變換于1917年由Radon提出^[11]，20世紀(jì)70年代被引入到地球物理領(lǐng)域。1986年，Hampson^[4]利用拋物線(xiàn)Radon變換在動(dòng)校正后的共中心點(diǎn)（common middle point， CMP）道集上壓制表面多次波，并提出利用最小二乘反演法來(lái)優(yōu)化Radon域的數(shù)據(jù)，此反演方法減少了拋物線(xiàn)Radon域的模糊現(xiàn)象。然而影響Radon變換分辨率的主要因素包括空間假頻、截?cái)嗾`差、保幅性問(wèn)題等，通常可以應(yīng)用稀疏反演方法來(lái)解決地震數(shù)據(jù)的稀疏性問(wèn)題以提高Radon域中的數(shù)據(jù)分辨率。所以，1995年，Sacchi等^[12]將Radon變換考慮為稀疏反演過(guò)程，提高了Radon變換域的分辨率，并將拋物線(xiàn)Radon變換在頻率域?qū)嵭?，提高Radon變換的效率。2000年，Herrmann等^[13]利用一種高分辨的拋物線(xiàn)Radon變換區(qū)分一次波和多次波，以非迭代的方式處理混疊和分辨率問(wèn)題，使空間假頻的問(wèn)題得到一定的解決。2009年，熊登等^[14]提出利用混合域的Radon變換壓制多次波，提高了計(jì)算效率和分辨率。2013年，Li等^[15]利用λ-f域（其中，λ=fq）的Radon變換對(duì)多次波壓制后的數(shù)據(jù)進(jìn)行保幅。2016年，鞏向博等^[16]提出混合域雙曲線(xiàn)Radon變換多次波壓制方法，提高了分辨率和計(jì)算效率。

本文重點(diǎn)研究基于稀疏反演的高分辨率Radon變換。稀疏反演問(wèn)題在求解時(shí)，會(huì)存在解不唯一的情況，為解決此問(wèn)題，可以在高分辨率Radon變換基礎(chǔ)上加入L范數(shù)約束的稀疏約束條件，得到稀疏約束的Radon變換。L₀范數(shù)正則化約束具有稀疏性，但其最優(yōu)化問(wèn)題求解較困難，可以使用L₁范數(shù)替代^[17]。L₁正則化是L₀正則化的最優(yōu)凸近似，它比L₀正則化容易求解并且更具稀疏性，但已有研究表明L₁正則化并不是最佳的稀疏正則化^[18]。2010年，Xu等^[19]提出了L_1/2正則化，證明了L_1/2正則化可以得到比L₁正則化更稀疏、更穩(wěn)定的解，是L_p（01/2正則化應(yīng)用于各種研究中，如2012年Zeng等^[20]用于稀疏雷達(dá)成像中、2013年Sun等^[21]用于遙感數(shù)據(jù)高光譜解混中、2019年康治梁等^[22]用于地震稀疏反褶積中。

本文將L_1/2正則化方法引入拋物線(xiàn)Radon變換壓制多次波中，并與L₁正則化方法進(jìn)行對(duì)比，通過(guò)模型數(shù)據(jù)及實(shí)際數(shù)據(jù)驗(yàn)證方法的有效性。

1 基本理論

1.1 拋物線(xiàn)Radon變換

Radon變換多次波壓制方法基于一次波和多次波的時(shí)差差異，在采用介于一次波與多次波之間的速度對(duì)地震數(shù)據(jù)的CMP道集進(jìn)行動(dòng)校正后，同相軸的剩余時(shí)差接近拋物線(xiàn)形狀，因此經(jīng)過(guò)拋物線(xiàn)Radon變換，可以將一次波和多次波具有不同曲率的同相軸映射到Radon域中的不同區(qū)域，以實(shí)現(xiàn)一次波與多次波的分離。拋物線(xiàn)Radon變換的求和軌跡沿著拋物線(xiàn)進(jìn)行，其正反變換的表達(dá)式分別為：

式中：d（t，x_k）為地震數(shù)據(jù)；m（τ，q_j）為拋物線(xiàn)Radon變換后的Radon域數(shù)據(jù)；x_k為偏移距；N_x為地震道數(shù)；t為時(shí)間；q_j為曲率參數(shù)；N_q為曲率參數(shù)的數(shù)量。拋物線(xiàn)Radon變換使在時(shí)-空域具有拋物線(xiàn)形狀的同相軸在拋物線(xiàn)Radon域映射成一個(gè)點(diǎn)。

為提高拋物線(xiàn)Radon變換的計(jì)算效率，通常對(duì)地震剖面沿時(shí)間方向進(jìn)行Fourier變換，其正反變換公式分別為：

式中：f為頻率；M（q_j，f）為頻率域Radon數(shù)據(jù)；D（x_k，f）為頻率域地震數(shù)據(jù)。利用式（3）使每個(gè)頻率成分分別反演，獨(dú)立進(jìn)行Radon變換。式（3）和（4）表示為矩陣形式可分別描述為：

m=A^Hd；??? （5）

d=Am。??? （6）

其中，算子A及其共軛轉(zhuǎn)置A^H分別定義為：

Hampson^[4]采用拋物線(xiàn)Radon變換在CMP道集中壓制多次波，并且利用最小二乘反演法，得到最小二乘解為

m=（A^HA）^-1A^Hd。 ???（9）

通常，為了反演過(guò)程的穩(wěn)定性，在矩陣主對(duì)角線(xiàn)上引入阻尼參數(shù)α，即

m=（A^HA+αI）^-1A^Hd。 ???（10）

α一般選擇為矩陣A^HA的百分之幾。

最小二乘拋物線(xiàn)Radon變換提高了Radon域數(shù)據(jù)的分辨率以及變換的精度和聚焦能力，使Radon域假象減少，限制部分能量擴(kuò)散，且確保了重構(gòu)數(shù)據(jù)與原始數(shù)據(jù)的一致性，然而，結(jié)果仍然不是最佳的。

影響Radon變換分辨率的主要因素有以下幾點(diǎn)^[23]：一是q的取值范圍。若q的取值范圍較大，則Radon域的能量不能較好地聚焦，產(chǎn)生空間假頻問(wèn)題。q的變化量Δq和最大值q_max應(yīng)滿(mǎn)足^[24]：

式中：x_max為最大跑檢距；f_max為信號(hào)的最高頻率。二是由于炮檢距的有限性，能量出現(xiàn)發(fā)散現(xiàn)象，使Radon域中的映射不是一個(gè)點(diǎn)，而是具有拖尾現(xiàn)象的剪刀狀。因此，希望Radon域的解更加稀疏，使x-t域的同相軸在Radon域投影成一個(gè)帶限的脈沖，應(yīng)用稀疏反演方法解決地震數(shù)據(jù)的稀疏性問(wèn)題，而后優(yōu)化反演中的稀疏化參數(shù)，使同相軸在Radon域有更好的局部化特征，對(duì)反演結(jié)果不斷地進(jìn)行迭代修改，并優(yōu)化求解算法，提高Radon域中的數(shù)據(jù)分辨率。三是噪聲。應(yīng)在保證有效信號(hào)不丟失的情況下盡可能地壓制噪聲，故對(duì)信號(hào)進(jìn)行保幅性處理。

1.2 基于L_1/2正則化的拋物線(xiàn)Radon變換

利用L范數(shù)作為懲罰項(xiàng)來(lái)約束反演目標(biāo)函數(shù)的正則化方法，可以使結(jié)果滿(mǎn)足稀疏化，保留待反演參數(shù)的特征信息，并且求得唯一的稀疏解，從而獲得稀疏約束的高分辨率拋物線(xiàn)Radon變換。常用的正則化方法有L₀正則化、L₁正則化、L₂正則化、L_p正則化，用單位球幾何形狀表示如圖1所示。L₀范數(shù)表示向量中非零元素的數(shù)目，其正則化約束區(qū)域沒(méi)有角（圖1a），很難出現(xiàn)零解，所以該方法具有稀疏性，但其求解困難，應(yīng)用較少；L₁正則化的約束區(qū)域是正方形（圖1b），該正則化具有對(duì)異常值不敏感的特性，是凸優(yōu)化問(wèn)題，可以使用增廣拉格朗日乘子法^[25]、迭代軟閾值法^[26]、迭代重加權(quán)最小二乘法^[27]、梯度投影法^[28]等進(jìn)行求解，能實(shí)現(xiàn)較好的稀疏效果；L₂正則化可以避免模型發(fā)生過(guò)擬合現(xiàn)象，但其約束區(qū)域沒(méi)有角（圖1c），很難出現(xiàn)零解，可以避免模型發(fā)生過(guò)擬合現(xiàn)象，但沒(méi)有L₁正則化稀疏。

Chartrand等^[18]在壓縮感知方面對(duì)L₁正則化優(yōu)化問(wèn)題進(jìn)行研究的實(shí)驗(yàn)數(shù)據(jù)表明，L_p（01正則化能從更少的觀測(cè)數(shù)據(jù)中恢復(fù)稀疏信號(hào)，信號(hào)恢復(fù)能力更強(qiáng)、更穩(wěn)定、有更低的預(yù)測(cè)誤差。Xu等^[19]系統(tǒng)研究了L_p（01/2正則化問(wèn)題。如圖1d所示，L_1/2正則化的解出現(xiàn)在拐角處，表明比L₁正則化更稀疏。研究表明，當(dāng)1/2≤p<1時(shí)，L_1/2正則化是L_p正則化中最稀疏的；當(dāng)0p正則化具有與L_1/2正則化相似的性質(zhì)，p的取值對(duì)稀疏解的影響不大^[19]。所以說(shuō)L_1/2正則化是L_p（0

L_1/2正則化的最優(yōu)化問(wèn)題可表示為

L_p（0

根據(jù)GISA，式（14）的迭代解可表示為：

m_k+1=GST（θ_k+1，ξμ，p）; ???（15）

θ_k+1=m_k-μA^H（Am_k-d）。 ???（16）

式中：m_k為迭代k次的迭代解；ξ為迭代步長(zhǎng)；GST（·）為廣義軟閾值算子，定義為

式中，sgn（·）為符號(hào)函數(shù)。

為了加速算法的收斂性，采用Beck等^[32]提出的快速收斂策略，將Nesterov加速算法引入到GISA中：

此方法提高了反演的收斂速度以及計(jì)算效率。最終得到m︿，獲得了拋物線(xiàn)Radon變換更具稀疏性、且高分辨率的結(jié)果，達(dá)到理想的多次波壓制效果。

2 正演模型測(cè)試

為驗(yàn)證以上理論，利用模擬地震記錄進(jìn)行模型測(cè)試。模擬數(shù)據(jù)理論模型為八層水平層狀速度模型（圖2）。時(shí)間方向有500個(gè)采樣點(diǎn)，采樣間隔為4 ms，空間方向上為126道，偏移距為0～3 150 m，道間距25 m。圖3為動(dòng)校正前后的CMP道集。從原始CMP道集（圖3a）中可以觀察到多條一次波和一系列多次波，多次波集中在中下部0.6～2.0 s。從經(jīng)過(guò)動(dòng)校正后的CMP道集（圖3b）中可以看出，由于多次波的速度較相同法線(xiàn)時(shí)間一次波的速度更低，所以在CMP道集中按照一次波速度動(dòng)校正后，一次波被拉平，而多次波因校正不足，呈現(xiàn)下拉的狀態(tài)，使多次波同相軸剩余時(shí)差更接近拋物線(xiàn)形態(tài)。在圖3b中觀察到的下拉狀態(tài)的同相軸均為多次波。

利用傳統(tǒng)最小二乘拋物線(xiàn)Radon變換方法和基于L₁正則化的拋物線(xiàn)Radon變換方法作為對(duì)比組，以驗(yàn)證本文方法在多次波壓制方面的效果。從圖4可以看出，經(jīng)過(guò)拋物線(xiàn)Radon變換后，一次波分布在q=0附近和q＜0區(qū)域，多次波分布在q＞0區(qū)域。利用最小二乘拋物線(xiàn)Radon變換得到的Radon域數(shù)據(jù)，同相軸在Radon域中沒(méi)有較好地聚成一點(diǎn)，產(chǎn)生較多假象（圖4a）；利用基于L₁正則化拋物線(xiàn)Radon變換得到的Radon域數(shù)據(jù)假象減少，但效果并不理想（圖4b）；使用基于L_1/2正則化拋物線(xiàn)Radon變換得到的Radon域數(shù)據(jù)（圖4c），比圖4b的“拖尾”效應(yīng)得到明顯的改善，能量更匯聚，數(shù)據(jù)變換精度和聚焦能力得到進(jìn)一步提高。

接下來(lái)，在Radon域中切除多次波所在的q＞0區(qū)域，再變換回時(shí)-空域，得到只有一次波的動(dòng)校正后的CMP道集，再經(jīng)過(guò)反動(dòng)校正，得到壓制多次波后的CMP道集。

將多次波壓制結(jié)果與不含多次波的波場(chǎng)進(jìn)行對(duì)比，結(jié)果如圖5所示?？梢钥闯觯鹤钚《藪佄锞€(xiàn)Radon變換（圖5b）和基于L₁正則化的拋物線(xiàn)Radon變換（圖5c）的多次波壓制結(jié)果在0.7、1.2以及1.6 s處仍有一些殘余的多次波能量；本文方法的多次波壓制結(jié)果在相同位置幾乎沒(méi)有殘余的多次波（圖5d）。

通過(guò)速度譜對(duì)結(jié)果進(jìn)行進(jìn)一步對(duì)比。圖6為本文方法壓制多次波前后的速度譜，從中可以明顯看到，時(shí)間在0.8～2.0 s、速度在2 000 m/s的多次波能量得到了較好的壓制，而一次波的能量得以保留，利于更準(zhǔn)確地識(shí)別速度信息。

綜上，在基于3種算法的拋物線(xiàn)Radon變換壓制多次波中，本文方法能獲得最稀疏的解，使多次波得到極好的壓制，并且也能有效保護(hù)一次波的振幅，具有一定程度的保幅性。

為測(cè)試本文方法的抗噪能力，以S_NR=m/（m-m-）作為信噪比（m為不含噪聲信號(hào)，m-為含噪聲信號(hào)）。分別以信噪比為5、-5、-15 dB加入到某一模型數(shù)據(jù)中，得到不同信噪比的正演模型數(shù)據(jù)。利用本文方法對(duì)數(shù)據(jù)進(jìn)行重建，結(jié)果如圖7所示，不同信噪比數(shù)據(jù)的重構(gòu)誤差如表1所示。

由圖7、表1可知，即使噪聲較強(qiáng)，對(duì)重建數(shù)據(jù)也并沒(méi)有過(guò)大影響；可以證明本文方法的抗噪能力較強(qiáng)，能夠?qū)崿F(xiàn)具有抗噪性能的多次波壓制。

3 實(shí)際數(shù)據(jù)處理

為進(jìn)一步驗(yàn)證本文方法，采用某海上實(shí)際數(shù)據(jù)動(dòng)校正后的CMP道集（圖8）進(jìn)行多次波壓制處理。時(shí)間方向上有1 305個(gè)采樣點(diǎn)，截取3.2～4.8 s多次波較多的部分進(jìn)行多次波壓制，該部分?jǐn)?shù)據(jù)在時(shí)間方向上有400個(gè)采樣點(diǎn)，采樣間隔為4 ms，空間方向上為92道，道間距為175 m，偏移距為16 100 m。實(shí)際數(shù)據(jù)中較深層位置3.8～4.8 s處存在大量的多次波，導(dǎo)致一次波數(shù)據(jù)被掩蓋，影響一次波成像效果。

圖9為利用3種算法得到的Radon域數(shù)據(jù)，可以看出，本文方法處理的Radon域數(shù)據(jù)（圖9c）與最小二乘拋物線(xiàn)Radon變換（圖9a）和基于L₁正則化的拋物線(xiàn)Radon變換（圖9b）相比具有更高的分辨率，能量聚焦明顯且假象明顯較少。

圖10為多次波壓制結(jié)果。將3種方法的多次波壓制結(jié)果（圖10b、c、d）與動(dòng)校正后的原始CMP道集（圖10a）進(jìn)行對(duì)比，可以看到本文方法（圖10d）的壓制結(jié)果優(yōu)于最小二乘拋物線(xiàn)Radon變換（圖10b）和基于L₁正則化的拋物線(xiàn)Radon變換（圖10c）：在3.8～4.2 s和4.4～4.7 s處，圖10b和圖10c明顯有大量的多次波殘留，而且由于近偏移距處的多次波和一次波速度相近，在Radon域切除q＞0區(qū)域時(shí)，也會(huì)影響一次波的投影區(qū)域，使近偏移距損失了部分一次波能量；而本文方法的一次波能量得到了一定的保護(hù)（圖10d），并且在3.8～4.2 s和4.4～4.7 s處沒(méi)有明顯的多次波殘留。

圖11a、b為分別利用基于L₁正則化和L_1/2正則化Radon變換后又經(jīng)Radon反變換重構(gòu)的地震數(shù)據(jù)，對(duì)比圖11a、b可以直觀地看出，本文方法的重構(gòu)數(shù)據(jù)連續(xù)性更強(qiáng)，且與原始數(shù)據(jù)（圖10a）具有更強(qiáng)的一致性。

利用重構(gòu)誤差定量分析重構(gòu)數(shù)據(jù)^[33]與原始數(shù)據(jù)的一致性：

式中：m為原始數(shù)據(jù)；m′為重構(gòu)數(shù)據(jù)；s為重構(gòu)誤差。經(jīng)計(jì)算，求得基于L₁正則化和L_1/2正則化的拋物線(xiàn)Radon變換重構(gòu)誤差分別為：s_L1=15%，s_L1/2=8%。實(shí)際數(shù)據(jù)的波場(chǎng)較為復(fù)雜，含有其他噪聲，本文方法的重構(gòu)數(shù)據(jù)比基于L₁正則化拋物線(xiàn)Radon變換的方法更逼近原始數(shù)據(jù)。

圖12為壓制多次波前后3種方法的單道對(duì)比，可以進(jìn)一步對(duì)比多次波壓制效果。截取了3.7～4.8 s處多次波較豐富的部分，可以直觀地看到，利較小，直觀地證明了本文方法的有效性。

4 討論

隨著三維地震勘探技術(shù)的日益發(fā)展，三維地震數(shù)據(jù)處理也被廣泛關(guān)注。為驗(yàn)證本文方法對(duì)三維資料的實(shí)用性，我們模擬生成了三維地震數(shù)據(jù)的CMP道集?？v測(cè)線(xiàn)偏移距范圍為-5 000～5 000 m，聯(lián)絡(luò)測(cè)線(xiàn)方向偏移距范圍為-2 500～2 500 m，多次波剩余時(shí)差分別為600、400、200、60 ms，振幅分別為-1.0、0.5、0.1、-0.9。圖13a為從三維模擬數(shù)據(jù)中抽取的5個(gè)含多次波的三維小面源CMP道集，圖13b為對(duì)應(yīng)的按標(biāo)量偏移距大小排列的三維CMP道集，可以清晰地看到多次波的存在（箭頭處）。

利用最小二乘拋物線(xiàn)Radon變換、基于L₁正則化的拋物線(xiàn)Radon變換以及文中提出的基于L_1/2正則化的拋物線(xiàn)Radon變換對(duì)以上數(shù)據(jù)進(jìn)行多次波壓制，結(jié)果如圖14所示。最小二乘拋物線(xiàn)Radon變換（圖14a）和基于L₁正則化的拋物線(xiàn)Radon變換（圖14b）的Radon數(shù)據(jù)中可以看到較多的假象，同相軸沒(méi)有很好地聚成一點(diǎn)；而本文方法的Radon域數(shù)據(jù)（圖14c）能量聚焦明顯，假象減少。

壓制多次波后的最終結(jié)果如圖15所示，可以清晰地看到圖15a、b中都未見(jiàn)多次波殘留，并且一次波同相軸未被影響，得到了較好的多次波壓制結(jié)果。

三維實(shí)際數(shù)據(jù)的多次波壓制問(wèn)題將會(huì)在后續(xù)的研究中繼續(xù)進(jìn)行。

5 結(jié)論與展望

本文系統(tǒng)闡述了拋物線(xiàn)Radon域壓制多次波的方法，提出了利用基于L_1/2正則化的拋物線(xiàn)Radon變換多次波壓制方法，得到以下結(jié)論：

1）由于L_1/2正則化可以獲得比L₁正則化更稀疏的解，所以基于L_1/2正則化的稀疏反演Radon變換具有更高的分辨率。

2）理論模型和實(shí)際數(shù)據(jù)結(jié)果證明本文方法具有較高的分辨率，重構(gòu)的地震數(shù)據(jù)與原始數(shù)據(jù)有更高的一致性，提高了多次波壓制的效果，減少了有用信號(hào)能量受到不必要的損失，為后續(xù)的一次波成像提供高質(zhì)量的數(shù)據(jù)基礎(chǔ)。

但在復(fù)雜地質(zhì)條件下，由于波場(chǎng)中微弱的有效信號(hào)容易被噪聲覆蓋，所以參數(shù)的設(shè)置需要更加謹(jǐn)慎，否則可能會(huì)因?yàn)閴褐圃肼暥鴵p失有效信號(hào)。另外，由于本文方法應(yīng)用L_1/2范數(shù)，雖然提高了地震信號(hào)的稀疏性，但使計(jì)算時(shí)間增加，所以此方法的計(jì)算效率相比于L₁范數(shù)方法沒(méi)有較大優(yōu)勢(shì)。而且由于提高了反演精度，會(huì)對(duì)地震數(shù)據(jù)的保幅性有較高的要求，所以針對(duì)以上局限性我們將會(huì)從以下三方面進(jìn)行改進(jìn)：

1）將考慮地震數(shù)據(jù)的振幅與偏移距的特征（amplitude versus offset， AVO），將本文方法與正交多項(xiàng)式變換相結(jié)合。

2）對(duì)求解算法進(jìn)行進(jìn)一步改進(jìn)，提高算法的計(jì)算效率。

3）將本文方法與其他多次波壓制方法相結(jié)合，對(duì)復(fù)雜地質(zhì)條件的地震資料進(jìn)行高效率高質(zhì)量的處理。另外可以與機(jī)器學(xué)習(xí)算法相結(jié)合，提高計(jì)算效率。

參考文獻(xiàn)（References）：

[1]Weglein A B. Multiple Attenuation： An Overview of Recent Advances and the Road Ahead[J]. The Leading Edge， 1999， 18 （1）： 40-44.

[2]李鵬，劉伊克，常旭，等.多次波問(wèn)題的研究進(jìn)展[J].地球物理學(xué)進(jìn)展，2006，21（3）： 888-897.

Li Peng， Liu Yike， Chang Xu， et al. Progress on the Multiple Problems[J]. Progress in Geophysics， 2006， 21（3）： 888-897.

[3]Robinson E A. Predictive Decomposition of Seismic Traces[J]. Geophysics， 1957， 22（4）： 767-778.

[4]Hampson D. Inverse Velocity Stacking for Multiple Elimination[J]. Canadian Journal of Exploration Geophysicists， 1986， 22（1）： 44-55.

[5]Ryu J V. Decomposition （DECOM） Approach Applied to Wave Field Analysis with Seismic Reflection Records[J]. Geophysics，1982， 47 （6）： 869-883.

[6]胡天躍， 王潤(rùn)秋，White R E. 地震資料處理中的聚束濾波方法[J]. 地球物理學(xué)報(bào)， 2000， 43（1）： 105-115.

Hu Tianyue， Wang Runqiu， White R E. Method in Seismic Data Processing[J]. Chinese Journal of Geophysics， 2000， 43（1）： 105-115.

[7]Weglein A B， Boyse W E， Anderson J E. Obtaining Three-Dimensional Velocity Information Directly from Reflection Seismic Data： An Inverse Scattering Formalism[J]. Geophysics， 1981， 46（8）： 1116-1120.

[8]Berryhill J R， Kim Y C. Deep-Water Peg Legs and Multiples：Emulation and Suppression[J]. Geophysics， 1986， 51（12）： 2177-2184.

[9]Berkhout A J， Verschuur D J. Estimation of Multiple Scattering by Iterative Inversion： Part I： Theoretical Considerations[J]. Geophysics， 1997， 62（5）： 1586-1595.

[10]王睿，王德利，胡斌，等.基于f-x EMD波場(chǎng)分離的繞射多次波壓制方法[J].吉林大學(xué)學(xué)報(bào)（地球科學(xué)版），2021，51（2）：597-606.

Wang Rui， Wang Deli， Hu Bin， et al. Diffracted Multiple Elimination Based on f-x EMD Wavefield Separation[J]. Journal of Jilin University （Earth Science Edition）， 2021， 51 （2）：597-606.

[11]Radon J. ber die Bestimmung von Funktionen durch ihre Integralwerte lngs gewisser Mannigfaltigkeiten， Ber[J].Verh Schs Akad， 1917， 69： 262-277.

[12]Sacchi M D， Ulrych T J. High-Resolution Velocity Gather and Offset Space Reconstruction[J]. Geophysics， 1995， 60（4）： 1169-1177.

[13]Herrmann P， Mojesky T， Magesan M， et al.De-aliased， High-Resolution Radon Transforms [C]//SEG 70th Annual International Meeting. Calgary： SGE， 2000： 1953-1956.

[14]熊登， 趙偉， 張劍鋒. 混合域高分辨率拋物Radon變換及在衰減多次波中的應(yīng)用[J]. 地球物理學(xué)報(bào)， 2009， 52（4）： 1068-1077.

Xiong Deng， Zhao Wei， Zhang Jianfeng.Hybrid-Domain High-Resolution Parabolic Radon Transform and Its Application to Demultiple[J]. Chinese Journal of Geophysics， 2009， 52（4）： 1068-1077.

[15]Li Z N， Li Z C.Multiple Attenuation Using λ-f Domain High-Resolution Radon Transform[J]. Applied Geophysics， 2013， 10（4）： 433-441.

[16]鞏向博， 韓立國(guó)， 王升超．混合域高分辨率雙曲Radon變換及其在多次波壓制中的應(yīng)用[J]. 石油物探， 2016， 55（5）： 711-718.

Gong Xiangbo， Han Liguo， Wang Shengchao. High-Resolution Hyperbolic Radon Transform and Its Application in Multiple Supression[J]. Geoghysical Prospecting for Petroleum，2016， 55（5）： 711-718.

[17]Chen S S， Donoho D L， Saunders M A. Atomic Decomposition by Basis Pursuit[J]. SIAM Review， 2001， 43（1）： 129-159.

[18]Chartrand R， Staneva V. Restricted Isometry Properties and Nonconvex Compressive Sensing[J]. Inverse Problems， 2008， 24（3）： 657-682.

[19]Xu Zongben， Zhang Hai， Wang Yao， et al.L_1/2Regularization[J]. Science China （Information Science）， 2010， 53： 1159-1169.

[20]Zeng Jinshan， Fang Jian， Xu Zongben.Spare SAR Imaging Based on L_1/2Regularization[J]. Science China （Information Science）， 2012， 55（8）： 1755-1775.

[21]Sun Le， Wu Zebin， Xiao Liang， et al. A Novel L_1/2Sparse Regression Method for Hyperspectral Unmixing[J]. International Journal of Remote Sensing， 2013， 34（20）： 6983-7001.

[22]康治梁，張雪冰.基于L_1/2正則化理論的地震稀疏反褶積[J].石油物探， 2019， 58（6）：855-863.

Kang Zhiliang， Zhang Xuebing. Seismic Sparse Deconvolution Based on L_1/2Regularization[J]. Geophysical Prospecting for Petroleum， 2019， 58（6）： 855-863.

[23]高棒棒.基于高分辨率稀疏Radon變換保幅處理方法研究[D]. 西安：西安石油大學(xué)， 2021：13-15.

Gao Bangbang. Research on Amplitude Preserving Processing Method Based on High-Resolution Sparse Radon Transform [D]. Xian： Xian Shiyou University， 2021： 13-15.

[24]張軍華，呂寧，雷凌，等. 拋物線(xiàn)拉冬變換消除多次波的應(yīng)用要素分析[J].石油地球物理勘探，2004，39（4）： 398-405.

Zhang Junhua， Lü Ning， Lei Ling， et al. Analysis of Applied Factors for Using Parabolic Radon Transform to Remove Multiple[J]. Oil Geophysical Prospecting， 2004， 39（4）： 398-405.

[25]Yang A Y， Sastry S S， Ganesh A， et al. Fast l₁-Minimization Algorithms and an Application in Robust Face Recognition： A Review[C]//2010 IEEE International Conference on Image Processin. Hong Kong： IEEE， 2010： 1849-1852.

[26]Daubechies I， Defrise M， De Mol C. AnIterative Thresholding Algorithm for Linear Inverse Problems with a Sparsity Constraint[J]. Communications on Pure and Applied Mathematics， 2004， 57（11）： 1413-1457.

[27]Beaton A E， Tukey J W. The Fitting of Power Series， Meaning Polynomials， Illustrated on Band-Spectroscopic Data[J].Technometrics， 1974， 16（2）： 147-185.

[28]Figueiredo M A T， Nowak R D， Wright S J. Gradient Projection for Sparse Reconstruction： Application to Compressed Sensing and Other Inverse Problems[J]. IEEE Journal of Selected Topics in Signal Processing， 2007， 1（4）： 586-597.

[29]Xu Zongben， Chang Xiangyu， Xu Fengmin， et al. L_1/2Regularization： A Thresholding Representation Theory and a Fast Solver[J]. IEEE Transactions on Neural Nwtworks and Learning Systenms， 2012， 23（7）： 1013-1027.

[30]Zuo Wangmeng， Meng Deyu， Zhang Lei， et al.A Generalized Iterated Shrinkage Algorithm for Non-Convex Sparse Coding[C]//2013 IEEE International Conference on Computer Vision. Sydney： IEEE， 2013： 217-224.

[31]李允怡.基于非凸正則化的圖像恢復(fù)技術(shù)研究[D]. 南京：南京郵電大學(xué)， 2020： 30-33.

Li Yunyi. Research on Image Recovery Technologyvia Nonconvex Regularization Method [D]. Nanjing： Nanjing University of Posta and Telecommunications， 2020： 30-33.

[32]Beck A， Teboulle M. A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems[J]. SIAM Journal on Imaging Sciences， 2009， 2（1）： 183-202.

[33]張巖，劉小秋，李杰，等. 基于時(shí)頻聯(lián)合深度學(xué)習(xí)的地震數(shù)據(jù)重建[J].吉林大學(xué)學(xué)報(bào)（地球科學(xué)版），2023， 53（1）：283-296.

Zhang Yan， Liu Xiaoqiu， Li Jie， et al. Seismic Data Reconstruction Based on Joint Time-Frequency Deep Learning[J]. Journal of Jilin University （Earth Science Edition）， 2023， 53（1）： 283-296.