基于EMD 的井中雷達信號預處理

鄒 寧，金楊超，郭 成，陶 杉，宋 海

(1. 中國石油化工集團公司碳酸鹽巖縫洞型油藏提高采收率重點實驗室 烏魯木齊 830011；2. 中國石油化工股份有限公司西北油田分公司 烏魯木齊 830011；3. 電子科技大學資源與環境學院 成都 611731)

井中雷達系統尚未有統一的國際化標準[1-4]，井中雷達數據傳輸速度有限且數據很容易因為復雜的環境而被“污染”，如何從原始數據中提取出有效數據，并且對井中儀器進行實時控制及監控[5-7]，是井中雷達系統數據處理所面臨的首要問題。

在信號預處理方面，井中雷達信號因為工作環境復雜，遇到的干擾十分劇烈，同時由于其非線性、非穩態、超寬帶的特點，導致其使用傅里葉變換、小波變換等常規分析方法難以實現信號分析[8-10]。井中雷達信號預處理工作需要將雷達對地下異常目標體的回波信號進行提取分析，由于目標反射較弱，或是探測深度較深導致回波信號衰減較大，導致信號強度很弱，同時又由于其反射環境復雜，容易受系統、井壁、井液等影響，導致雜波干擾非常劇烈[11-12]。因此，如何將復雜的回波數據分離出噪聲和有用信號，獲取真實雷達反射回波，是信號預處理甚至雷達硬件系統設計工作所面臨的最大問題，具有重要的研究意義。

文獻[13]將低維離散余弦變換(discrete cosine transform, DCT)特征向量與支持向量機(support vector machines, SVM)分類方案結合使用，采用Wiener濾波技術對探地雷達信號進行了自適應預處理，研究結果表明該方法能有效增強信號的分類精度，在計算效率和準確性方面能取得良好的性能。文獻[14]發現探地雷達信號中有較強的散射波以及由色散引起的傳播噪聲，因此采用二維頻域響應濾波器和基于傳播常數的濾波來抑制其散射，就信號幅度水平而言，該濾波器在特定應用下信噪比有一定提高。文獻[15]提出了一種反濾波方法用于消除探地雷達信號的子波衰減和頻散，該方法認為地層介質的邊界條件系數是隨機數，充分利用了最小相位濾波器作為其等效濾波器，有效提高了探地雷達的探測深度，并且能夠進一步提高探測分辨率。

另一方面，通過EMD (empirical mode decomposition method)對信號進行模態分解，這種分解方式不同于其他濾波和小波變換等方法，EMD 分解對信號具有自適應性，它的優點是不會運用任何已經定義好的函數作為基底，而是根據所分析的信號自適應的生成固有模態函數[16]。文獻[17-18]先后提出了EEMD (ensemble empirical mode decomposition)算法、高維集合EMD 分解(modified ensemble EMD, MEEMD)，通過向原始數據中加入噪聲，能夠有效減小經驗模態分解中存在的混疊問題，同時也將其推廣到高維處理。文獻[19-25]充分利用了模態經驗分解的特性，根據EMD 分解進行了多分辨率以及時頻分析方法的研究，進一步拓展了EMD 分解技術在探地雷達信號處理中的降噪方法和應用前景。

以上方法僅運用到探地雷達的信號處理中，并沒有應用到井中雷達。本文基于EMD 分解的信號預處理方法，對希爾伯特-黃變換算法進行了闡述，通過構造函數的方式驗證算法，并對算法進行了井中雷達的信號預處理。

1 井中雷達信號特征

1.1 井中雷達信號模型

井中雷達的基本工作原理如圖1 所示。其接收信號可看作由發射信號與系統各個模塊、反射信號以及周圍環境等作用產生的，通過式(1)進行線性模型描述：

