基于混合優化算法的地震數據匹配追蹤分解

蔡涵鵬 ，賀振華, ，高剛，黃德濟

(1. 中國石油川慶鉆探工程有限公司 地球物理勘探公司，四川 成都，610213；2. 成都理工大學 油氣藏地質及開發工程國家重點實驗室，四川 成都，610059；3. 成都理工大學 地球物理學院，四川 成都，610059；4. 成都理工大學 地球探測與信息技術教育部重點實驗室，四川 成都，610059)

匹配追蹤算法是Mallat等[1]提出作為信號稀疏分解中應用較為普遍的一種算法。匹配追蹤算法克服了時窗傅里葉變化不能描述信號結構變化和小波變換不能提供精確估算子波原子頻率的缺點[2]。匹配追蹤算法具有較高的時頻分辨率、暫態結構的局部自適應性、信號結構的參數表示等優良特性，能夠較好地適應信號分析[3-6]。最近，國內外一些學者將匹配追蹤算法應用于地震信號分析[7-10]。匹配追蹤算法從一個過完備的子波原子庫出發，采用某種策略每次選取與地震道信號最匹配的基(也稱為原子域或時頻原子)，最終將地震道信號表示為若干時頻原子的線性組合。應用匹配追蹤算法，創建超完備子波庫所使用的子波原子基與原始信號越相似，則分解效果越好，因此，需要選擇合適的原子基。本文作者以Morlet小波作為基本匹配子波，其原因在于Morlet小波適合具有能量吸收和速度頻散的地震信號的時頻分析[11-12]。由于地震信號具有非平穩特征，為了能夠對地震道信號更好地進行時頻譜分析，需要選擇具有幅度、時間延遲、尺度、頻率和相位移參數的小波原子。常規的匹配追算法是一種不斷迭代、尋求最佳匹配的貪婪算法，子波庫中的可供選擇用于匹配的子波是固定的，1次迭代只能確定1個匹配子波，對子波庫的全方位掃描會需要大量的計算機運行時間，導致搜索效率很低，不適合于大規模三維地震資料處理。為了提高計算效率，Liu等[9-10]將地震信號的瞬時特征引入到匹配追蹤算法中，Wang[13]提出三步法實現匹配追蹤，張繁昌等[14]提出了雙參數快速匹配追蹤算法，張顯文等[15]提出兩步法匹配追蹤算法，提高了匹配追蹤算法的效率，并用于河道砂體識別、隨機噪聲消除、含氣儲層烴類檢測和時頻分析等。每次迭代尋找最優匹配小波的過程本質上是一個尋找非凸的，高度非線性問題的全局最優解的過程。本文作者以5個參數控制的Morlet小波作為基本匹配子波，提出全局優化算法與局部優化算法相結合的混合優化算法，實現一次性計算描述Morlet小波的5個參數，使得計算效率和精確度得到提高。匹配追蹤算法具有較強的多解性。為了減少多解性和使得分析信號更快地被分解，利用殘差信號能量與分解子波能量之比最小作為目標函數。對基于混合優化算法快速匹配追蹤算法的迭代次數進行有效控制，不僅有效地避免了毫無意義的計算，在一定程度上也能減少解的多解性。由于匹配追蹤算法具有自適應的分解特點，具有較高的時頻分辨率，考慮到各個子波原子之間的獨立性，通過計算每次分解得到子波原子的Wigner-Ville 分布，將所有原子的時頻分布疊加到一起構成整個地震道信號的高分辨時頻分布，以滿足高精度時頻分析、流體檢測和時頻屬性提取的要求，為烴類檢測和儲層描述等提供有效的手段。

1 方法原理

1.1 匹配追蹤基本思想

構建信號的時頻原子基函數集，它由窗函數經平移、伸縮、相位移和頻率調制而成，其形式如下：

