基于探地雷達雜波抑制與偏移成像的樹木根系定位方法

李光輝 徐 匯 劉 敏

(江南大學人工智能與計算機學院， 無錫 214122)

0 引言

根系檢測是樹木健康狀況評價的重要手段，在古樹名木和果樹養護管理等方面都能發揮重要作用。現有的根系檢測方法大致分為破壞性檢測和無損檢測[1](Non-destructive testing, NDT)兩類。傳統的破壞性檢測方法包括根鉆法、土芯法、土壤剖面法和全根挖掘法[2-3]等，這些方法勞動強度大且費時費力，不適合大規模應用，還可能對周圍的根系造成不可逆轉的破壞。而X射線斷層掃描、核磁共振方法、聲學方法和電阻率層析成像[4-8]等檢測方法能夠實現根系參數無損測量。探地雷達(Ground-penetrating radar, GPR)作為一種新興的地球物理探測方法，已經被應用于樹木根系和樹干的無損檢測[9-12]。與其它無損檢測方法相比，探地雷達具有易于操作，野外便攜，可以掃描大型根系并且能夠快速、經濟地進行重復測量等優點[13-15]。

探地雷達沿掃描路徑向地下發射電磁波，通過接收到的反射回波來繪制掃描區域的地下高分辨率雷達圖(B-scan圖像)。但是，探地雷達B-scan圖像通常會因為發射和接收天線之間的直接耦合、地表反射和非平坦地形的散射響應引起的雜波造成模糊[16]。因此，探地雷達B-scan圖像應用于任何地下物體檢測算法之前，必須先進行雜波抑制。最簡單常用的方法是均值減法(Mean subtraction, MS)，但是目標響應的強度會受到均值減法的影響而減小。子空間方法，如奇異值分解(Singular value decomposition, SVD)、主成分分析(Principal component analysis, PCA)和獨立成分分析(Independent component analysis, ICA)，將B-scan圖像視為幾個子空間的并集，通過去除雜波子空間來減弱雜波[17-20]。但是，子空間方法需要對雜波子空間進行精確測量和劃分才能獲得理想的雜波抑制效果。魯棒主成分分析(Robust principal component analysis, RPCA)通過優化低秩和稀疏矩陣表示問題來分離雜波和目標響應，在真實數據上表現優于子空間方法[21]。

通過偏移技術對雜波抑制后的B-scan圖像進行聚焦以實現目標定位。探地雷達常用的偏移算法包括后向投影算法(Back projection, BP)、克希霍夫(Kirchhoff)積分法、逆時偏移(Reverse time migration, RTM)、頻率-波數域偏移(F-K偏移)[22-25]等。F-K偏移在頻率-波速域中采用快速傅里葉變換(Fast Fourier transform, FFT)進行求解，相對其他3種偏移方法，具有更高的計算效率，實際應用廣泛。然而偏移成像聚焦效果與偏移速度密切相關，對于點目標來說，當偏移速度與實際速度一致時，能量匯聚到同一繞射點上，偏移能夠完全聚焦，如果偏移速度小于實際速度，偏移就不能完全聚焦，目標形成的雙曲線不能完全消除，如果偏移速度大于實際速度，雙曲線開口方向有反轉的趨勢。為有效利用B-scan圖像定位樹木根系，本文提出基于RPCA和深度自動編碼器結合的雜波抑制方法，使用直接最小二乘法擬合目標回波雙曲線估算介質相對介電常數求取偏移速度，最后通過改進插值的F-K偏移實現根系定位。

1 研究方法

1.1 魯棒深度自動編碼器

從數學上講，觀測到的探地雷達B-scan數據X可以由一個表示直達波等背景信號的低秩矩陣L，一個表示目標回波信號的稀疏矩陣S和一個隨機噪聲矩陣N聯合表示[26]，即

X=L+S+N

(1)

魯棒主成分分析利用凸優化求解式(1)。

(2)

式中 ‖·‖*——核范數

‖·‖1——l1范數

‖·‖F——F范數

λ——控制稀疏矩陣稀疏性的懲罰系數

ε——任意正數

核范數用于對矩陣的奇異值求和來獲取低秩分量，RPCA用迭代奇異值的線性方法來求解核范數，l1范數用于計算矩陣元素的絕對值之和來獲取稀疏分量。本文應用RPCA和深度自動編碼器相結合的魯棒深度自動編碼器(Robust deep autoencoder, RDAE)來抑制探地雷達B-scan圖像的雜波。典型的自動編碼器是一個3層結構的神經網絡。如圖1所示，自動編碼器由包含n個節點的輸入層、h個節點的隱藏層和n個節點的輸出層組成。自動編碼器的目標是最小化輸入數據x與輸出數據x*之間的重構誤差，因此可以將最小化重構誤差作為自動編碼器的損失函數，即