式中，d(t)表示發射和接收天線間的直達波分量；fs(t)是經過脈沖源產生的信號在天線上產生的激勵信號，會由于信號拖尾、波形畸變等對回波信號產生一定的影響；fa1(t)是該脈沖信號經過發射天線時所產生的沖激響應；fa2(t)是接收端接收天線對直達波信號的沖激響應；fc(t)是由于收發天線間的互耦效應所產生的互耦響應，主要包含了天線互耦作用fc1(t)以及井內泥漿和雷達儀器周圍沖洗帶影響fc2(t)。一般而言，井中雷達在實際工作中，雷達收發天線間的相對位置固定，直達波僅與儀器周圍、收發天線之間的介質情況有關。因此，在探測過程中，直達波分量僅受井中填充井液影響，在一定井深范圍內直達波分量是相對穩定的。

c(t)主要是電路耦合分量，僅針對收發一體的井中雷達儀器，由于收發一體的井中雷達接收天線與發射天線合而為一，通過 TR 組件對天線收發狀態的切換，導致源與采集之間隔離較為困難使電路耦合分量，隨電路工作穩定性呈一定波動。s(t)為雷達波信號經過地層傳輸，遇到目標反射而被接收端獲取的信號回波分量，也是井中雷達探測所感興趣的信號回波分量：

式中，fg(t)是反射、折射、散射所引起的地層沖激響應；ft(t)為探測地層中目標的沖擊響應，一般是由地層中孔洞、裂縫、埋藏物等產生，受目標的大小、形狀、電磁特性以及傳播地層的電磁特性影響，從而表現出不穩定的波形，產生目標物的相應波形，是井中雷達探測的目標分量，而雷達信號預處理所要達到的目的便是獲取該信號的所有分量。

n(t)為加性噪聲分量，其來源較為復雜，而因為井中雷達的特殊性，其主要來源是系統熱噪聲。不同介質由于電磁特性不同，所得到的回波信號也不同。井中雷達回波信號受到的干擾來源復雜，極易受到井中及井周復雜的介質環境影響，且由于井中雷達器件的局限性，其接收到的信號信噪比較低，受到井周地層介質的干擾強烈，波形畸變更為嚴重。因此，通過合適且針對性較強的信號預處理手段來進行波形修正，降低環境、器件干擾，提高系統信噪比，進行井中雷達信號預處理從而獲得更好的雷達探測結果。

1.2 井中雷達數據類型

井中雷達數據由3 種形式構成：A-掃、B-掃、C-掃，分別對應雷達探測一維、二維、三維數據，表現形式分別是回波波形、波形矩陣灰度圖、灰度圖陣列。

井中雷達A 掃數據由一次回波采樣構成，能夠描述在某一井深處，垂直井向的特定水平方向內的回波信號。其數學表達式使用柱坐標系描述為：

式中，j、k為固定值，指示該回波信號的方向角與深度；i為時間采樣點數；N為時窗長度。

B 掃數據是由一定井段范圍內的多個A 掃數據組成，描述的是在井段范圍內，某一特定水平方向內的回波信號。其數學表達式使用柱坐標系描述為：

式中，j為固定值，指示特定水平方向；k為深度采樣點數；M為某深度范圍內采樣波數。

C 掃數據是井中雷達多個通道、不同方向，天線在相同井深處，多個方向內的回波信號，由雷達儀器在固定深度處由固定方向的天線旋轉采樣多次回波信號獲取，或者由多個固定方向接收天線同時接收獲得。其數學表達式使用柱坐標系描述為：

式中，j為雷達通道編號；P為雷達儀器通道數。

2 基于EMD 的信號預處理

井中雷達的數據后期處理流程主要分為數據預處理、信號預處理、成像與反演3 個階段。其中，信號預處理主要是將信號進行去噪、濾波、去直流分量、增益控制處理，目標是剔除所有干擾信號，僅留下有用的回波信號。常規的方法有平均濾波、帶通濾波、滑動窗口濾波、自動增益控制、小波變換等。

