簡介gprMax是一套用于探地雷達(dá)數(shù)值建模的開源軟件基于有限差分時(shí)域FDTD方法在三維空間內(nèi)求解麥克斯韋方程面向地球物理勘探、道路檢測、建筑結(jié)構(gòu)評(píng)估及天線設(shè)計(jì)等領(lǐng)域的科研人員和工程技術(shù)人群。壓縮包整體約33.1MB共收錄373個(gè)文件以Python源碼.py與Cython加速模塊.pyx為核心配合豐富的仿真輸入文件.in、結(jié)果輸出數(shù)據(jù).out、.npz、.vti、文檔說明.rst、.md、.pdf及圖片.png目錄模塊劃分清晰便于按需檢索。該資源已有3154人學(xué)習(xí)內(nèi)容涵蓋CPU并行求解器、基于CUDA的GPU求解器、GSSI與MALA等典型天線模型、后處理腳本以及入門示例讀者可通過源碼級(jí)代碼快速搭建仿真環(huán)境完整運(yùn)行典型GPR探測案例理解FDTD建模與計(jì)算流程并為后續(xù)算法改進(jìn)或二次開發(fā)提供可復(fù)用的參考實(shí)現(xiàn)。 做探地雷達(dá)項(xiàng)目那陣子我?guī)缀趺刻於家獙?duì)著儀器采集回來的B-scan發(fā)愁地下那個(gè)目標(biāo)到底長什么樣、埋了多深、介電常數(shù)差多少才會(huì)出現(xiàn)這種雙曲線繞射直接在實(shí)驗(yàn)場地上挖坑驗(yàn)證成本太高換場地又費(fèi)時(shí)間。后來我把目光轉(zhuǎn)到數(shù)值模擬上用gprMax這款開源軟件做正演把“地下有什么”變成可控變量一次能生成幾十組模擬數(shù)據(jù)再回頭跟實(shí)測波形對(duì)比很多疑問就清晰了。gprMax基于有限差分時(shí)域FDTD方法求解麥克斯韋方程組模擬電磁波在介質(zhì)中的傳播過程。軟件完全開源項(xiàng)目托管在GitHub核心用Python編寫數(shù)值內(nèi)核用Cython加速。無論是研究GPR天線特性、評(píng)估不同地下目標(biāo)的回波特征還是為反演算法準(zhǔn)備訓(xùn)練數(shù)據(jù)它都是目前社區(qū)里最常用的GPR數(shù)值建模工具之一。這篇文章寫給準(zhǔn)備用gprMax做建模、但不想把官方文檔從頭啃到尾的人我會(huì)從原理講到能跑的模型最后把幾個(gè)常見的坑也一并說掉。1. 先把問題講清楚GPR模擬到底在做什么1.1 數(shù)值建模要解決的實(shí)際問題探地雷達(dá)的觀測方式并不復(fù)雜發(fā)射天線向地下輻射高頻電磁波電磁波在介質(zhì)分界面處發(fā)生反射接收天線記錄回波。但真正做數(shù)據(jù)分析時(shí)麻煩就來了。同樣一個(gè)空洞在混凝土里和土層里回波幅度和相位特征完全不同同樣一段鋼筋埋深不同雙曲線繞射的開口大小也不一樣。實(shí)測數(shù)據(jù)沒法告訴你“這個(gè)波形到底對(duì)應(yīng)什么物理結(jié)構(gòu)”因?yàn)榈叵陆橘|(zhì)的組合幾乎是無限的。數(shù)值建模的價(jià)值就在這里。你可以把目標(biāo)埋深、尺寸、填充介質(zhì)、背景介電常數(shù)全部設(shè)成已知參數(shù)跑一次正演模擬得到理論上“應(yīng)該出現(xiàn)”的波形。然后把實(shí)測波形和模擬波形對(duì)比反復(fù)調(diào)整模型參數(shù)就能反推地下的真實(shí)情況。這個(gè)思路用在管線探測、道路病害檢測、橋墩無損評(píng)估上都成立也能幫雷達(dá)設(shè)備廠商做天線設(shè)計(jì)和指標(biāo)驗(yàn)證。1.2 為什么偏偏是FDTD電磁場數(shù)值方法有好幾種矩量法MoM、有限元FEM、時(shí)域有限差分FDTD各有各的擅長場景。gprMax選擇FDTD我個(gè)人的理解是它和GPR這類寬帶、有耗、非均勻介質(zhì)問題特別搭。FDTD直接在時(shí)間域離散麥克斯韋方程組一次計(jì)算就能得到寬頻帶響應(yīng)而GPR用的恰恰是ns級(jí)脈沖頻譜覆蓋從幾十MHz到幾GHz用頻域方法需要逐頻點(diǎn)求解計(jì)算量會(huì)翻好幾倍。另一個(gè)優(yōu)點(diǎn)是FDTD處理非均勻介質(zhì)很方便混凝土、土壤、鋼筋、空洞這種非規(guī)則分布的結(jié)構(gòu)只需給每個(gè)網(wǎng)格點(diǎn)賦予不同電磁參數(shù)即可不需要像有限元那樣重新劃分網(wǎng)格。代價(jià)是它需要滿足CFL穩(wěn)定性條件時(shí)間步長不能太大而且三維模型網(wǎng)格多了內(nèi)存會(huì)緊張這個(gè)后面細(xì)說。1.3 gprMax到底能輸出什么結(jié)果gprMax輸出的物理量主要是接收點(diǎn)處的電場分量對(duì)應(yīng)實(shí)測GPR的A-scan單道波形。把多個(gè)接收點(diǎn)按測線排列把A-scan拼起來就能得到B-scan二維剖面圖也就是我們常說的雷達(dá)圖像這是GPR數(shù)據(jù)分析最?？吹母袷?。三維情況下還能生成C-scan也就是某個(gè)深度切片上的回波強(qiáng)度分布。此外gprMax還可以輸出幾何模型視圖geometry view用來檢查建模時(shí)設(shè)置的介質(zhì)分布是否符合預(yù)期。這點(diǎn)非常重要模型里物體位置形狀設(shè)置錯(cuò)了跑出來的波形再漂亮也是廢的。2. 環(huán)境準(zhǔn)備環(huán)境裝好后面才順2.1 獲取gprMax的幾種方式gprMax的開發(fā)維護(hù)主要靠GitHub倉庫最新版本通常在那里最先發(fā)布。我一般推薦直接克隆源碼倉庫git clone https://github.com/gprMax/gprMax.git cd gprMax python setup.py install如果你不想動(dòng)源碼也可以用pip直接裝pip install gprMax兩種方式我都試過。pip方式簡單但有時(shí)版本更新不及時(shí)源碼方式能直接看到底層代碼后續(xù)如果要做二次開發(fā)或調(diào)試建議用源碼方式。另外gprMax的官方文檔網(wǎng)站提供了詳細(xì)的API說明和示例模型建議下載源碼時(shí)把examples目錄一并保留里面的tutorials系列是很好的入門素材。2.2 用Anaconda搭環(huán)境順便換個(gè)國內(nèi)源gprMax是Python程序底層依賴numpy、cython、h5py這些庫。我習(xí)慣先裝Anaconda再為gprMax單獨(dú)建一個(gè)虛擬環(huán)境避免跟其他項(xiàng)目的Python包沖突conda create -n gprmax_env python3.9 conda activate gprmax_env pip install gprMax國內(nèi)網(wǎng)絡(luò)環(huán)境下直接從官方源下載包經(jīng)常很慢。清華開源軟件鏡像站是我一直在用的方案pip源和conda源都能加速換源命令也不復(fù)雜pip config set global.index-url https://pypi.tuna.tsinghua.edu.cn/simple裝好之后驗(yàn)證一下能不能正常導(dǎo)入python -c import gprMax; print(gprMax.__version__)如果沒報(bào)錯(cuò)說明環(huán)境基本就緒。老版本的gprMax還需要額外安裝Cython并編譯內(nèi)核新版在pip install時(shí)會(huì)自動(dòng)處理這一點(diǎn)省了不少事。2.3 快速跑通一個(gè)官方示例環(huán)境配好后別著急自己寫模型先把官方示例跑一遍。gprMax的examples目錄下有很多. in文件找一個(gè)最簡單的跑python -m gprMax user_models/Bscan_2D.in -n 1這會(huì)在同目錄下生成Bscan_2D.out等多個(gè)文件。終端會(huì)打印模擬進(jìn)度和耗時(shí)。第一次跑通這個(gè)示例說明安裝沒問題也讓你對(duì)軟件的使用流程有個(gè)直觀感受。3. 第一次建模仿真混凝土空洞檢測3.1 模型文件從哪開始寫gprMax的模型文件是.in后綴的文本文件格式相對(duì)簡單。核心思路就是把你要模擬的場景用“定義域、網(wǎng)格、材料、波源、接收點(diǎn)、目標(biāo)幾何體”這幾個(gè)要素描述出來。我寫一個(gè)實(shí)際用過的模型模擬混凝土板內(nèi)部存在空氣空洞的場景。設(shè)模型為2D板厚0.3米長1.2米發(fā)射天線放在表面測線沿x方向布置。#title: concrete_void_2d #domain: 1.2 0.3 0.02 #dx_dy_dz: 0.002 0.002 0.002 #time_window: 1.5e-8 #material: 6 0.005 1 0 concrete #material: 1 0 1 0 air #waveform: ricker 1 1.5e9 rick #hertzian_dipole: z 0.10 0.15 0.01 rick #rx: 0.30 0.15 0.01 #rx: 0.50 0.15 0.01 #rx: 0.70 0.15 0.01 #box: 0.55 0.10 0.005 0.65 0.16 0.015 air #geometry_view: 0 0 0 1.2 0.3 0.02 0.002 0.002 0.002 void_geo geo簡單解釋一下。#domain定義模型尺寸單位是米#dx_dy_dz是空間步長決定網(wǎng)格細(xì)度#time_window是模擬時(shí)長窗口#material定義材料參數(shù)第一個(gè)數(shù)字是相對(duì)介電常數(shù)第二個(gè)是電導(dǎo)率#hertzian_dipole表示偶極子源理想化的點(diǎn)源z方向極化#rx是接收點(diǎn)坐標(biāo)#box在指定區(qū)域內(nèi)填充材料這里把一大塊混凝土區(qū)域變成空氣模擬空洞#geometry_view則用來導(dǎo)出幾何模型圖。3.2 材料與波源參數(shù)怎么定材料參數(shù)是GPR模擬里最容易隨意設(shè)置、也最容易出錯(cuò)的地方?；炷恋南鄬?duì)介電常數(shù)通常取6到9之間干燥混凝土偏低濕混凝土偏高。這里取6.0電導(dǎo)率取0.005 S/m屬于常見經(jīng)驗(yàn)值?？諝獾慕殡姵?shù)是1電導(dǎo)率是0直接作為材料定義即可。波源方面gprMax內(nèi)置了多種波形最常用的是Ricker子波。#waveform: ricker 1 1.5e9 rick表示用Ricker子波振幅為1中心頻率1.5GHz。這個(gè)頻率在探地雷達(dá)里屬于中等偏高的頻段對(duì)厘米級(jí)空洞的分辨率較好但穿透深度有限。如果目標(biāo)埋深超過1米我通常會(huì)把頻率降到500MHz甚至更低不然高頻能量衰減太快。3.3 跑起來并讀取結(jié)果在模型文件所在目錄執(zhí)行python -m gprMax concrete_void_2d.in -n 1-n參數(shù)是并行線程數(shù)單核模擬時(shí)間可能較長多核能明顯加快速度。模擬結(jié)束后會(huì)生成concrete_void_2d.out文件新版gprMax的輸出文件是HDF5格式可以用h5py直接讀import h5py import matplotlib.pyplot as plt f h5py.File(concrete_void_2d.out, r) rx1_data f[/rxs/rx1/Ez][:] plt.plot(rx1_data) plt.show()讀出來的數(shù)據(jù)就是接收點(diǎn)處的電場時(shí)域波形也就是A-scan。多個(gè)接收點(diǎn)的數(shù)據(jù)拼在一起就能得到B-scan剖面圖。我第一次看到空洞位置在B-scan上出現(xiàn)明顯的雙曲線繞射時(shí)算是真正理解了GPR數(shù)據(jù)形態(tài)是怎么來的。4. 提升精度與性能空間步長、時(shí)間窗與內(nèi)存4.1 空間步長的“每波長10個(gè)網(wǎng)格”法則網(wǎng)格步長是FDTD模擬里最關(guān)鍵的參數(shù)之一。步距太大數(shù)值色散嚴(yán)重波前會(huì)變形步距太小網(wǎng)格數(shù)暴增內(nèi)存和計(jì)算時(shí)間都頂不住。業(yè)內(nèi)常用經(jīng)驗(yàn)是每個(gè)最小波長內(nèi)至少要有10個(gè)網(wǎng)格。介質(zhì)中的波長計(jì)算公式是lambda c / (f * sqrt(eps_r))以1.5GHz、混凝土介電常數(shù)6為例lambda 3e8 / (1.5e9 * sqrt(6)) ≈ 0.082米如果按10個(gè)網(wǎng)格算網(wǎng)格步長取0.008米就夠了。我在上面模型里用的是0.002米每波長大約40個(gè)網(wǎng)格精度有富余代價(jià)是計(jì)算時(shí)間變長。實(shí)際項(xiàng)目中可以根據(jù)目標(biāo)尺寸和計(jì)算資源折中我一般控制在每波長15到20個(gè)網(wǎng)格之間。4.2 時(shí)間窗踩到兩個(gè)極限#time_window設(shè)的是模擬總時(shí)長。設(shè)太短目標(biāo)回波還沒走完邊界就截?cái)嗔嗽O(shè)太長白白增加時(shí)間步迭代次數(shù)。該怎么判斷設(shè)多少合適簡單估算方法電磁波從發(fā)射點(diǎn)到目標(biāo)再返回接收點(diǎn)的總路徑長度除以介質(zhì)中的波速再留一點(diǎn)余量。比如目標(biāo)深度0.15米發(fā)射接收都在表面往返路徑約0.3米混凝土中波速約1.22e8 m/s那么回波到達(dá)時(shí)間大約是2.5ns。我把時(shí)間窗設(shè)成15ns一方面把多次反射也包含進(jìn)去另一方面也為地下更深處的微弱回波留足空間。還要注意FDTD的內(nèi)穩(wěn)性。gprMax會(huì)根據(jù)空間步長自動(dòng)計(jì)算時(shí)間步長確保滿足CFL條件但用戶設(shè)置的時(shí)間窗如果過大會(huì)顯著增加迭代次數(shù)。我的經(jīng)驗(yàn)是先跑一個(gè)短時(shí)間窗的快速測試確認(rèn)波形合理后再拉長時(shí)間窗做精細(xì)模擬。4.3 三維大模型的內(nèi)存優(yōu)化思路很多人一上來就建三維模型跑著跑著內(nèi)存就爆了。一個(gè)三維模型的總網(wǎng)格數(shù)是三個(gè)方向網(wǎng)格數(shù)的乘積網(wǎng)格密度稍微高一點(diǎn)網(wǎng)格數(shù)輕松破億。我通常這樣控制內(nèi)存。第一能降維就降維。很多問題在二維模型下就能說明白精度足夠速度卻快兩個(gè)數(shù)量級(jí)。第二用對(duì)稱性。如果模型結(jié)構(gòu)和波源位置關(guān)于某個(gè)平面對(duì)稱可以只建一半配合對(duì)稱邊界條件。第三適當(dāng)增大網(wǎng)格步長犧牲少量精度換取內(nèi)存空間。第四如果條件允許用多線程模式跑現(xiàn)代處理器的多核優(yōu)勢能發(fā)揮出來。5. 我踩過的坑波形發(fā)散、直接波掩蓋、噪聲異常5.1 波形發(fā)散最常見的原因模擬出來的波形隨時(shí)間增大而振幅爆炸式增長這種情況我遇到過好多次。排查下來大部分原因是材料參數(shù)不合理。具體來說某層介質(zhì)的電導(dǎo)率如果設(shè)成一個(gè)異常大的值比如10以上FDTD迭代時(shí)場值可能出現(xiàn)指數(shù)增長。另外如果兩種相鄰材料的介電常數(shù)差異巨大而網(wǎng)格又太粗分界面處會(huì)產(chǎn)生非物理反射看起來像發(fā)散。遇到這類問題第一步先檢查材料參數(shù)是否在合理范圍內(nèi)第二步把網(wǎng)格加密一檔再看波形是否收斂。5.2 直接波太強(qiáng)目標(biāo)回波被淹沒在GPR實(shí)測中發(fā)射天線和接收天線之間的直達(dá)波、地表反射波往往幅度非常大目標(biāo)回波反而很弱。模擬中也一樣靠近源位置的接收點(diǎn)數(shù)據(jù)目標(biāo)回波可能完全被直接波壓住。處理辦法有兩種。一種是把發(fā)射天線和接收天線拉開一定距離模擬分離式天線的場景減少直耦影響。另一種是跑兩次模擬第一次不放目標(biāo)物體得到“背景場”第二次放目標(biāo)得到“總場”兩者相減就能把目標(biāo)散射場單獨(dú)提取出來。這個(gè)思路也可以用在實(shí)測數(shù)據(jù)的背景扣除處理上。5.3 輸出數(shù)據(jù)有異常噪聲先看幾何視圖有幾次我跑出的B-scan上出現(xiàn)了完全不符合物理直覺的強(qiáng)反射條帶檢查模型文件參數(shù)也覺得沒問題但問題就出在目標(biāo)位置和邊界太近吸收入射波的PML沒有完全發(fā)揮作用。這時(shí)候一定要先看#geometry_view導(dǎo)出的幾何模型圖確認(rèn)目標(biāo)、邊界、材料分布和預(yù)想一致。gprMax的幾何視圖是灰色的介質(zhì)灰度圖不同介質(zhì)用不同灰度區(qū)分很容易看出設(shè)置錯(cuò)誤。另外一個(gè)容易忽略的點(diǎn)是接收點(diǎn)的位置。如果接收點(diǎn)正好落在網(wǎng)格奇點(diǎn)上或離波源太近讀數(shù)可能異常。我的習(xí)慣是避開整數(shù)網(wǎng)格位置接收點(diǎn)坐標(biāo)稍微移動(dòng)半個(gè)網(wǎng)格結(jié)果會(huì)更穩(wěn)定。最后分享一點(diǎn)個(gè)人體會(huì)如果你剛開始用gprMax我的建議是先別急著改自己的模型參數(shù)把官方tutorials里的幾個(gè)例子完整跑一遍挨個(gè)打開輸出的A-scan和B-scan看看。我當(dāng)初就是跳過了這一步直接上手改復(fù)雜模型結(jié)果花了大量時(shí)間排查一個(gè)本來很簡單的語法錯(cuò)誤。等到你把“建模型、跑模擬、看結(jié)果”這個(gè)閉環(huán)跑順了再回頭去調(diào)材料參數(shù)、做參數(shù)掃描、加并行計(jì)算就都順理成章了。數(shù)值模擬這個(gè)東西前期慢一點(diǎn)后面反而快很多。本文還有配套的精品資源點(diǎn)擊獲取