基于OpenMP并行計算的匹配追蹤時頻分析方法

鄧世廣 王淑艷 趙文津 劉志偉 何 潤(①中國地震臺網中心，北京 100045； 中國地質科學院，北京 100037； 中電科海洋信息技術研究院有限公司，北京 100041； 中國石油勘探開發研究院西北分院，甘肅蘭州 730020)

1 引言

匹配追蹤算法由Mallat等[1]率先提出，其本質是對信號做稀疏分解。與其他時頻方法相比，匹配追蹤方法具有更高的時頻分辨率[2]，目前在油氣檢測、薄層識別以及噪聲壓制等方面有著廣泛應用[3-8]。鑒于傳統的匹配追蹤貪婪迭代算法計算量過大，Liu等[9，10]先后提出了基于Ricker子波和Morlet子波庫的匹配追蹤算法，通過Hilbert變換求取地震信號的瞬時屬性(振幅、相位、頻率)，再引入匹配追蹤算法，顯著提高了計算效率。張繁昌等[11，12]分別運用雙參數快速匹配追蹤算法、動態最優搜索和時頻原子正交變換方法，均不同程度地提高了計算效率。邵君[13]對基于匹配追蹤的信號稀疏分解算法進行了研究。基于Liu等[9，10]的成果，Wang[14]通過引入控制Morlet小波的尺度參數，提出了三步法匹配追蹤算法，提高了算法的效率。范興利等[2]探討了以Morlet小波作為時頻原子的匹配追蹤算法中尺度因子對信號與時頻原子匹配特征的控制作用。

隨著計算機硬件的快速發展，配備多核處理器的計算機已十分普遍，而以往的匹配追蹤算法基本都是串行執行，未能發揮多核處理器的計算優勢。OpenMP(Open Multi-Processing)采用共享內存的編程模式，方便各個線程對數據的讀取與調用，可充分發揮多核處理器的性能[15]。通過可并行性分析，本文認為匹配追蹤時頻分析算法適合并行運算，并基于OpenMP實現了以Morlet小波作為時頻原子的匹配追蹤算法的并行運算，另外對計算多個Morlet子波的平滑偽Wigner-Ville分布(SPWVD)的過程也采用并行運算。實際數據的計算結果表明，基于OpenMP并行計算的匹配追蹤時頻分析方法在保持計算精度的前提下能充分發揮多核處理器的性能，有效提高計算效率。

2 匹配追蹤時頻分析原理

2.1 匹配追蹤算法思想

匹配追蹤算法是將信號投影到一系列時頻原子之上，通過這些時頻原子的線性組合精確表示復雜原始信號。這里所謂的時頻原子在實現過程中可具體化為由頻率、相位、時間延遲、尺度等參數確定的小波，在地震信號處理中Ricker 小波和Morlet 小波是兩種常用時頻原子形式[2]。

匹配追蹤算法思想以數學的形式可表示為：設H表示希爾伯特空間，定義H中的一個時頻原子庫D，gγ∈D為時頻原子庫中的一個時頻原子，原子庫中的每個原子都應該是歸一化的，即滿足‖gγ‖=1。之后，從原子庫中選擇與信號f(f∈H)最佳匹配的原子gγ 0，這時信號f就可表示為

(1)

(2)

(3)

按此過程不斷迭代分解，直到殘差能量小于所設閾值，這樣信號f最終被分解為

(4)

將信號和時頻原子抽象表示為高維空間中的矢量，則匹配追蹤的算法思想可由圖1表示[2，16]。

圖1 匹配追蹤示意圖

2.2 基于Morlet子波的匹配追蹤算法

Morlet等[17]提出的Morlet子波在時間域的表達式為

exp{i[ωm(t-u)+φ]}

(5)

式中：ωm是子波的中心(角)頻率；u是時間延遲；σ是尺度參數；φ是相位參數。通過這4個參數可以確定唯一的Morlet子波。基于Morlet子波的匹配追蹤算法是一個迭代計算的過程，每次迭代都會估算出一個最佳匹配的子波mn，其中n是迭代次數。經過N次迭代，一個地震信號f(t)可表示為N個Morlet子波的線性組合