井中雷達預信號處理過程中，滑動平均濾波有時能更好地去噪，有效抑制周期性和漸變狀的干擾，如由系統溫度引起的采樣時鐘偏移，但滑動平均濾波不適用于由脈沖干擾所造成的噪聲，即由時域波形的脈沖采樣值誤差，如采樣噪點、環境干擾等，因此不適用于隨機強噪點干擾嚴重的場合。同時，與去背景噪聲算法一樣，對平行于井向的平面或線性信號是無能為力的，需要針對不同的場景來應用不同的處理方法。一般而言，井中雷達信號在傳播過程中，幾乎不會產生頻率偏移，因此，根據雷達的發射源的頻率情況，對接收數據進行相應的帶通濾波，濾波系數主要取決于發射信號的頻譜。

但自然信號通常都是非線性和非平穩的，因此，在處理這些非線性和非平穩時域信號時，無論是傅里葉變換還是小波變換都會對信號引入線性分解，不是最佳方法。相較而言，對于井中雷達信號而言，由于信號的強非線性以及非平穩性，EMD無疑是非常有效的分析方法。通過對信號的分解及時頻域分析，對信號根據其信號自身的特點，自適應地分解信號，為信號的時頻域分析提供核心思路。

希爾伯特-黃變換算法是目前熱門的信號分析算法，屬于一種經驗方法，其主要過程分為模態經驗分解EMD 算法和希爾伯特譜分析兩個部分，通過對信號進行時域分解，從時頻域對信號進行詳盡的解析。

2.1.1 EMD 算法

EMD 的基本思想主要在于對所有復雜信號而言，其時間序列都可以分解成具有不同特征尺度的有限個本征模函數(intrinsic mode function, IMF)，即IMF分量。每一個IMF 具有局部正交、自適應的特點，用于分析非平穩信號時，有效避免了不同分量之間隨時間變化的相同頻率部分互相干擾，大大增加了非平穩信號的分析預測和解釋能力。

EMD 的原理來源于一個簡單的假設，即任何信號都由不同的IMF 組成，每個IMF 在一個分離的時間尺度上代表一個嵌入的特征振蕩。EMD 分解具有的局部自適應性，適用于非線性以及非平穩信號的分析。EMD 的出發點是考慮局部振蕩水平上的振蕩信號，并對作為新信號的慢振蕩分量進行迭代，如式(7)所示：

2.1 希爾伯特-黃變換算法

計算流程如圖2 所示。

圖2 EMD 分解迭代流程圖

該方法旨在對非平穩和非線性信號進行分解，以展示其隱含的準周期性和特征。它是一種局部自適應方法，以數據作為驅動，因此更適合于非線性和非平穩數據的分析。

常用的數據分析方法中，大部分基于傅里葉變換進行，但用傅里葉變換進行數據分析時，通常需要滿足兩個條件：信號必須是平穩的，并且在給定的時間跨度內是線性的，然而一般來說，物理信號不滿足這些條件。EMD 分解方法將信號分解為振幅和頻率調制函數，也就是固有模態函數IMF。EMD 分解可以在自適應時間-頻率-幅度空間中處理非線性和非平穩信號的數據，將信號進行分解為IMF 和殘差，給定方程如下：

式中，I(n) 是 多分量的待分解信號； IMFm(n)是第m個固有模態函數； R esM(n)為對應于M個IMF 的殘差，或者稱為趨勢(mean trend)。

EMD 算法計算過程概述如下，以s(t)作為輸入信號：

1) 找出s(t)局部極大值和極小值點，獲得極大值點和極小值點的序列，分別將極大值點和極小值點進行3 次樣條插值，此時，獲得s(t)的上包絡線和下包絡線，兩個包絡線互相不對稱；

2) 對極大值包絡線和極小值包絡線取平均值，此時獲得了包絡均值，將原始信號減去均值包絡線，此時獲得第一個分量：

