基于自適應分數階傅里葉變換的地質雷達直達波抑制方法

周麗軍

（山西省交通科技研發有限公司，山西 太原 030032）

0 引言

直達波抑制問題最早出現在探地雷達應用中[1]，當探地雷達在檢測被掩埋目標時，根據接收到的目標反射信號特點判斷被掩埋目標，發射信號從空氣中入射到掩埋目標時在地面會產生較強的反射信號，這種反射信號將疊加在目標反射信號中形成直達波信號，由于直達波信息強度大，對目標信號檢測造成干擾，因此在提取目標反射信號前需要去除直達波。針對直達波抑制問題，目前常用的有背景相消法、平均對消法、門限法、收發隔離法等[2-4]，利用平均對消法在去除直達波信號的同時，也去除了一部分目標的反射信號，若反射信號較弱則難以保留有效信息。應用背景相消法需提前采集沒有目標的背景直達波信號，再采集有目標的回波信號，從接收回波中減去背景直達波，但考慮到實際環境的不同，在應用中采集背景直達波信號比較困難。還有文獻則根據信號在時間、空間上的不同特征，提出時域濾波、空域濾波、時- 空域濾波相結合的方法抑制直達波[5]。還有學者利用目標與環境在空間與頻域上的差異，采用距離- 方位- 多普勒三域聯合處理方法抑制直達波[6]。

自適應技術在抑制直達波中應用也很廣泛[7-8]，其思想是利用直達波信號與反射信號相關性很低的特征來分離直達波，如有文獻提出利用兩道相鄰波形中直達波的相關性進行自適應對消濾除直達波，但是考慮到接收的直達波信號若受到信道噪聲污染或有較多干擾源，抑制效果將不理想[9]。還有文獻建立了稀疏模型，使用其形態成分分析技術重構有效目標信號與干擾信息[10]。有文獻提出基于經驗模態分解（EMD）或者集成經驗模態分解（EEMD）模型的方法，將信號分成多個單頻信號和殘余信號，從中分離出直達波，此方法適合不平穩信號的處理[11-12]。

分數階傅里葉變換是將信號映射到時頻域空間，實現信號時頻域展開，同時不同階次分數階傅里葉變換對應不同的時頻域，展示出信號從時域逐步變化到頻域的所有變化特征，可以為信號的時頻分析提供更大的選擇余地[13]。有文獻利用分數階傅里葉變換進行動目標檢測與識別，通過選擇合適的變換階數，將最佳分數階傅里葉變換域的信號幅值作為檢測統計量，與門限進行比較后判斷目標有無[14]。有文獻在分數階傅里葉變換域進行濾波從而抑制噪聲域干擾[15]。由于變換角度的可變性，使得信號在不同分數階傅里葉變換中存在多種結果，因此需要尋求一個最優的變換角度，使得變換結果最能體現需求。目前解決這一問題的思想主要集中于最小時間帶寬積原理、最大幅值原理、最大峰度系數原理等[13]。

本文提出基于自適應分數階傅里葉變換的地質雷達直達波抑制方法，根據直達波與目標回波曲線形狀的差異將其映射到分數域空間，逐步調整分數階實現直達波與目標信息在分數域的分離，并建立目標區域與直達波區域能量比最大化的準則，獲取最佳分數階的直達波抑制效果。

1 分數階傅里葉變換

一維信號x(t)的p 階分數傅里葉變換定義為[13]：

式中：exp （-iπpn/ ）2 是分數階傅里葉變換的特征值；φn(t)是歸一化的Hermite-Gaussian 函數。由一維分數階傅里葉變換的定義，可以得到二維信號的分數階傅里葉變換的定義為：

式中：Kp1,p2(s，t，u，v)=Kp1(s，u)Kp2(t，v)是可分離的變換核函數。

隨著p1，p2的選擇不同，可以利用多個階次的分數傅里葉變換對信號進行多分數域分析。然而如何選擇一個最佳的分數值，使得探地雷達回波的直達波在此分數域被抑制，而突出目標的信息，是接下來要解決的問題。

2 最佳分數階選擇

