信號自適應識別吸收補償方法

馬雄，李國發，劉立彬，桂志先

（1. 長江大學地球物理與石油資源學院，湖北武漢 430100；2. 油氣資源與勘探技術教育部重點實驗室（長江大學），湖北武漢 430100；3. 中國石油大學（北京）地球物理學院，北京 102249；4. 勝利油田分公司物探研究院，山東東營 257022）

0 引言

地下地層多為黏彈性介質，地震波在黏彈性介質中傳播會受到地層吸收作用的影響。該吸收效應是降低地震分辨率的主要原因，可見地層吸收補償是提高分辨率的根本途徑［1-3］。然而，在開展地層吸收補償的工程實踐中，不穩定性和對高頻噪聲的放大效應嚴重影響了該技術的應用及效果［4］。目前，國內外已推出多種地層吸收補償技術，如常用的非穩態反褶積［5］、反Q濾波［6-13］等。

反Q濾波技術的主要目的是補償地震記錄的能量損失，校正地震數據的相位畸變，提高地震資料的垂向分辨率。由于反Q濾波補償技術是一個能量指數放大過程，因此該技術本身存在嚴重的數值不穩定性問題［14］。對此，目前主要有兩種應對策略：一是對指數補償算子進行穩定化處理，如Wang［15］提出一種穩定化的反Q濾波算法，其核心是對振幅補償算子做穩定化處理，增強了反Q濾波的穩定性和抗噪性；二是將吸收補償問題進行數學轉化，利用反演框架實現反Q濾波，如Zhang等［16］在貝葉斯反演框架下，利用非穩態反演實現反Q濾波，較好地恢復了被吸收之前的地震信號。

以上反Q濾波方法的共性是基于單道吸收補償策略，其補償算法中沒有也無法考慮地震信號在空間上的約束關系。因此，反Q濾波結果的空間連續性較差，補償精度和穩定性都有所欠缺。

為了進一步提高地層吸收補償算法的可靠性，有人提出將地震信號空間結構特征作為正則化約束條件［17-18］，引入吸收補償算法中。地層的傾角和傾向、同相軸斜率、道間時差等是很簡單的空間結構特征，基于這些特征可發展基于空間反射結構正則化的地層吸收補償方法。例如，Ma 等［19］將傳統的單道吸收補償算法推廣到多道，并引入橫向約束項，實現基于橫向約束的地層吸收補償。該算法隱含的假設條件是地層（近似）水平，以保證橫向約束項（0°傾角約束）的有效性。在實際資料中，陡傾角也較常見，為此Ma 等［20］進一步發展了該算法，提出基于傾角約束的地層吸收補償方法。由于該改進算法需預先求取空間每個點的傾角信息，并將其作為傾角約束項引入多道吸收補償算法中，因此它具有較強的適應能力。

地震反射結構的傾向和傾角等幾何信息一般通過地震數據的空間梯度進行映射和估算，計算結果的抗噪性較差，且對地震數據品質要求較高。實際上，傾角和傾向等空間連續性信息隱含了地震反射的空間可預測性。因此，從理論上講，基于地震反射的空間可預測性，可構建一種對空間結構具有更強刻畫能力的空間反射結構表征算子，該算子既能描述地震反射幾何結構，也能描述反射振幅的空間變化。進一步地，將該空間反射結構表征算子作為一種先驗信息引入地層吸收補償的目標函數，賦予地層吸收補償反演系統一種信號自適應識別能力，使該吸收補償系統在一定程度上避免對噪聲干擾的放大效應，提高吸收補償算法的穩定性和抗噪性。

1 原理和方法

1.1 地震信號自適應識別的原理

空間可預測性是信號與隨機噪聲的本質差異。一般而言，地震信號在空間方向具有一定的連續性和相干性，這種相干性可通過空間反射結構算子進行預測和描述；隨機噪聲在空間分布上并不具備相關性，無法用空間預測算子表征。因此，兩者可通過空間反射結構算子進行甄別和區分。

假設有一個由20 個系數構成的地震信號空間預測算子bi，j，其形式［21］為

式中：bi，j中i、j分別指沿時間和空間方向的序號；“0”為被預測點位置；“×”表示不參與計算的點。

從bi，j的形式可知，被預測樣點值與其周圍樣點值的關系可表示為

式中：sk，l為地震數據樣點值；k和l分別為地震數據在時間（縱向）和空間（橫向）方向的索引號。