(6)

第一步，首先通過復地震道分析技術估算Morlet子波的三個參數(un,ωn,φn)。對于一個實際的地震道，可通過Hilbert變換獲得其復地震道，并以此進一步計算得到復地震道屬性[18]。以復地震道最大包絡處的時間為子波的時間延遲參數un，瞬時頻率作為中心頻率ωn，瞬時相位作為相位參數φn。然后利用下式搜尋尺度參數σn

(7)

第二步，利用式(7)進一步優選四個參數。圍繞某一參數ξ的搜索范圍是[ξ-Δξ,ξ+Δξ]，其中Δu是時間采樣間隔、Δσ是尺度參數σ的采樣間隔、Δω是頻率采樣間隔、Δφ為5°。值得注意的是該參數優選的過程需大量的運算。

第三步，通過下式計算最佳匹配子波mγn的振幅

(8)

(9)

經過N次迭代之后的殘差R(N)可認為是數據中的噪聲。

2.3 Morlet子波的平滑偽Wigner-Ville分布

通過短時傅里葉變換得到地震信號的時頻分布是目前常用的時頻分析方法，但通過短時傅里葉變換方法分析信號時頻特性時無法同時滿足較高的時間分辨率與頻率分辨率[19]。Wigner-Ville分布無需選取窗函數，避免了線性時頻表示中時間分辨率與頻率分辨率的相互矛盾，因而具有更好的時頻分辨率。對于地震信號s(t)的Wigner-Ville分布可以寫成

(10)

式中：t是時移[20]。但當信號s(t)=s1(t)+s2(t)時，根據式(10)，其Wigner-Ville分布可寫為

(11)

可進一步簡寫為

Ws(t,ω)=Ws1(t,ω)+Ws2(t,ω)+

2Re[Ws1,s2(t,ω)]

(12)

式中： Re[Ws1,s2(t,ω)]為交叉項；Ws1(t,ω)和Ws2(t,ω)分別是信號s1(t)和s2(t)的Wigner-Ville分布，稱為自相關項。多分量復雜信號的Wigner-Ville分布中交叉項的存在會嚴重影響信號的時頻分布(圖2)。為了抑制交叉項干擾，人們提出平滑偽Wigner-Ville分布(SPWVD)[21]，即對WVD分別增加時域和頻域窗函數g(u)和h(τ)，其表達式為

式中g(u)和h(τ)是兩個實對稱窗函數，且h(0)=g(0)=1。可通過選擇窗函數g(u)和h(τ)來控制時域和頻域的平滑尺度。

圖2 Wigner-Ville 分布的交叉項干擾

基于匹配追蹤的時頻分析的思路及流程(圖3)是： ①通過匹配追蹤將復雜信號分解為多個Morlet子波； ②計算各個子波的平滑偽Wigner-Ville分布[22]； ③疊加各子波的時頻分布得到原信號的時頻分布。這樣就能保證多分量復雜信號s(t)的時頻譜在具有較高時頻分辨率的同時又有效避免了交叉項的干擾。 圖4顯示了對地震信號(圖4a)進行匹配追蹤運算得到了一系列的Morlet子波(圖4b)。如果直接計算該地震信號的平滑偽Wigner-Ville分布，其結果受到了嚴重的交叉項干擾，已經難以識別信號在不同時刻對應的頻率分布(圖4c)。采用基于匹配追蹤的時頻分析方法，分別計算分解得到的Morlet子波的平滑偽Wigner-Ville分布，然后疊合得到整個信號的時頻分布，這樣有效避免了交叉項的干擾，信號對應的時頻特征十分清晰(圖4d)。

圖3 基于匹配追蹤的時頻分析方法的計算流程