圖1 自動編碼器結構圖Fig.1 Architecture of autoencoder

(3)

式中 ‖·‖2——l2范數

E(·)——編碼器D(·)——解碼器

具有多個隱藏層的自動編碼器被稱為深度自動編碼器[27]，每增加一個隱藏層需要增加一對編碼器和解碼器，通過多個編碼器和解碼器，深度自動編碼器可以有效地表示輸入數據上的復雜分布。RDAE獲取稀疏分量的方式類似于RPCA，采用l1范數，但是RDAE繼承了深度自動編碼器獲取低秩分量的非線性能力，在不規則雜波條件下優于RPCA，RDAE將式(2)修改為優化問題

(4)

1.2 直接最小二乘法

一般地，當探地雷達天線掃描經過樹根上部土壤時，樹根反射雷達脈沖形成的回波雙曲線可以表示為

(5)

式中x——天線位置

t——雷達波的雙向傳播時間

c——真空中光速

εr——土壤的相對介電常數

(x0,t0)——雙曲線頂點坐標

在式(5)中，樹根直徑被忽略，這是可以接受的，因為本文中的樹根直徑比雷達波長小。土壤的相對介電常數采用直接最小二乘法(Direct least square, DLS)擬合雙曲線計算[28]。相較于霍夫變換(Hough transform, HT)，DLS對不規則雙曲線擬合效果更好。圖2為執行DLS算法的流程圖。為了說明該過程，給出了圖3所示的DLS方法處理仿真B-scan圖像的示例。圖3a、3b為進行零點校正和雜波抑制后的B-scan圖像，轉為灰度圖以獲取聯通區域，采用Canny算子提取的聯通域邊緣擬合圖3c中標紅的上邊緣。在DLS循環的每次迭代中，從上邊緣隨機選取10個坐標點，記錄每次求解的參數(x0,t0,εr)，使用DLS擬合的雙曲線見圖3d。本文已利用仿真數據對DLS的準確性進行了測試，該仿真涉及具有不同相對介電常數的均質土壤。在用DLS擬合雙曲線的循環迭代過程中，當待擬合的點到雙曲線的誤差平方和(Sum of the squared errors, SSE)小于給定閾值時，終止迭代。

圖2 基于DLS的雙曲線參數估計算法流程圖Fig.2 Flow chart of DLS-based algorithm for hyperbola parameter estimation

圖3 基于DLS算法擬合雙曲線示例Fig.3 Examples of hyperbola fitting with DLS algorithm

1.3 F-K偏移

根據惠更斯原理，如果電磁波的波前(Wavefront)在t時刻的位置是已知的，將t時刻波前上的每個點看作t+Δt時刻波前的子波源，如圖4所示，t+Δt時刻的波前可由t時刻波前疊加重建。在探地雷達檢測中，電磁波從空氣傳播進入地下界面，將地下界面的各個反射點看作子波源，對接收到的記錄剖面作時間反演，即為探地雷達的偏移成像。

圖4 惠更斯原理圖Fig.4 Schematic illustration of Huygens principle

設探地雷達B-scan記錄剖面為p(x,z=0,t)，z是垂直坐標，向下為正，P(kx,z=0,ω)是偏移后剖面p(x,z,t=0)的二維傅里葉變換，即

P(kx,z=0,ω)=?p(x,z=0,t)e-j(kxx+ωt)dxdt

(6)

將式(6)在頻率-波數域做波場外推，深度z處的波場可以表示為

P(kx,z,ω)=P(kx,z=0,ω)ejkzz

(7)

由于探地雷達記錄電磁波的雙程走時，因此可以把波速視為介質中波速v的一半，根據色散方程，kx、kz和ω的關系可以表示為

(8)

設p(x,z,t)是P(kx,z,ω)關于kz、ω的二維傅里葉逆變換，代入色散方程式(8)，同時令t=0可得偏移后的圖像。

(9)

為了提高偏移圖像的分辨率，在F-K偏移的實現上改進了時頻插值方法[29]。首先，對雜波抑制后的B-scan數據在t方向進行零填充來提高頻率域的分辨率，使得ω到kz域的插值誤差減小，然后在快速傅里葉逆變換前對kx域中數據進行零填充以提高x方向上的空間采樣率。改進后的F-K偏移流程如圖5所示。

圖5 改進的F-K偏移流程圖Fig.5 Flow chart of improved F-K migration