從式（2）可知地震信號空間預測過程為：將空間預測算子作為權系數，對被預測點周圍樣點值做加權求和，所得求和結果作為被預測點的值。這一數學運算隱含的物理意義是：地震信號的各樣點值并不是獨立存在的，它與周圍的樣點值應滿足預測濾波器所表征的數學關系；而隨機噪聲不具備特定的空間結構特征，故不能用上述空間預測算子進行數學表征。因此，利用上述空間預測算子可對地震信號與噪聲進行自適應識別和區分。

地震信號空間預測算子可通過求解如下最小二乘問題得到［21］

式中M和N分別表示地震數據在時間和空間方向的采樣點數。

1.2 基于反演的單道吸收補償方法

對干擾噪聲的放大效應是反Q濾波技術在實際應用中的主要弊端。其根本原因在于，雖然地震信號經歷了地層吸收衰減的影響，但相同時刻的噪聲干擾并未經歷與信號同樣的衰減作用，故反Q濾波技術在補償地震信號的同時，也會放大地震噪聲，造成補償結果的不穩定。傳統方法的解決思路是，在反演框架下，通過引入先驗信息或正則化約束條件，構建非穩態反演（吸收補償）的目標函數，從而在一定程度上緩解吸收補償結果的不穩定性。

黏彈性介質對地震波的吸收作用是一個典型的非穩態濾波過程，可由非穩態褶積模型描述［22］

式中：t和τ均為時間變量，t為數據記錄時間，τ為子波時間；d(t)表示包含吸收衰減的單道地震記錄；r(τ)為地震反射系數；w(t，τ)表示非穩態地震子波，其表達式為

式中：W(ω)為震源子波的頻譜，ω為角頻率；a(ω，τ)為吸收衰減函數（表征地下介質的吸收特性），其表達式為

式中：ωr為參考角頻率；Q為品質因子，且γ= 1/πQ。

將式（4）離散，可改成如下矩陣方程

簡記為

式中：d為單道地震記錄向量；W為時變子波矩陣；r為反射系數向量。

一般而言，可假設地震反射系數滿足稀疏概率分布，為此可采用L1范數對反射系數施加正則化約束，并構建如下的目標函數［14］

式中λ1為縱向調節因子，用于調節稀疏正則化約束項的權重。利用迭代重加權算法求解該目標函數，得到吸收補償結果為

式中：rn+1表示第n+1次迭代的反射系數結果；上標T 表示轉置運算；Ωn表示第n次迭代的輔助矩陣，且其中ε2為穩定化因子。

常規地層吸收補償方法是單道運行模式，即對地震數據中的各道依次進行非穩態反演。反演過程中未考慮、也無法考慮不同地震道補償結果之間的空間依賴關系。由此可見，傳統反Q濾波本身并不具備空間信號識別能力，補償過程中無法對有效信號與地震噪聲進行有效區分和選擇性補償。為了提高補償算法的穩定性和抗噪性，可將地震信號的空間可預測性作為一種先驗信息引入吸收補償的目標函數，賦予反演系統一種信號自適應識別能力，從而實現信號與噪聲的自適應識別和選擇性補償，在系統結構層面上消除噪聲干擾對吸收補償的影響。

1.3 信號自適應識別多道吸收補償方法

為了將空間預測算子描述的不同地震道之間的依賴關系引入吸收補償反演系統，首先需將單道非穩態褶積模型（式（8））推廣到多道形式，該多道非穩態褶積模型可表示為［22］

簡記為

上兩式中：dl、Wl、rl對應為第l（l=1，2，…，N）道的地震數據向量、時變子波矩陣及反射系數向量；s為多道地震記錄超級向量；m為多道反射系數超級向量；G為多道時變子波構建的塊狀矩陣。

在傳統吸收補償反演系統（式（9））中，只考慮了反射系數的縱向稀疏特性，并未考慮地震信號在空間方向的反射結構特征。為此，可利用空間預測算子bi，j對地震信號和隨機噪聲的自適應區分特性，將其作為一種先驗信息引入吸收補償的目標函數，構建一種信號自適應約束多道吸收補償反演系統，消除隨機噪聲對吸收補償的影響。

為實現反射結構自適應約束，首先對空間預測算子進行改造，得到反射結構表征算子

顯然，該算子與空間預測算子的唯一區別是中心位置（待預測點位置）的系數為-1。這意味著將該算子與待約束地震記錄（補償后記錄）進行褶積的結果不再是預測后的地震記錄，而是預測結果與實際測量結果的殘差，可表示為

式中：?為（二維）褶積運算符；bi，j此時為從補償前地震數據中提取的空間預測算子（或反射結構表征算子）；為補償后地震記錄；ek，l為預測誤差。若補償前、后地震數據反射結構特征類似，則預測誤差較小；反之，預測誤差較大。這也是利用反射結構表征算子約束補償后地震數據的核心思想。