式中：σ為尺度因子；τ為時間平移因子；fm為頻率調制因子；φ為相位移因子；記 γ=γ(σ，τ，fm，φ)，γ∈Γ；Γ=R+×R×R；g(t)為窗函數。

利用式(1)構建1個時頻原子基函數集合M{gγ(t)}γ∈Γ，匹配追蹤算法通過迭代的方式，每次從時頻原子集合中選擇一個時頻原子gγi(t)∈M，將需要分解的信號 s(t)投影到某個時頻原子 gγi(t)上，并計算信號的剩余部分Rs，使得

式中：An是第n次提取的時頻原子基gγn(t)的幅度；RsN是原始信號經過N次迭代分解后的殘余信號。

匹配追蹤算法成功的關鍵之一是選取適合于信號分析的時頻原子。對于地震信號，選擇Morlet小波作為時頻原子，其原因在于Morlet小波與地震子波具有良好的相似性，可以通過Morlet小波的伸縮、頻移、相位移和頻率調制來匹配實際地震數據，得到較好的分解結果。由5個參數表征的Morlet小波的時間表達式為：

式中：A為小波的幅度。

匹配追蹤算法成功的關鍵之二是精確地獲取用于描述時頻原子的參數。

1.2 參數優化算法

匹配追蹤算法是一種不斷迭代尋找最佳匹配的貪婪算法[1]，從巨大的時頻原子集合形成的解空間中反復迭代搜索最佳匹配的計算效率低，是一個NP完全問題[16-17]，為了保證在實際地震數據處理中的實用性和應用性，需要提高最優時頻原子的搜索運算效率。Liu等[9-10,13-14]描述時頻原子的參數獲取方法均是基于復數道分析的思想，應用復數道分析獲得時頻原子的時間位移、相位移、頻率等參數。然而，本文作者在研究中發現，針對于一個單波而言，復數道分析獲得瞬時振幅包絡最大值以及對應位置的瞬時相位和瞬時頻率能夠準確地確定時頻原子的時間位移、相位移和頻率，而對于像地震信號的調諧復波，復數道分析獲得參數并不能精確地用于描述時頻原子的參數，尤其是頻率屬性(瞬時頻率可能為負值)。匹配追蹤算法的整個分解過程實質上是一個多參數尋優過程，而粒子群優化算法是一種簡單、有效的隨機全局優化算法，可以解決大量非線性和多極值的復雜問題，具有概念簡單，易于實現，搜索速度快的優點，但是，由于粒子群優化算法中所有粒子都向最優解的方向飛去，粒子趨向同一化，群體多樣性逐漸消失，致使后期算法的收斂速度明顯變慢，甚至處于停滯狀態為了克服這一缺點，引進在局部算法中有著完善的數學理論基礎，采用不精確線性搜索時的超線性收斂性和處理實際問題能力的BFGS算法。本文將粒子群優化算法與BFGS算法相結合的混合優化算法引進到匹配追蹤算法中，解決匹配追蹤過程中的參數尋優問題，提高搜索最優解的精確度和效率。

1.2.1 粒子群優化算法

粒子群優化算法是一種基于迭代尋優的優化算法[18-19]。鳥被作為沒有質量和體積的微粒，并延伸到N維空間。粒子在 N維空間的位置表示為 Xi=(x1，x2，…，xN)的向量，飛行的速度表示為Vi=(v1，v2，…，vN)的向量。每一個粒子都有一個目標函數的適應度，并且已知粒子本身到目前為止發現的最好位置 Pbest(粒子自己的飛行經驗)、現在的位置Xi和到目前為止整個群體中所有粒子發現的最佳位置Pbest,g(粒子同伴的飛行經驗)。粒子就是通過自己的經驗和同伴中最好的經驗來確定下一次飛行。粒子群優化算法首先將粒子初始化為一群隨機粒子，然后，通過下列公式更新自己的速度位置，迭代尋找最優解：