2 實驗與結果分析

2.1 實驗數據采集

實測數據采集使用美國GSSI公司的SIR-3000型探地雷達，雷達天線包含400、900 MHz 2種頻率，如圖6所示。實驗數據分成實測數據和仿真數據2組。實測場地位于江南大學西北操場沙坑，選擇人工填埋的香樟粗枝代替真實根系，開展根系定位的模擬實驗，檢測區域的土壤類型主要為細沙，細沙材質均一，易于模擬且相對介電常數與土壤較為接近，實測實驗期間未發生降雨，細沙的含水率穩定。仿真數據使用基于時域有限差分(FDTD)方法求解麥克斯韋方程的開源軟件gprMax[30]獲得，gprMax允許CPU并行計算，并且支持GPU加速運算。gprMax中大多數命令是可選的，但必須指定模型尺寸，空間離散化步長和時窗大小。本文實驗數據處理平臺搭載Intel Core(TM) i7-9700處理器，主頻3.0 GHz，內存16 GB，Nvidia GTX1660顯卡，6GB顯存。開發環境為Matlab 2019a、Python 3.8和深度學習框架Pytorch 1.5。

圖6 GSSI SIR-3000型探地雷達Fig.6 GSSI SIR-3000 model ground-penetrating radar

2.2 雜波抑制

首先使用仿真模型對提出的方法進行測試，仿真模型模擬了城市路面的分層結構，如圖7所示，仿真的分層模型由15 cm厚的瀝青層、15 cm厚的水泥層和50 cm深的干燥黏土組成，水泥層的粗糙上下表面高度在3 cm內變化，埋在土壤中的根系半徑為4 cm，探地雷達沿掃描方向進行探測，仿真基本參數見表1。

圖7 城市路面結構的仿真模型Fig.7 Simulation model of urban road structure

表1 仿真基本參數Tab.1 Simulation parameters

使用改善因子(Improvement factor, IF)定量分析雜波抑制前后圖像信雜比(Signal-to-cutter ratio, SCR)的改善程度，IF值越大，雜波抑制方法性能越好。IF定義為

(10)

其中，F為改善因子，Sbefore和Safter分別是應用雜波抑制方法前后B-scan圖像的信雜比，信雜比定義為

(11)

式中Ns——目標信號區域Rs中的像素數

Nc——雜波區域Rc中的像素數

I(p)——B-scan圖像第p個像素的信號強度

如圖8所示，由于水泥層粗糙不平表面的影響，許多不均勻的雜波無法通過MS和SVD很好地與目標回波分離。RPCA和RDAE通過低秩和稀疏矩陣分解，能更好地分離雜波并保留目標回波。RDAE中深度自動編碼器共包含7層，每層的節點數為512、128、32、8、32、128、512，激活函數為Sigmoid，以便深度自動編碼器學習輸入的非線性表示。表2給出了4種方法的IF值，去除2個最大主成分的SVD性能比MS略好，RPCA和RDAE能更好地分離目標信號與雜波，但RDAE具有更高的IF值。

圖8 不同方法對仿真數據的雜波抑制結果Fig.8 Clutter suppression results for simulated data with different methods

表2 仿真和實測數據的IF結果Tab.2 IF results of simulated and real data

實測數據通過預埋樹根的方式采集，香樟樹枝預埋深度為0.3 m，填埋介質為沙子。天線頻率設置為400 MHz，天線移動步長0.01 m，沿著4 m長的測線進行掃描。獲得的探地雷達B-scan圖像包含313×409個像素點。掃描區域的原始探地雷達B-scan圖像經零點校正后如圖9a所示。觀察到目標回波信號出現在10 ns內，因此在雜波抑制前對原始數據進行時變增益，降低其在10 ns往后的信號強度。圖9b～9e分別為使用MS、SVD、RPCA和RDAE方法對圖9a進行雜波抑制后的結果。由于不平整表面的影響，MS和SVD呈現相似的結果，雜波抑制后的目標圖像仍包含一些雜波部分，但去除了2個最大奇異值的SVD降低了目標信號的強度，IF值低于MS。RPCA方法將目標與雜波分離良好，基于RDAE的雜波抑制方法可以更有效地從原始B-scan圖像中提取目標響應，并且在4種雜波抑制方法中具有最清晰的背景、最高的IF值。

圖9 不同方法對真實數據的雜波抑制結果Fig.9 Clutter suppression results for real data with different methods

2.3 相對介電常數估計