圖4 基于匹配追蹤的平滑偽Wigner-Ville分布計算結果(a)地震信號； (b)匹配追蹤分解得到的Morlet子波； (c)直接計算的平滑偽Wigner-Ville分布； (d)基于匹配追蹤計算的Wigner-Ville分布

3 基于OpenMP的并行實現與分析

隨著計算機硬件的迅速發展，配置多核處理器的計算機已越來越普遍。在多核處理器的硬件基礎上，以共享物理內存的形式的并行運算得到了快速發展。OpenMP以一種面向共享存儲的并行程序標準于1997年正式發表，之后又陸續發表了面向Fortran、C/C++的四個標準。經過幾年的推廣，OpenMP已成為關于共享存儲并行計算的新的工業標準，目前在中小型的共享存儲并行系統上得到了廣泛應用[23]。

Fork-Join類型的OpenMP程序機構是最基本的OpenMP程序方式，也是應用最廣泛的一種OpenMP程序形式[23]。Fork-Join(分叉—連接)模式的含義是：當程序開始時，由主線程執行程序的串行部分，在接到并行任務時，主線程執行Fork操作，派生出其他線程共同完成并行部分，完成后各線程執行Join操作，回到主線程執行串行部分[23-26](圖5)。通過這種Fork-Join模式，編譯器可以充分利用多核處理器的優勢完成并行計算任務。

匹配追蹤算法需要龐大的運算量，Wang[14]提出的三步法有效提高了匹配追蹤的計算效率，但所需計算量依然較大。為提高匹配追蹤算法的計算效率，首先對該算法進行了可并行性分析，之后利用OpenMP實現了基于匹配追蹤時頻分析方法的并行運算，從而在多核處理器環境下有效提高其運算效率。

圖5 Fork-join模式示意圖

基于匹配追蹤的時頻分析算法中大量的運算主要集中在搜尋最佳匹配子波的過程以及計算所有子波平滑偽Wigner-Ville分布的過程(圖3中灰框所示的步驟)。因此分析這兩個步驟是否適合并行運算。

3.1 子波尋優的并行計算

匹配追蹤搜尋最佳匹配子波的過程本質上是一個多重嵌套循環過程，其示意代碼如下

for (時間延遲參數)

{

for (頻率參數尋優范圍)

{

for (相位參數尋優范圍)

{

for (尺度參數尋優范圍)

{

尋優運算

}

}

}

}

這里的尋優運算實際上就是計算Morlet子波與待分解信號之間的投影能量，在整個參數范圍內尋找使投影能量最大的一組參數。具體來說,即：每次循環都將此次運算的投影值與上次循環計算的投影值作比較，選取并存儲較大的投影值及其對應參數，下次循環計算的投影值將與此次比較出的較大值作比較。最終整個循環結束即可得到參數范圍內的最優子波。這種多次嵌套循環結構的運算非常適合進行并行運算，我們選擇的并行思路是：對時間延遲參數的循環執行并行運算，開辟多個內存空間用來存儲各個線程計算的最優子波參數，這樣整個循環結束后每個線程都得到了各自最優的子波參數，最后對比這幾個子波參數，從中選取最終的最佳匹配子波(圖6)。本質上就是把一項工作分成了n份，然后分配給n個線程同時處理。通過這樣的任務分配充分利用了多核處理器的計算能力，有效提高了匹配追蹤的尋優效率。

3.2 多子波時頻分布的并行計算

完成匹配追蹤之后，一個地震信號會被分解為多個Morlet子波，需要計算每個子波的平滑偽Wigner-Ville分布，之后將其疊加得到該信號的時頻分布。然而，對于一個記錄時間較長的地震信號，通過匹配追蹤將會分解得到龐大數目的子波，對所有子波計算平滑偽Wigner-Ville分布的過程將會消耗大量的時間。此過程是一個針對所有子波的循環計算過程，如果基于串行運算，即只有一個線程執行所有的運算，那么所有的子波只能按順序依次計算其時頻分布。然而，基于OpenMP能夠非常方便地開辟多個線程同時計算多個Morlet子波的平滑偽Wigner-Ville分布，這將成倍提高該步驟的計算效率。