式中：Vin為粒子i在第n次迭代的速度；ω為慣性因子；c1和c2為學習因子，為正常數，分別代表粒子的局部和全局學習能力；r1和 r2為[0，1]的隨機數。另外，為了粒子在解空間內飛行，在迭代過程中，設置飛行速度和粒子位置的邊界。但是，由于在搜索過程中，粒子的搜索方位局限在 Pbest和 Pbest,g附近，因此不能保證以概率1收斂到全局最優解。

1.2.2 BFGS方法

BFGS方法是一個擬牛頓方法，具有二次終止性、整體收斂性和超線性收斂性，算法產生共軛的搜索方法[20]。算法的收斂性速度比僅為線性收斂的梯度最速下降法的速度快。在實際計算過程中，由于舍入誤差的存在及一維搜索的不精確性，導致BFGS算法效率的影響也比梯度最速下降法和共軛梯度法小，且精度均較后兩者的高。BFGS方法是一個有效的局部優化算法。BFGS算法求解無約束優化問題 min(f(x))的主要步驟如下：

第1步：確定變量維數N和BFGS算法收斂精度ε，給初始解賦值x0∈Rn，k=0；

第3步：沿著方向dk做線性搜索求取αk＞0，且令 xk+1=xk+αkdk；

第4步：校正Hk產生Hk+1使得擬牛頓條件成立。其中：

第5步：令k=k+1，轉入第2步。

在上述的算法中，初始Hesse逆近似H0通常取單位矩陣。

1.2.3 混合算法描述

為了克服具有良好全局優化的粒子群優化算法在優化后期階段收斂速度慢，甚至停滯狀態的問題，在粒子群優化算法后期中引入BFGS方法，利用BFGS方法的整體收斂性和超線性收斂性來加速收斂速度。BFGS方法屬于局部優化算法，其優化結果的好壞取決于初始位置的選擇，因此，利用具有全局搜索能力的粒子群算法提供給BFGS方法初始解位置。將兩種優化算法結合的混合算法就能取長補短，改善優化的結果和速度。簡單的混合使用2種優化算法，并不能取得理想的效果，原因在于粒子群優化算法可能過早陷入局部最優，使得粒子群優化算法變異能力差。由于BFGS方法中采用不精確一維搜索，在搜索過程中，設定最大搜索次數，若超過最大搜索次數，則表明BFGS方法已經陷入局部最優，同時也表明上次粒子群優化算法結果并不是BFGS算法的一個良好初值，此時，應該提前結束BFGS搜索，在BFGS算法所得的最后一次解的較大范圍內生成新的種群空間。此外，若BFGS算法搜索過程中gk+1≤ε，則比較目標函數與預設值之間的大小，若目標函數值大于預設值，則在當前解x*較小的空間范圍內生成新的種群空間xk，然后，從新的種群空間中開始進行粒子群搜索，重復此過程直到目標函數小于預設值為止。生成新種群空間的公式為

式中：c為調節常數；rand( )∈(0，1)的隨機數。

1.3 目標函數及迭代終止準則

上述的優化算法均是求取目標函數極小值。優化的參數包括描述時頻原子的幅度、頻率、相位移、尺度因子和時間位移5個變量。為了表達方便，將式(2)的目標函數改寫為

即使得每次分解后剩余信號能量與分解子波能量比值最小。

由于實際的地震記錄中包含隨機噪聲，在迭代初期，殘差能量急劇減少，但是在一定水平時，殘差能量基本保持一個相對穩定的水平。其原因在于當分解至噪聲水平時，隨機噪聲使得從殘差信號中減去匹配子波時出現“此消彼長”的顯現，使得剩余信號能量變化非常小，將此時剩余信號能量作為收斂能量閥值。一旦信號分解后剩余能量達到收斂能量閥值，繼續迭代將耗費機時，且毫無意義。因此，在實際地震信號匹配追蹤分解時，可以隨機抽取若干地震道進行試驗分析，確定收斂能量閥值以及分解次數，可以提高匹配追蹤的效率。在實際的程序中也可以設置分解后剩余能量相對于分解前信號能量的百分比進行終止分解條件，從而達到能夠自適應地終止匹配追蹤分解。