3)將h1(t)作為輸入重復進行步驟1)、步驟2)獲得h2(t)，依此類推；

4)重復上述步驟，使用相鄰兩次分解分量的標準差作為停止準則，一般取值為 0 .2 ≤S≤0.3，S計算如下：

此時獲得hp(t)，當再次重復步驟1)、2)時，將會獲得均值包絡為x=0的一條直線或接近零的直線。此時，hp(t)的上下包絡線具有良好的對稱性，繼續篩選下去，信號不會有太多的變化，因此停止篩選。

5)當滿足停止準則時，hp(t)滿足IMF 分解條件，此時第一次EMD 分解過程完成，hp(t)為第一個IMF，記為c1(t)。 將原始輸入信號減去c1(t)，獲得剩余量r1(t)：

6)將r1(t)作為新的信號輸入，重新執行步驟1)～步驟5)，獲得新的余量r2(t)以及第二個模態c1(t)， 如此重復n次，依次獲得c3(t),c4(t),···,cn(t)，此時剩余的余量rn(t)成為單調函數，無法再進行IMF 分解，整個EMD 分解過程完成。

由上述步驟可得：

此時，原始信號被分解成為n個模態分量ck(t)以及剩余一個單調趨勢函數分量rn(t)。至此，模態經驗分解完成。

2.1.2 EMD 時頻域分析

通過EMD 分解后，將原信號分解為一組單分量信號IMF 的組合，同時對算法進行模擬信號的分解驗證，接下來對其分解結果進行進一步的希爾伯特譜分析。

將各IMF 分量分別進行希爾伯特變換，記為Hck，此時就得到了各個IMF 的瞬時特征量，因此原信號可以表示為：

也即希爾伯特譜，可以看出，希爾伯特譜與傅里葉譜相比，傅里葉譜的推廣更具有一般化的特點，希爾伯特譜更具有頻率的瞬時特性。

此時對每一個IMF 分量進行解析信號的構造，構造函數為：

最終，通過EMD 分解的解析函數的構造提取，即可繪制出原始信號的“時間-頻率”分布，并且在每個時間點上，都能找到每個平穩態分量的頻率值。

2.2 BEMD 及其快速算法

EMD 分解的二維表現形式為二維模態經驗分解(BEMD)，也叫2D EMD。BEMD 同樣需要尋找B 掃數據的局部極大值點以及極小值點，分別對極大值點和極小值點進行后續插值，獲得極大值曲面和極小值曲面，而提取每一個BIMF 也需要多次迭代，其計算流程與一維EMD 分解算法大致相同。

由于井中雷達特性所造成的非線性、不平穩、頻帶廣的干擾噪聲，對于B 掃數據使用局部適應的二維EMD 分解方法，可以有效去除非線性干擾，進一步保證數據質量。二維EMD 分解主要用于雷達B 掃數據濾波去噪，通過對雷達B 掃數據進行二維EMD 分解得到不同尺度特征的BIMF 分量。由于在提取過程中引入了極值點距離，從而定義了局部尺度，同時在分量提取時，由于篩選條件的限定，由BEMD 方法獲取的圖像，分解成為一組由“精細”至“粗糙”的圖像組。首先，分解出的是比較細節的部分，包含了較為高頻的分量圖像信息，之后分解出的圖像逐漸尺度變大，顯示了較為低頻的分量圖像信息，最終剩余部分是最大尺度的趨勢分量圖像信息。BEMD 方法對雷達B 掃數據具有一定的自適應性，計算過程大致與一維相似。

二維EMD 分解涉及到的插值方式與一維較為不同，本文使用的二維插值方式為多重二次曲面(multi-quadric, MQ)插值，MQ 插值方案非常適用于二維包絡插值，該插值方式主要用于對地理表面以及重力和磁場表面進行插值，已被廣泛應用于地球科學和其他領域。MQ 方法是一種使用徑向基函數(radial basis function, RBF)的散點數據全局插值方案。對于散點數據(xi,yi)極值，MQ 方法使用方程式：