圖6 子波尋優并行思路示意圖

4 實際資料試算與分析

為了驗證并行算法的性能，選擇了實際地震資料進行試算并對計算結果進行分析。試算所用計算機的硬件參數為：Intel(R) Core(TM) i7-4710MQ CPU @ 2.50GHz 4核8線程處理器，8GB內存，64位Win7操作系統。為節省試算時間，所用的實際地震數據為截取的部分疊后數據，記錄時長為1000ms，采樣點數為501，采樣間隔為2ms。

4.1 并行前后計算結果對比

利用上述地震數據進行基于匹配追蹤的時頻分析運算。圖7展示的是串行計算的結果。在圖7a中，通過對比原始地震信號(黑色)與通過Morlet子波重構的地震信號(紅色)，可見通過匹配追蹤算法能夠有效地重構原始信號，這也證實了子波分解的有效性； 圖7b為分解得到的Morlet子波； 圖7c為計算得到的時頻譜。為了驗證并行運算與串行計算的計算精度，圖8展示了通過并行運算得到的結果。對比表明，并行運算與串行運算的計算結果一致。為了進一步直觀地展示兩種計算結果的一致性， 將兩種計算結果做差并繪圖(圖9)，從圖中可見，無論是子波分解結果(圖9a)還是時頻譜結果(圖9b)都并無差異。

圖7 串行計算結果(a)原始的(黑色)和通過Morlet子波重構的(紅色)地震信號； (b)分解得到的Morlet子波； (c)計算得到的時頻譜

圖8 并行計算結果(a)原始的(黑色)和通過Morlet子波重構的(紅色)地震信號； (b)分解得到的Morlet子波； (c)計算得到的時頻譜

圖9 串行與并行計算結果差異圖(a)子波； (b)頻譜

4.2 并行加速性能分析

為了描述并行程序的性能，一般采用加速比指標來度量[27]。加速比定義為

式中：n表示并行線程數；t1表示最優串行算法的運行時間；tn表示n個線程并行算法的運行時間。加速比越大，則并行算法加速性能越好。

統計所用地震數據中某一地震道在不同線程數的情況下匹配追蹤和平滑偽Wigner-Ville分布各自的時耗，其中匹配追蹤的子波分解數目為50，同時計算了不同線程數的加速比(表1)。基于OpenMP的并行運算可通過num_threads()函數設定并行計算使用的線程個數。在并行加速性能分析過程中關閉了所有其他程序，保證CUP占用率穩定在3%以內以提高分析的準確性。不過當調用全部8個線程進行并行計算時，操作系統占用的約3%的CPU資源還是會對并行加速性能有輕微的影響，也就是說該算法在8線程并行運算時如果去除掉操作系統的影響還能夠達到更高的加速比，盡管此時受到了操作系統的影響，該算法依然體現出了良好的加速性能。圖10顯示了并行算法的加速比曲線，從中可直觀看到，匹配追蹤與SPWVD具有相似的加速比曲線， 且隨著線程數的增加加速比呈近線性增長。這表明了該并行算法具有良好的加速性能，在計算機硬件條件(可調用線程數)提升時，并行加速比也會隨之以線性增長。

表1 不同線程數對應的計算時間及加速比

圖10 不同線程數的加速比

為了進一步測試并行算法性能，選擇相同的地震道，計算所用線程數固定為8，選擇不同的子波分解數目進行運算時間的統計并計算加速比(表2)。圖11展示了在8線程并行運算時，不同子波分解數目時對應的并行加速比。從該圖可知，在子波分解數目小于20時，并行運算加速比隨著運算量的增大加速比也逐漸增大，當子波分解數超過20后，并行運算的加速比達到飽和，穩定在5.5左右。另外可看到，并行計算Morlet子波的平滑偽Wigner-Ville分布的加速效果相對匹配追蹤的并行運算更加明顯。