1.4 構建時頻譜

信號 s(t)分解成為一組時頻原子的線性組合后，再求出每個時頻原子的 Wigner-Ville分布，既可以保持 Wigner-Ville分布具有高的時頻分辨率特點，又可以避免直接求原信號Wigner-Ville分布的交叉項干擾，能夠得到高質量的時頻能量分布。針對Morlet小波作為時頻原子基，時頻振幅譜的解析表達式為：

式中：fm,n為第n個匹配子波的平均頻率；||gγn||為標準化因子，其解析表達式為

利用式(10)能夠提高式(8)和(9)的計算速度。

2 數值試驗分析

圖1所示為由不同主頻、相位、時間位移、尺度因子及幅度的4個子波合成的地震記錄。利用本文提出的方法將合成地震記錄進行分解，并利用式(3)進行重構。經過自適應快速時頻分解后，得到的時頻譜分布如圖 2(a)所示，由于在時頻原子基的搜索過程中采用了具有全局尋優的微粒群優化算法及局部尋優的BFGS算法，使得搜索速度加快，運算效率得到一定程度提高。比較圖2(a)與基于短時傅里葉變換(圖2(b))的時頻譜可見：基于匹配追蹤分解的各個信號分量在時頻平面的位置精確，具有良好的時頻聚集性，它們的時頻分辨率幾乎達到極限。圖2(b)所示的時頻聚集性能明顯比前者的差。短時傅里葉變換聚集性能差的原因在于強烈依附于窗函數的選擇。圖3所示為合成地震記錄與匹配追蹤重構地震記錄的對比，由圖3可見：地震道可以完全由匹配子波重構，其重構地震記錄與原始合成地震記錄的殘差與原始信號在同一數量級上比較，其強度和能量很弱，說明分解是可靠的。其合成地震記錄的分解次數和能量衰減關系如圖4所示。由圖4可見：本文提出的基于混合優化算法的匹配追蹤方法能夠精確地獲得每次分解子波原子的控制參數，且能量的衰減趨勢表明分解過程是收斂的。

圖1 不同主頻、相位、尺度、時移、幅度的Morlet子波合成的地震記錄Fig.1 Synthesis signal records using Morlet wavelets with different dominant frequencies, phases, scale factors, time delays, and magnitudes

圖2 合成信號的時頻譜Fig.2 Time-frequency spectral for synthesis signal

圖3 合成地震記錄的匹配追蹤結果Fig.3 Results of matching pursuits for synthetic seismic records

為了驗證算法的穩定性，在合成地震記錄中(圖1)加入最大幅度為 0.2的高斯隨機噪聲。利用本文提出的匹配追蹤算法將其進行分解，并用各個匹配子波進行信號重構。圖5所示為加入隨機噪聲后的匹配追蹤結果。由圖5可見：匹配追蹤結果較好地將有效信號進行重構，并能夠有效消除隨機噪聲。值得注意的是，對含有噪聲的信號，首先分析信號的信噪比，也能夠有效地確定匹配追蹤分解的終止準則。

圖4 分解次數與能量衰減的關系Fig.4 Relationship between decomposition times and energy attenuation

圖5 合成地震記錄加隨機噪聲后的匹配追蹤結果Fig.5 Results of matching pursuits for synthetic seismic records added random noise

在控制morlet時頻原子的5個參數中，尺度因子σ有著重要的意義。尺度因子σ控制著子波的延續長度。σ越大，子波延伸長度越長，反之，子波長度延伸越短。在地震波的有效信號中，子波的長度是有一定限度的，因此，在應用匹配追蹤算法重構地震信號時，可以通過控制尺度因子對含噪地震信號進行噪聲壓制。圖6所示為含噪合成地震信號應用匹配追蹤分解后，消除尺度因子σ大于4的子波以及噪聲的效果圖。從匹配追蹤結果中消除σ較小的子波，能夠有效地壓制具有尖峰特征的噪聲。通過從匹配追蹤結果中消除σ特別大的子波，能夠有效壓制像正弦或者余弦特征的噪聲，例如，地震測線上空高壓輸電線電流引起的地震記錄中從淺到深層的1個50 Hz左右的單頻干擾。