由此可獲取上下限的包絡值。

3 基于EMD 的仿真驗證

3.1 EMD 算法試驗

為驗證上述EMD 分解算法以及直觀反映其分解步驟和效果，嘗試對構造函數進行EMD 分解，模擬非平穩信號EMD 分解具體流程和結果。構造函數：

n取值范圍為0 ≤n≤500，采樣點數為500，其對應頻率約為5.00 Hz、7.14 Hz、83.00 Hz、此時信號函數圖像及其頻譜如圖3 所示。

圖3 構造函數信號波形圖

對上述構造函數進行EMD 分解，如圖4 所示。

圖4 構造函數EMD 分解

模態 3 分量的幅度值較小，頻率更低，表現出的特征更接近于信號的趨勢性分量，對信號重構影響較小。這樣，通過時域自適應分解算法，將信號分解為不同的頻率分量，有效分解為不同的特征尺度，便于進一步分析信號。

對該信號進行希爾伯特變換，并解析信號的構造，分解為瞬時頻率，繪制時頻分布如圖5a 所示。

圖5 構造函數的EMD 分解時頻分布圖

與圖3b 對比，可以看出經過EMD 分解，成功將擬合信號分解為4 個分量，其中，imf1 對應原信號83.00 Hz 分量；imf2 對應原信號7.14 Hz 分量；imf3 對應原信號5.00 Hz 分量；res.為信號殘差趨勢分量。各分量能夠較好地與原信號頻域對應，并且獲取了各個分量的時頻譜，能夠進行各個分量的時變頻譜分析。

3.2 BEMD 快速算法試驗

為驗證BEMD 快速算法的合理性及處理效果，建立了二維FDTD 模型進行仿真計算，其中，網格區域為8 m×8 m， 地層介電常數為εr=7，電導率為 σ=1×10-4S/m，接收天線和發射天線間距0.2 m，發射頻率為80 MHz 零階高斯脈沖，雷達探測深度為8 m，分辨率為0.7 m，目標為半徑0.3 m金屬圓形目標，時窗選取T=2×10-7s。通過FDTD正演仿真技術，建立FDTD 網格進行正演仿真模擬，對仿真結果進行加噪處理，以分析BEMD 對其處理效果。加入具有標準正態分布，即以0 為均值、1 為標準差的二維高斯噪聲，經過分析可得，通過BEMD 能夠分解B 掃數據，更清晰地將目標回波信號反映出來。將imf2、imf3 數據進行疊加，可以看出，通過BEMD 分解和對其分量的重組，使得目標的回波雙曲線更加明顯，特征更加突出，提高了B 掃數據的信噪比，有助于進一步的雷達信號成像與反演。通過對仿真數據進行分解，驗證了BEMD 對于井中雷達B 掃數據處理的降噪能力。

4 基于EMD 的實測數據處理

4.1 井中雷達數據EMD 分解

基于理論和仿真計算，對實測數據進行EMD分解，采用的原始數據為油井測試數據，此時并沒有對原始數據設定分解信號的基函數，通過信號的自適應包絡曲線以及迭代，將信號分解為如圖6 所示的多個不同特征的imf 分量和殘差。

圖6 井中雷達信號EMD 分解

其中各IMF 分量標注如圖，該信號共分解為7 個IMF 模態分量以及一個單調趨勢函數分量res，能看到每一次imf 分解完畢，獲得的模態分量在時窗范圍內是平穩的，其分解得到的特征模態分量從上至下是由高頻到低頻的順序進行排列，即imf1是信號的高頻主成部分，依次頻率減小，res.為頻率接近0 的殘余分量。

