基于GprMax 的地下冬筍正演模擬研究

余 鋮，周凱琪，于 航，溫春雨，樊夢成，楊科成

（浙江農(nóng)林大學(xué)光機(jī)電工程學(xué)院，浙江 杭州 311300）

0 引言

我國是世界上竹林資源最豐富的國家[1-2]，2022年竹林面積達(dá)701萬hm2，冬筍[3]是竹林收益的主要來源之一。由于冬筍埋藏于地下，目前多依賴于農(nóng)民的經(jīng)驗采收，效率低下，盲目挖掘也會破壞植物根莖。

探地雷達(dá)[4]（Ground Penetrating Radar, GPR）是指利用收發(fā)天線發(fā)出并接收高頻電磁波，來探測地下物質(zhì)分布規(guī)律的一種無損檢測方式，具有精度高、效率高、檢測特性平穩(wěn)、對環(huán)境適應(yīng)性好等優(yōu)勢，被廣泛應(yīng)用于植物根系探測領(lǐng)域。利用探地雷達(dá)探測冬筍是一種比較好的方法，但是探地雷達(dá)的缺陷就在于，要求研究人員必須具備解譯其信息的知識，也因為探地雷達(dá)信息的復(fù)雜性，解譯探地雷達(dá)信息時往往要求具有豐富的專業(yè)知識和經(jīng)驗[5]。尤其是地下環(huán)境復(fù)雜，使探地雷達(dá)信息的復(fù)雜度更高。

通過正演模擬可以高效積累雷達(dá)圖像解譯的經(jīng)驗，減少實測成本。本文基于GprMax[6-7]仿真軟件，以地下冬筍為對象，研究了不同干擾物、不同大小冬筍與不同深度冬筍的探地雷達(dá)正演模型，總結(jié)了冬筍的雷達(dá)探測經(jīng)驗。

1 探地雷達(dá)及正演原理

1.1 探地雷達(dá)原理介紹

探地雷達(dá)主要由雷達(dá)主機(jī)、顯示單元、接收發(fā)天線組成，是利用10 MHz 以上的高頻電磁波束對地下介質(zhì)成像的無損探測技術(shù)。探地雷達(dá)的工作原理如圖1 所示，工作時通過發(fā)射天線（TX）向地下發(fā)射高頻電磁波，當(dāng)遇到具有電性差異的界面時發(fā)生反射，由接收天線（RX）收集反射波，并經(jīng)采樣電路采集信號傳輸至雷達(dá)主機(jī)，計算出反射波的相位信息，擬合成探地雷達(dá)的二維剖面圖。

圖1 探地雷達(dá)工作原理

1.2 正演仿真原理及GprMax

在宏觀尺度上，所有的電磁現(xiàn)象都可由麥克斯韋方程組描述。以下為麥克斯韋方程組描述基本電磁場能量與其對源的影響情況的一階偏微分方程：

式中，E、H、B、D分別是電場強(qiáng)度、磁場強(qiáng)度、磁場感應(yīng)強(qiáng)度和電位移矢量；J是電流密度；q是電荷密度；t為時間。從模型的幾何形式和初始條件出發(fā)，通過利用麥克斯韋方程組，可以進(jìn)行探地雷達(dá)數(shù)據(jù)模擬，但是在模擬仿真過程中，麥克斯韋方程連續(xù)性的特點，使得它在計算機(jī)中永遠(yuǎn)無法求解。

Y e e[8]在1 9 9 6 年提出了時域有限差分法（FDTD），通過利用離散空間與時間的連續(xù)性，并利用在時域的遞推模擬了波的擴(kuò)散行為，由此得到了場分布圖。Yee 網(wǎng)絡(luò)如圖2 所示，是離散型FDTD 網(wǎng)格的代表形式。

圖2 FDTD Yee網(wǎng)絡(luò)

此外，在進(jìn)行電磁場仿真模擬時，若在空間無限遠(yuǎn)時將場傳遞至零值，則仍不能利用有限的計算空間來求解無限的仿真問題。因此，在距離發(fā)射源足夠遠(yuǎn)的地方應(yīng)用完全匹配層（PML）邊界吸收條件對電磁場計算進(jìn)行截斷，從而限制計算空間，邊界吸收條件通過吸收傳播過來的場，模擬出一個無邊界的空間。

GprMax 是一個可以進(jìn)行電磁波傳輸模擬的開源軟件，在采用了時域有限差分法（FDTD）和PML[9]（Perfectly Matched Layer）的邊界吸收條件下，以Yee元胞為電磁場的離散單元，將Maxwell 旋度微分方程轉(zhuǎn)化為差分方程，以達(dá)到正演模擬的效果。因在一定時間上迭代求解，具有很高的可靠性和收斂性。GprMax 設(shè)計之初主要用作探地雷達(dá)（GPR）的模型，同時也可能用作一些其他的應(yīng)用模型。GprMax 還具有內(nèi)置的天線模型庫，可以進(jìn)行粗糙地基的模擬搭建、復(fù)雜各向異性結(jié)構(gòu)的模擬等，可以適應(yīng)大部分環(huán)境的搭建要求。