圖6 尺度因子濾波結果Fig.6 Filtering results of scale factor

3 實例分析

基于混合優化算法完成的匹配追蹤方法應用于實際地震數據的噪聲壓制和時頻譜分析。圖 7(a)所示為原始地震振幅剖面(圖中從上到下，第 1個層位為E3d2U，第2個層位為E3d2M，第3個層位為E3d2L，豎線為鉆井位置，鉆井資料顯示氣藏位于白色箭頭標注，即E3d2M層之上)。原始地震數據中包含較大的隨機噪聲，導致資料的信噪比較低。從原始地震數據中減去匹配追蹤獲得的尺度因子σ＜0.3和σ＞8的子波分量以及匹配追蹤過程終止時剩余信號，得到重構的地震振幅剖面(圖 7(b))。經過匹配追蹤濾波重構的地震振幅剖面中信噪比明顯得到提高，同相軸的連續性得到加強，斷點更加清晰。由于基于控制參數σ和剩余能量信息的濾波與常規的截頻濾波(如F-K濾波)有顯著的差異，本文方法能夠有效地保留地震數據的有效頻率信息。

圖7 地震振幅剖面Fig.7 Seismic profile

圖8 基于匹配追蹤的共頻率剖面Fig.8 Isofrequency profile using matching pursuits

由于大地濾波作用和地層中流體的黏滯性導致地震波的高頻吸收快，地震子波逐漸變長。Morlet小波適當地描述了孔隙介質中波傳播時能量吸收和相位失真特征，因此，Morlet小波適合于含氣儲層的聲學特性分析。Castagna等[8]指出，在含氣儲集層的下方存在低頻陰影的效應，表現為在低頻切片上具有較高的能量異常。圖8所示為一組基于本文方法獲得的不同頻率的共頻率剖面。從圖8可見：在12 Hz(圖8(a))時，油氣層(白色箭頭標注)顯示為強能量，其下部瞬時譜能量也強(黑色圈注)，即存在低頻“上強下強”的低頻伴影特征。隨著頻率增至20 Hz和28 Hz(圖8(b)和(c))，油氣層為強能量顯示，但下部伴影能量明顯減弱，而當頻率增至36 Hz時(圖8(d))，油氣層仍顯示為強能量，但下部的伴影能量消失。

4 結論

(1) 所提出的基于具有全局優化能力的粒子群優化算法和具有局部優化能力的BFGS算法的自適應匹配追蹤算法，使得匹配追蹤算法不再依賴于復數道分析確定時頻原子的振幅、頻率和相位的初值，能夠在整個解空間中搜索最優的時頻原子控制參數，提高了匹配追蹤的精確度，并且使得匹配追蹤算法的計算效率提高。應用局部函數的解析表達式和根據信噪比確定匹配追蹤的停止準則，能夠進一步提高計算效率。

(2) 通過選擇不同的子波尺度因子和剩余能量占原始信號能量的百分比，能夠有效地壓制噪聲，并且該算法結合 Wigner-Ville分布可以對地震信號進行快速時頻譜分析，刻畫不同頻率的響應特征，為烴類檢測和儲層描述提供有效的手段。

[1] Mallat S, Zhang Z. Matching pursuit with time- frequency dictionaries[J]. IEEE Trans Signal Process, 1993, 41(12):3397-3421.

[2] Qian S, Chen D. Signal representation via adaptive normalized Gaussian functions[J]. Signal Processing, 1994, 36(1): 1-11.