表2 不同子波分解數對應的計算時間及加速比

圖11 不同子波分解數時8線程并行運算加速比

5 結論

基于理論分析與試算，認為基于匹配追蹤的時頻分析方法具有可并行性，本文基于OpenMP在多核計算機上實現了該算法的并行運算，有效提高了算法的計算效率，并獲得了與串行算法相同精度的計算結果。

對并行算法的運算時間統計分析表明： 本文并行算法具有良好的并行加速性能，匹配追蹤和多子波時頻譜計算的并行運算都具有近線性的加速比曲線。另外，加速比隨計算量的增加而增加，意味著計算量越大，越能充分發揮并行運算的優勢，當加速比增加到當前線程數所能產生的最大值之后，隨著計算量的增加，加速比將穩定在此峰值附近。

基于匹配追蹤的時頻分析方法適用于并行運算。為進一步提高其運算效率，可嘗試在計算機集群上實現MPI+OpenMP混合的并行算法。另外，基于圖形處理器(GPU)的數值計算日益發揮重要作用，探究基于GPU加速的匹配追蹤時頻分析方法也是有待深入探究的重要課題。

參考文獻

[1] Mallat S G，Zhang Z．Matching pursuits with time-frequency dictionaries．IEEE Transactions on Signal Processing， 1993，41(12)：3397-3415.

[2] 范興利，成谷．基于Morlet小波尺度參數尋優的匹配追蹤時頻分析．中山大學學報自然科學版，2014，53(6)：85-92．

Fan Xingli，Cheng Gu．Matching pursuit time-frequency analysis based on Morlet wavelet scale parameter optimization．Acta Scientiarum Naturalium Universitatis Sunyatseni，2014，53(6)：85-92.

[3] Castagna J P，Sun S，Siegfried R W．Instantaneous spectral analysis：detection of low-frequency shadows associated with hydrocarbons．The Leading Edge，2003，22(2)：120-127.

[4] 孫萬元，張會星，杜藝可．匹配追蹤時頻分析及其在油氣檢測中的應用．山東科技大學學報(自然科學版)， 2011，30(4)：51-57．

Sun Wanyuan，Zhang Huixing，Du Yike．Matching trace time-frequency analysis and its application in oil gas detection．Journal of Shandong University of Science and Technology (Natural Science)，2011，30(4)：51-57.

[5] Partyka G，Gridley J，Lopez J．Interpretational applications of spectral decomposition in reservoir characterization.The Leading Edge，1999，18(3)：173-184.

[6] 鄧志文，趙賢正，陳雨紅等．自適應波形多道匹配追蹤斷層識別技術．石油地球物理勘探，2017，52(3)：532-537，547．

Deng Zhiwen，Zhao Xianzheng，Chen Yuhong et al．Fault identification based on multichannel adaptive waveforms matching pursuit．OGP，2017，52(3)：532-537，547.

[7] 張在金，張軍華，李軍等．煤系地層地震強反射剝離方法研究及低頻伴影分析．石油地球物理勘探，2016，51(2)：376-383．

Zhang Zaijin，Zhang Junhua，Li Jun et al． A method for stripping coal seam strong reflection and low-frequency shadow analysis． OGP， 2016， 51(2)： 376-383.

[8] 段文勝，王鵬，黨青寧等．應用匹配追蹤傅里葉插值技術實現OVT域連片處理．石油地球物理勘探，2017，52(4)：669-677．

Duan Wensheng，Wang Peng，Dang Qingning et al．5D data regularization based on matching pursuit Fourier interpolation for the OVT domain data merging processing．OGP，2017，52(4)：669-677.

[9] Liu J，Wu Y，Han D et al．Time-frequency decomposition based on Ricker wavelet．SEG Technical Program Expanded Abstracts，2004，23：1937-1940.

[10] Liu J，Marfurt K J．Matching pursuit decomposition using Morlet wavelets．SEG Technical Program Expanded Abstracts，2005，24：786-789.