為了評估使用直接最小二乘法估計介質相對介電常數的準確性，在gprMax中模擬了12種不同的均質土壤，相對介電常數從2(干燥土壤)到13(潮濕土壤)。通過計算均方根相對誤差(Root-mean-square relative error, RMSRE)來評估準確性。

按照圖2所示流程，土壤相對介電常數的估計結果如圖10所示，擬合真實值和估計值得到擬合曲線為y=0.965x-0.018，均方根相對誤差為3.84%，相對介電常數估計有較高的準確性。

圖10 用DLS估算土壤相對介電常數的結果Fig.10 Results of soil relative permittivity estimation using DLS method

2.4 根系定位

圖11 均質土壤單根偏移成像結果Fig.11 Migration imaging results of single tree-root in homogeneous soil

對圖7中的分層模型進行偏移成像，在使用DLS估計模型相對介電常數時，為方便計算將其視為均一介質計算模型的等效介電常數[31]，得到整個模型相對介電常數為7.29，偏移速度為0.111 m/ns，將其作為F-K偏移的輸入，得到圖12a所示偏移成像結果，圖像中雙曲線未完全消除，偏移速度小于實際速度，點狀目標未能很好地聚焦。當速度增大為0.117 m/ns時，偏移成像如圖12b所示，能很好地聚焦于一點。

圖12 分層模型單根偏移成像結果Fig.12 Migration results of single tree-root in stratified model

圖13a為2.2節中實測模型的實物圖，將半徑0.052 m的香樟樹枝R3埋在0.3 m深的沙子中，對雜波抑制后的圖像(圖9e)使用DLS估算出的沙子相對介電常數為12.46，偏移速度為0.085 m/ns，圖13b為偏移成像結果。在相同的沙子介質中，依次放置R4、R5、R6這3個半徑不同的香樟樹枝于0.15、0.25、0.35 m深處，如圖13c所示。圖13d為雜波抑制后的偏移成像結果。

圖13 實測數據偏移成像結果Fig.13 Migration results of real data

2.5 根系定位結果及分析

表4 偏移前后根系參數對比Tab.4 Comparison of tree-root parameters before and after migration

圖14 偏移前后根系位置和尺寸對比Fig.14 Comparison of tree-root position and size before and after migration

表3 偏移前后根系參數對比Tab.3 Comparison of tree-root parameters before and after migration m

本文利用探地雷達檢測儀在檢測區域內掃描根系，收集數據，相較于根鉆法和全根挖掘法等傳統破壞性檢測方法，檢測快速，勞動成本低并且不會對周圍根系造成不可逆轉的破壞。為了進一步說明本文方法的優勢，以文中仿真城市路面結構的B-scan圖像和根R3的實測B-scan圖像為樣本，針對雜波抑制、雙曲線擬合及偏移成像3個步驟，分別對比了與RPCA、HT和Kirchhoff偏移的計算耗時。仿真的B-scan圖像尺寸為3 393×180，即由180個A-scan跡線組成，每個A-scan跡線包含3 393個采樣點，實測的B-scan圖像尺寸為313×409，表5為仿真數據和實測數據在根系定位3個步驟與其它方法的計算耗時對比，雖然RDAE對單個B-scan圖像耗時大于RPCA，但是RDAE方法具有更高的信雜比和IF值，并且可以使用預訓練的深度自動編碼器降低對多個B-scan圖像的雜波抑制時間。DLS和F-K偏移分別在雙曲線擬合和偏移成像步驟中計算耗時低于HT和Kirchhoff偏移，總體來說，本文所提的方法在計算耗時上優于其它3種方法。

表5 計算耗時對比Tab.5 Comparison of computational time s

3 結論

(1)將探地雷達B-scan圖像建模成表示雜波的低秩矩陣和表示目標回波的稀疏矩陣，雜波抑制問題轉換為低秩和稀疏矩陣分解問題。實驗表明，本文提出的RDAE雜波抑制方法在仿真和實測數據的雜波抑制效果均優于RPCA等傳統方法，雜波抑制后的圖像擁有更清晰的背景和更高的IF值。

(2)DLS擬合雙曲線能更精確地估計土壤的相對介電常數，在數值模擬實驗中估算的土壤相對介電常數的均方根相對誤差為3.84%，為F-K偏移提供了可靠的輸入。

(3)對于根系定位問題，本文提出基于RDAE雜波抑制的F-K偏移根系定位方法，在仿真和實測數據上都有較高的定位精度，偏移位置與真實位置的最大圓心距離為3.5 cm，最大半徑相對誤差和最大深度相對誤差分別為8.5%、8.7%，能夠應對復雜的分層介質模型，滿足實際應用中樹木根系無損檢測的需求。