[3] Rebollo-Neira L, Lowe D. Optimized orthogonal matching pursuit approach[J]. IEEE Signal Process Letters, 2002, 9(4):137-140.

[4] Capobianco E. Independent multi-resolution component analysis and matching pursuit[J]. Computational Statistics & Data Analysis, 2003, 42(6): 385-402.

[5] Andrle M, Rebollo-Neira L, Sagianos E. Backward-optimized orthogonal matching pursuit approach[J]. IEEE Signal Process Letters, 2004, 11(9): 705-708.

[6] Andrle M, Rebollo-Neira L. A swapping-based refinement of orthogonal matching pursuit strategies[J]. Signal Processing,2006, 86(3): 480-495.

[7] Wang B, Pann K. Kirchhoff migration of seismic data compressed by matching pursuit decomposition[C]//66th Annual International Meeting. Denver, Colorado: SEG, Expanded Abstracts, 1996: 1642-1645.

[8] Castagna J P, Sun S J, Siegfried R W. Instantaneous spectral analysis: Detection of low-frequency shadows associated with hydrocarbons[J]. The Leading Edge, 2003, 22(2): 120-127.

[9] Liu J, Marfurt K J. Matching pursuit decomposition using Morlet wavelet:[C]//75th Annual International Meeting. Houston, Texas:SEG, Expanded Abstracts, 2005: 786-789.

[10] Liu J, Wu Y, Han D, et al. Time-frequency decomposition based on Ricker wavelet[C]//74th Annual International Meeting.Denver, Colorado: SEG, Expanded Abstracts, 2004: 1937-1940.

[11] Morlet J, Arens G, Fourgeau E, et al. Wave propagation and sampling theory: Part I. Complex signal and scattering in multilayered media[J]. Geophysics, 1982, 47(2): 203-221.

[12] Morlet J, Arens G, Fourgeau E, et al. Wave propagation and sampling theory: Part II. Sampling theory and complex waves[J].Geophysics, 1982, 47(2): 222-236.

[13] Wang Y H. Seismic time-frequency spectral decomposition by matching pursuit[J]. Geophysics, 2007, 72(1): v13-v20.

[14] 張繁昌, 李傳輝, 印興耀. 基于動態匹配子波庫的地震數據快速匹配追蹤[J]. 石油地球物理勘探, 2010, 45(5): 667-673.ZHANG Fanchang, LI Chuanhui, YIN Xingyao. Seismic data fast matching pursuit based on dynamic matching wavelet library[J]. Oil Geophysical Prospecting, 2010, 45(5): 667-673.

[15] 張顯文, 韓立國, 王宇, 等. 地震信號譜分解匹配追蹤快速算法及應用[J]. 石油物探, 2010, 49(1): 1-6.ZHANG Xianwen, HAN Li-guo, WANG Yu, et al. Seismic spectral decomposition fast matching pursuit algorithm and its application[J]. Geophysical Prospecting for Petroleum, 2010,49(1): 1-6.

[16] Gribonval R. Fast matching pursuit with a multi-scale dictionary of Gaussian Chirps[J]. IEEE Trans Signal Processing, 2001,49(5): 994-1001.

[17] Davis G, Mallat S, Avellaneda M. Adaptive greedy approximations[J]. Constr Approx, 1997, 13(1): 57-98.

[18] Kennedy J, Eberhart R. Particle Swarm optimization[C]//Proceedings of IEEE International Conference on Neural Network. Perth, Western Australia, 1995: 1942-1948.

[19] Sun J, Xu W B. Particle swarm optimization with particles having quantum behavior[C]//IEEE Proceedings of Congress on Evolutionary Computation. Piscataway, New Jersey, 2004:325-331.

[20] 袁亞湘, 孫文瑜. 最優化理論與方法[M]. 北京: 科學出版社,2003: 219-234.YUAN Yaxiang, SUN Wenyu. Optimization theory and methods[M]. Beijing: Science Press, 2003: 219-234.