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

李光輝 徐 匯 劉 敏

(江南大學(xué)人工智能與計(jì)算機(jī)學(xué)院， 無錫 214122)

0 引言

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

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

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

1 研究方法

1.1 魯棒深度自動(dòng)編碼器

從數(shù)學(xué)上講，觀測到的探地雷達(dá)B-scan數(shù)據(jù)X可以由一個(gè)表示直達(dá)波等背景信號(hào)的低秩矩陣L，一個(gè)表示目標(biāo)回波信號(hào)的稀疏矩陣S和一個(gè)隨機(jī)噪聲矩陣N聯(lián)合表示[26]，即

X=L+S+N

(1)

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

(2)

式中 ‖·‖*——核范數(shù)

‖·‖1——l1范數(shù)

‖·‖F(xiàn)——F范數(shù)

λ——控制稀疏矩陣稀疏性的懲罰系數(shù)

ε——任意正數(shù)

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

圖1 自動(dòng)編碼器結(jié)構(gòu)圖Fig.1 Architecture of autoencoder

(3)

式中 ‖·‖2——l2范數(shù)

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

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

(4)

1.2 直接最小二乘法

一般地，當(dāng)探地雷達(dá)天線掃描經(jīng)過樹根上部土壤時(shí)，樹根反射雷達(dá)脈沖形成的回波雙曲線可以表示為

(5)

式中x——天線位置

t——雷達(dá)波的雙向傳播時(shí)間

c——真空中光速

εr——土壤的相對介電常數(shù)

(x0,t0)——雙曲線頂點(diǎn)坐標(biāo)

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

圖2 基于DLS的雙曲線參數(shù)估計(jì)算法流程圖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偏移

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

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

設(shè)探地雷達(dá)B-scan記錄剖面為p(x,z=0,t)，z是垂直坐標(biāo)，向下為正，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)在頻率-波數(shù)域做波場外推，深度z處的波場可以表示為

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

(7)

由于探地雷達(dá)記錄電磁波的雙程走時(shí)，因此可以把波速視為介質(zhì)中波速v的一半，根據(jù)色散方程，kx、kz和ω的關(guān)系可以表示為

(8)

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

(9)

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

圖5 改進(jìn)的F-K偏移流程圖Fig.5 Flow chart of improved F-K migration

2 實(shí)驗(yàn)與結(jié)果分析

2.1 實(shí)驗(yàn)數(shù)據(jù)采集

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

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

2.2 雜波抑制

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

圖7 城市路面結(jié)構(gòu)的仿真模型Fig.7 Simulation model of urban road structure

表1 仿真基本參數(shù)Tab.1 Simulation parameters

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

(10)

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

(11)

式中Ns——目標(biāo)信號(hào)區(qū)域Rs中的像素?cái)?shù)

Nc——雜波區(qū)域Rc中的像素?cái)?shù)

I(p)——B-scan圖像第p個(gè)像素的信號(hào)強(qiáng)度

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

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

表2 仿真和實(shí)測數(shù)據(jù)的IF結(jié)果Tab.2 IF results of simulated and real data

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

圖9 不同方法對真實(shí)數(shù)據(jù)的雜波抑制結(jié)果Fig.9 Clutter suppression results for real data with different methods

2.3 相對介電常數(shù)估計(jì)

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

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

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

2.4 根系定位

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

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

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

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

圖13 實(shí)測數(shù)據(jù)偏移成像結(jié)果Fig.13 Migration results of real data

2.5 根系定位結(jié)果及分析

表4 偏移前后根系參數(shù)對比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 偏移前后根系參數(shù)對比Tab.3 Comparison of tree-root parameters before and after migration m

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

表5 計(jì)算耗時(shí)對比Tab.5 Comparison of computational time s

3 結(jié)論

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

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

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