本文在抑制直達波的同時還需盡可能增強目標信號，因此可以考慮使用目標信號幅值能量與直達波幅值能量的比值作為最佳分數域的判定法則。那么對于信號 x(t,s)的（p1，p2）階傅里葉變換 Xp1,p2（u，v），對其進行逆變換F-1{·}操作，能得到此分數域的圖像為：

此時得到的圖像為復圖像，將上述復圖像寫成實部與虛部的形式：

因此，其幅度部分就能寫成：

對幅度信號進行測線方向累積，即：

則yp1,p2(t)表示(p1,p2)階分數域逆變換圖像在時間上的能量值。對于選定的探地雷達天線，能估計出直達波的返回時間，設為t0，計算目標信號幅值能量與直達波幅值能量的比值，即：

3 實驗結果與分析

利用長、寬、高分別為 1.8 m、1.2 m、0.6 m 的沙坑模擬地下環境，測得沙的相對介電常數為4.25，沙表面下方0.2 m 處埋入一水平放置的水瓶模擬充水病害，水管直徑8 cm，長度20 cm。探地雷達為美國GSSI 公司的SIR-20 系列，天線中心頻率為400 MHz，實驗場景與目標放置位置如圖1 所示。由于篇幅所限，本文取第10 道測線的數據進行處理。

第10 道測線的B 掃圖如圖2a 所示，在0～5 ns之間有較強的直達波，在6 ns 左右出現較弱的雙曲線圖像，由于400 MHz 天線分辨率較低，使得目標區域的回波與直達波部分重合。同時從圖2a 中可以看出，在8 ns 附近有來自沙坑底部的回波以及電磁波在水瓶內的多次反射回波形成的干擾信號。

圖2 第10道測線B掃圖與分數階傅里葉變換及其逆變換

圖3 最佳分數階的選擇與其傅里葉變換去直達波效果

圖4 文獻[11]與文獻[12]去直達波方法

對測得的原始圖像進行上述分數階傅里葉變換，當取P1=0.6，P2=0.45 時得到的傅里葉變換如圖2b 所示。在(0.6, 0.45)階傅里葉變換中，呈現出3 個區域，分別映射了原始圖像中的直達波、目標回波和沙底部背景回波。對此分數階傅里葉變換進行逆變換，得到的逆變換圖像如圖2c 所示，比較圖2c 與圖2a，可以看出經過(0.6, 0.45)階傅里葉變換逆變換得到的圖像中，第一層直達波(1 ns 附近)幅值能量降低，而目標區域的幅值能量有明顯的提升。對圖2c的圖像繼續做分數階傅里葉變換以及逆變換，式（7）的能量比值，對分數階P1與P2遍歷[0 1]區間值，取能量比值最大的分數階為最佳分數階傅里葉變換，如圖3a 所示，在(0.2, 0.95)階時獲得最大能量比，由此得到的自適應分數階傅里葉變換的逆變換圖如圖3b 所示，比較圖3b 與圖2a，能明顯看出直達波已被去除，同時也去除了沙底反射的背景回波，而且增強了目標區域能量。

將本文方法與平均對消法[11]、EMD 去除直達波方法[2]進行比較，如圖4 所示，利用平均對消法去除了0～4 ns 處的直達波，但同時也減弱了目標區域能量，并且增強了電磁波在目標內的多次反射回波。利用EMD 方法增強了目標區域能量，去除了直達波，但也增強了電磁波在目標內的多次反射回波。而本文方法則有效去除了多次反射回波，進一步增強了目標區域能量。

4 結論

本文提出基于自適應分數階傅里葉變換的地質雷達直達波抑制方法，將探地雷達接收回波中的水平線狀直達波與雙曲線狀目標回波映射到分數域空間，進行分數階傅里葉變換，實現信號在分數域的時頻展開，通過多次變換逐步分離直達波與目標信息。建立重構圖像中的目標區域與直達波區域能量比最大化的準則，獲取最佳分數階的直達波抑制效果。在最佳分數階傅里葉變換處理下，既抑制了直達波，同時也增強了目標回波。相比于平均對消法與EMD方法，本文提出的方法除了抑制直達波，增強目標回波，還能有效去除多次反射產生的干擾。