式（14）的二維褶積運算較繁瑣，可采用Helix 變換將其轉化為一維褶積運算［22］

式中：e為預測誤差向量；B為反射結構表征算子褶積矩陣（利用空間預測算子bi，j構造得到），其作用是在反演過程中對地震信號和隨機噪聲進行自適應的識別和選擇性的補償；是由地震子波構成的分塊矩陣，其作用是將待反演的反射系數m轉化成地震記錄?，以便更好實現地震信號空間約束。

通過最小化預測誤差，就可使反演結果體現的反射特征與原始數據的反射特征最相似，即實現對吸收補償結果進行反射結構特征約束的目的。該反射結構約束項可表示為

將反射結構自適應約束項引入多道吸收補償反演系統，構建信號自適應識別多道吸收補償目標函數

式中：等式右端第一項為數據匹配誤差項（表征模型數據與實際數據的匹配殘差），第二項為縱向稀疏約束項，第三項為空間反射信號自適應約束項；λ2為信號空間約束項調節因子，其作用是調節信號自適應約束項的權重。

針對目標函數?(m)，Ma 等［22］采用交替方向乘子法（Alternating direction method of multipliers，ADMM）進行求解，該算法在理論上十分有效，已廣泛應用于壓縮感知和機器學習領域。但在算法測試時，發現ADMM 算法收斂速度較慢，往往需迭代幾十甚至上百次才能得到較可靠結果。而且，在算法實現過程中需引入中間變量和調節參數，這些參數的取值也會影響收斂速度。在某些情況下，參數選擇不當甚至可能帶來不收斂的后果。為此，本文采用迭代重加權算法［23］求解該目標函數，得到如下迭代表達式得到多道反射系數m后，可將其轉化為二維反射系數模擬[r1，r2，…，rM]；然后，將反演結果與子波褶積，得到補償后的地震記錄

綜上，本文提出的具體計算過程如下：

（1）利用式（3）從原始數據中計算地震信號的空間預測算子bi，j；

（2）基于吸收衰減理論，推導單道非穩態褶積模型（式（8）），并將單道褶積模型推廣到多道模型（式（12））；

（3）利用步驟（1）中求出的空間預測算子，構建空間反射結構約束項（式（16））；

（4）同時考慮反射系數的稀疏約束和地震數據的空間約束，建立信號自適應識別多道吸收補償目標函數（式（17））；

（5）求解式（17），得到反射系數反演結果（式（18））；

（6）根據式（18），計算補償后的地震記錄（式（19））。

2 模型測試與分析

首先采用層狀模型數據測試、評判本文方法的有效性和優越性。圖1a 是未考慮地層吸收效應的合成地震記錄，模擬采用的地震子波為50 Hz 的零相位雷克子波。可見未考慮吸收衰減的地震數據分辨率較高，地震反射結構清晰，深層能量與淺層相當。該合成記錄可作為參考數據（圖1a），用于定性和定量評價地層吸收補償結果的優劣。圖1b 是考慮了地層吸收衰減（Q=50）的合成記錄，可見由于地層的吸收衰減作用，接收的地震數據分辨率降低，深層能量較弱，原有的反射特征弱化。在該合成記錄中加入不同強度的隨機噪聲，得到信噪比（SNR）分別為20（圖1c）和5（圖1d）的含噪聲衰減記錄。

圖1 模型數據

分別采用單道吸收補償和信號自適應識別多道吸收補償對SNR=20 含噪地震記錄（圖1c）進行地層吸收補償，提高地震資料分辨率。圖2 展示兩種方法的補償結果，可見兩種方法的吸收補償結果都明顯增強了深層能量，也提高了原始數據的垂向分辨率。但本文方法可有效地抑制地震噪聲的放大，在提高分辨率的同時，有效地保持了地震記錄的信噪比和空間連續性，地震反射結構更清晰（圖2b）。

圖3展示兩種方法對SNR=5含噪地震記錄進行吸收補償的結果。隨著噪聲強度的增大，吸收補償的不穩定性問題更凸顯。可見此時單道吸收補償結果雖恢復了深層地震信號，提高了地震分辨率，但補償結果的信噪比較低。然而，得益于引入了地震信號空間結構自適應約束，本文方法補償結果（圖3b）地震反射結構更清晰，連續性更好。

圖3 SNR=5 含噪數據單道吸收補償（a）和自適應識別多道吸收補償（b）結果對比