2 實驗研究與討論

冬筍通常和竹子的根系一起生長，周圍難免會有石頭等介質(zhì)的干擾，使用探地雷達(dá)探測地下冬筍時，需要分辨出這些干擾信息。為了在不破壞生態(tài)的情況下，挖掘出合適生長時期的冬筍，還需判斷出冬筍的埋藏深度與大小。

本實驗通過GprMax 模擬了不同形狀的干擾物，包括圓形截面的根系、矩形截面的石頭，同時模擬了不同深度、不同大小的冬筍的正演模擬圖像。模型參數(shù)如表1 所示[10]。

表1 模型參數(shù)

2.1 不同形狀介質(zhì)正演模擬

首先，探究不同形狀、相同埋藏深度的模型的探地雷達(dá)正演模擬，設(shè)計模型整體寬1 m，深0.5 m。由于在GprMax 的運(yùn)行中，發(fā)射天線和接收天線都需要占據(jù)模型中一定的空間，為了給發(fā)射天線和接收天線預(yù)留位置，介質(zhì)層的厚度設(shè)置為0.1 m，空間步長設(shè)置為0.1 m。探測深度在0.05 m～0.15 m之間，使用Ricker 子波源模擬Gpr 的天線，激發(fā)頻率選擇中心頻率為1 500 MHz。

首先對不同介質(zhì)模型進(jìn)行正演模擬，如圖3 所示，模型a 為圓形截面的根系模型，模型b 為矩形截面的巖石模型，模型c為冬筍模型。

圖3 不同介質(zhì)GprMax模型構(gòu)建示意圖

不同形狀介質(zhì)的GprMax 正演模擬結(jié)果如圖4 所示。從形狀類別上來看，可以根據(jù)圖像的雙曲線特征來判斷介質(zhì)的形狀，圓形介質(zhì)的B 掃描圖像平滑、規(guī)整且沒有交疊的現(xiàn)象；而矩形介質(zhì)的B 掃描圖像頂端出現(xiàn)了交疊和顏色加深現(xiàn)象，在交點處有繞射；冬筍介質(zhì)的內(nèi)部有多次反射，且內(nèi)部的反射波雜亂無章。以此特征可作為判斷冬筍是否存在的依據(jù)。

圖4 不同形狀介質(zhì)的GprMax正演模擬結(jié)果

2.2 不同大小冬筍正演模擬

與上述模型相同，構(gòu)建土壤與空氣層模型，然后設(shè)置埋藏深度相同、不同大小的冬筍，如圖5 所示，從左至右分別為長度15 cm、10 cm 與5 cm 的冬筍回波示意圖，埋藏深度均為5 cm。

圖5 不同大小冬筍GprMax模型構(gòu)建示意圖

不同大小冬筍的正演模擬結(jié)果如圖6 所示，從圖中可以看出，冬筍的上下界面反射出多道回波，冬筍的長度越長，則上下界面的雙曲線間隔距離越大。因此，可以通過兩條雙曲線出現(xiàn)的時間間隔估算冬筍的縱向尺寸。

圖6 不同大小冬筍的正演模擬結(jié)果

2.3 不同深度冬筍正演模擬

最后構(gòu)建大小相同、不同埋藏深度的冬筍模型，如圖7 所示，從左至右分別為埋藏深度5 cm、10 cm、15 cm的冬筍，冬筍長度均設(shè)置為15 cm。

圖7 不同埋藏深度冬筍GprMax模型構(gòu)建示意圖

相同長度、不同埋藏深度冬筍的GprMax 正演模擬結(jié)果如圖8 所示，從回波數(shù)據(jù)來看，雙曲線的頂點位置反映出冬筍的埋藏深度，通過分析此特征可以判斷出介質(zhì)的埋藏深度。冬筍的埋藏深度越大，雙曲線的出現(xiàn)時間則越晚，隨著深度增加，雙曲線的開口大小也有所增加。因此，埋藏越深的冬筍，雙曲線出現(xiàn)時間越晚、開口越大，埋藏越淺的冬筍則相反。

圖8 不同埋藏深度冬筍的GprMax正演模擬結(jié)果

3 結(jié)束語

為給探地雷達(dá)在冬筍探測領(lǐng)域的應(yīng)用提供經(jīng)驗，本文利用GprMax 軟件對不同形狀干擾介質(zhì)、不同大小及不同埋藏深度的冬筍等模型進(jìn)行正演模擬，并進(jìn)行了模型間的對比。結(jié)果表明，正演模擬的圖像能較好地反映各不同模型的成像特征，能有效反饋出介質(zhì)模型的物理特性如大小、形狀、深度等信息。截面形狀不同的物體的雙曲線回波形狀不同，雙曲線的上下界面距離反映出物體的大小，頂點位置則反映出介質(zhì)的埋藏深度。通過分析這些特征可以有效探測識別出地下冬筍的存在、生長時期（大?。┡c埋藏深度。本文所應(yīng)用的GprMax 正演模擬方法有效揭示了探地雷達(dá)的圖像規(guī)律，也證實了探地雷達(dá)在冬筍探測領(lǐng)域的可行性。