[11] 張繁昌，李傳輝，印興耀．基于動態匹配子波庫的地震數據快速匹配追蹤．石油地球物理勘探，2010，45(5)：667-673．

Zhang Fanchang，Li Chuanhui，Yin Xingyao．Seismic data fast matching pursuit based on dynamic matching wavelet library．OGP，2010，45(5)：667-673.

[12] 張繁昌，李傳輝． 基于正交時頻原子的地震信號快速匹配追蹤． 地球物理學報，2012，55(1)：277-283．

Zhang Fanchang，Li Chuanhui．Orthogonal time-frequency atom based fast matching pursuit for seismic signal．Chinese Journal of Geophysics，2012，55(1)：277-283.

[13] 邵君．基于MP的信號稀疏分解算法研究[學位論文]．四川成都：西南交通大學，2006.

[14] Wang Y． Seismic time-frequency spectral decomposition by matching pursuit． Geophysics，2007，72(1)：V13-V20.

[15] Greg S，Richard B，Yang Xiaoyun．Multicore image processing with OpenMP．IEEE Signal Processing Magazine，2010，27(2)：134-138.

[16] 邱娜．地震子波分解與重構技術研究[學位論文]．山東青島：中國海洋大學，2012.

[17] Morlet J，Arens G，Fourgeau E et al．Wave propagation and sampling theory part Ⅱ：Sampling theory and complex waves．Geophysics，1982，47(2)：222-236.

[18] Taner M T，Koehler F，Sheriff R E．Complex seismic trace analysis．Geophysics，1979，44(6)：1041-1063.

[19] 付勛勛，秦啟榮，徐峰等．基于Wigne-Ville分布和短時傅立葉變換時頻分布計算地震波衰減梯度．新疆石油地質，2012，33(3)：353-356．

Fu Xunxun，Qin Qirong，Xu Feng et al． Estimation of seismic wave energy attenuation gradient based on WVD and STFT．Xinjiang Petroleum Geology，2012， 33(3)：353-356.

[20] 胡廣書．現代信號處理教程(第2版)．北京：清華大學出版社，2015.

[21] Auger F，Flandrin P．Improving the readability of time-frequency and time-scale representations by the reassignment method．IEEE Transactions on Signal Processing，1995，43(5)：1068-1089.

[22] 張顯文，韓立國，王宇等．地震信號譜分解匹配追蹤快速算法及其應用．石油物探，2010，49(1)：1-6．

Zhang Xianwen，Han Liguo，Wang Yu et al．Seismic spectral decomposition fast matching pursuit algorithm and its application．GPP，2010，49(1)：1-6.

[23] 陳永健．OpenMP編譯與優化技術研究[學位論文]．北京：清華大學，2004.

[24] 余濤，朱自強，魯光銀等．基于OpenMP的重力張量并行正演．物探化探計算技術，2013，35(4)：446-451．

Yu Tao，Zhu Ziqiang，Lu Guangyin et al．Parallel forward modeling of gravity tensor based on OpenMP.Computing Techniques for Geophysical and Geoche-mical Exploration，2013，35(4)：446-451.

[25] Pacheco P S著，鄧倩妮譯.并行程序設計導論．北京：機械工業出版社，2013.

[26] 劉揚，王鵬，楊瑞等．基于OpenMP的遙感影像并行ISODATA聚類研究．計算機工程，2016，42(7)：238-243.

Liu Yang，Wang Peng，Yang Rui et al．Research on parallel ISODATA clustering for remote sensing image based on OpenMP．Computer Engineering，2016，42(7)：238-243.

[27] 巫小婷，鄧家先.基于OpenMP的壓縮感知并行處理算法． 計算機應用，2012，32(3)：617-619．

Wu Xiaoting，Deng Jiaxian．Compressed sensing parallel processing algorithm based on OpenMP．Journal of Computer Applications，2012，32(3)：617-619．