進一步，定量計算了兩種方法吸收補償結果與參考數據的相關系數。在SNR=20 情況下，單道補償和本文方法補償結果與參考數據的相關系數分別為0.75 和0.91；而SNR=5 時，兩者相關系數分別為0.63 和0.82。定性分析和定量比較表明：本文方法具有更高的補償精度，更強的穩定性和抗噪性，其吸收補償結果空間連續性更好，與參考結果更匹配。

同時，還比較了采用ADMM算法和迭代重加權算法求解多道補償目標函數的效率。測試設備的CPU參數為13 thGen Intel（R） Core（TM） i7-13700KF，主頻為3.4 GHz。在SNR=20 情況下，ADMM 算法和迭代重加權算法的計算時間分別為172.55 和158.37 s；在SNR=5 時，ADMM 算法和迭代重加權算法的計算時間分別為185.04和161.21 s。顯然，本文算法運算效率更高，具一定優勢。

圖4顯示SNR=5時兩種方法補償結果與無衰減數據的單道對比。可以發現，相較于單道補償結果（藍色），本文方法補償結果（紅色）與參考曲線（綠色）無論是振幅還是相位都更吻合，補償效果更好。

圖4 SNR=5 時兩種方法補償結果與無衰減數據的單道對比

圖5是SNR=5時兩種方法補償結果與原始數據的單道振幅譜對比。可以看出，兩種方法均能有效補償地震高頻成分的衰減，拓寬地震數據有效頻帶。另外，單道補償振幅譜（藍色）比本文方法振幅譜（紅色）的高頻成分能量略強，這可能緣于單道算法對高頻噪聲的過度放大。

3 實際資料應用

利用實際地震數據，進一步對本文方法的實用性和可靠性進行了測試和分析。圖6a 為吸收補償前的實際地震數據，可見該輸入地震記錄分辨率較低，深層能量較弱，且視覺上看不到太多噪聲干擾。分別用單道吸收補償和信號自適應識別多道吸收補償方法處理該數據，得到相應的補償結果（圖6b 和圖6c）。其中單道吸收補償算法雖提高了原始數據的分辨率，但補償后地震剖面的噪聲較明顯，地震記錄信噪比較低。在補償過程中，通過引入地震信號自適應約束（本文方法），該算法具備了信號與噪聲的甄別能力。因此，在對地震信號進行補償的同時，并未放大噪聲干擾的影響，其補償結果具有更高信噪比和更好的空間連續性。圖7 為實際地震數據的局部放大顯示，可見經本文方法吸收補償后，地震剖面的整體結構更明確，空間連續性更好，地質現象和構造細節也更清晰。

圖6 原始實際資料和地層吸收補償結果

為了更充分展示本文方法對地震反射結構特征的恢復能力，將利用本文方法補償前、后的井旁地震數據與合成地震記錄（單道情形的噪聲放大過于明顯，故未被采用）進行標定對比（圖8）。可以看出，利用本文方法補償處理后，測井數據合成地震記錄與井旁地震記錄在反射結構特征上能較好地吻合，表明本文方法對實際資料也頗具實用性和可靠性。

圖8 本文方法（吸收）補償前（a）、后（b）地震數據與合成地震記錄的標定對比

圖9展示了吸收補償前、后地震數據的振幅譜，可見吸收補償后的地震數據高頻成分被抬升，地震頻帶有一定的展寬，說明兩種吸收補償方法都能有效提高地震資料的分辨率。但這兩種方法補償結果的振幅譜在50 Hz 以上存在一定差異，可能緣于單道補償方法對高頻噪聲的放大效應，這也從側面說明本文方法能有效壓制噪聲干擾，保障補償后地震資料的信噪比。

圖9 吸收補償前、后地震數據的振幅譜

4 結論

（1）本文提出一種信號自適應識別的吸收補償方法。該算法的核心是，在吸收補償目標函數中引入地震反射結構自適應識別約束項，在反演過程中實現信號與噪聲的自適應識別和選擇性補償。

（2）空間可預測性是有效信號與隨機噪聲的本質差異，利用空間預測算子可對地震信號與噪聲進行自適應識別和區分，并以此構建信號自適應約束正則化項。

（3）模型數據測試和實際資料應用結果表明：本文所提自適應識別吸收補償方法不僅能提高地震資料的分辨率，還可保護補償結果的空間構造特征和信噪比，具有實用性和較強穩定性。

本文方法為多道同時反演算法，計算量和存儲量較大。后期可根據矩陣特點進行優化處理；同時，也考慮通過并行計算進一步提升計算效率。