模態特征函數分量從上至下是嚴格由高頻到低頻的順序排列，此處并非 imf1 頻率要高于imf2，而是允許后者局部頻率高于前者局部頻率，但在同一時間采樣點處，頻率由高到低排列。如圖7 所示，框內標注處的分量局部信號，包含有較低頻率的局部信號。

圖7 EMD 分解的局部信號

這也正好反映了EMD 分解的局部性強的特性，即相鄰的imf 可能包含有相同頻率范圍內的震蕩，但在相同的頻率范圍內的震蕩不會出現在同一時窗范圍內。

將上述信號進行時頻域分解，由于信號的不穩定性，分解之后的時頻譜起伏較為劇烈，隨時間采樣點數增加，頻率變化范圍較大，頻率起伏范圍可達300 MHz。通過時頻譜可以看出，imf1、imf2、imf3 主要包含的是頻率為100～300 MHz 分量，與井中雷達設計的系統工作頻率相符合，因此取imf1、imf2、imf3 這3 個分量進行疊加，疊加之后的信號形成B 掃圖像。與原信號B 掃圖像對比可得，對原B 掃數據按行模態經驗分解并疊加imf1、imf2、imf3 后，合成的B 掃信號能夠更多的對地質異常目標進行反映，異常目標更清晰明顯，異質構造特征更突出，對分析地層異常提供更合理有效的解釋。從分解過程也容易看出，EMD 分解還具有線性、完備性的特點。

4.2 井中雷達數據BEMD 分解

對實測井中雷達數據進行BEMD 分解，所選取的數據為中石化榆林大125 井約1 975～2 995 m深井中雷達通道1 數據，B 掃數據去除背景后通過對其進行BEMD 分解，同時疊加imf2、imf3 分量，能夠很好地剔除原始數據中由于硬件采集模塊帶來的采樣噪點以及其他各部分系統模塊帶來的系統噪聲，減小高頻干擾，同時通過分離低頻分量，剔除低頻干擾。因此，地層波動情況得以更清晰地反映出來，使得回波雙曲線更加突出，對于分析地層的不同情況提供更好的信息量支持。

由上述討論可知，基于EMD 分解的井中雷達數據處理通過模態經驗分解算法，根據信號本身的時間尺度特性對信號進行分解，有效地克服了傅里葉變換的局限性，更適應于井中雷達信號非線性、非平穩的特點，可以在任何時刻獲得信號的頻率分布，在時頻域上提供更高的頻率分辨率。可見，在井中雷達數據分析處理方面基于EMD 分解的方法優于小波變換，該方法很好地解決了雷達數據的模式混合問題，保證了分解對噪聲的穩定性。針對不同類型儲層，特別是縫洞型碳酸鹽巖儲層這種地層結構復雜的情況，基于EMD 分解的方法相對于其他常規信號處理技術，可以獲取到更多地質信息。但其模式的分裂問題仍然是一個開放的問題，更重要的是其數學分析并不完整，在理論分析和理解、性能增強、優化等方面仍存在很多問題。

5 結 束 語

本文簡要介紹了井中雷達信號模型和數據類型的數學描述以及常規信號預處理方法，對井中雷達信號預處理的理論進行了推導研究。對希爾伯特-黃變換方法進行了詳細推導，也進行了基于EMD 信號預處理方法的理論研究。通過仿真的方式驗證了將EMD 的一維方法推廣到了實測二維B 掃數據方法的正確性、有效性。通過實測數據對各算法進行驗證，其中，EMD 有效提高了采集原始數據的信噪比；BEMD 能夠有效剔除高頻、低頻的干擾，較一維的模態經驗分解有較大的改進。結合仿真和實測數據的處理結果可以看出，基于EMD 的時頻域方法能夠有效對井周地層解釋提供依據，具備應用到井中雷達數據處理的潛力。

本文研究工作得到了中國石油化工股份有限公司西北油田分公司項目(34400007-20-ZC0607-0050)的支持，在此表示感謝。