1、GprMax是愛丁堡大學的Antonis Giannopoulos于1996年推出來的一種基于時域有限差分（FDTD）算法和理想匹配層（PML）邊界吸收條件的探地雷達正演數值模擬軟件，用于探地雷達成像研究。其中，GprMax2D是二維正演，GprMax3D為三維正演。該軟件可以在Windows、Linux和MacOS三個平臺上使用。本文主要針對Windows平臺進行說明。一、軟件獲得該軟件為免費軟件，可以去GprMax官網下載。也可點此直接下載。二、準備工作軟件無需安裝，下載后用解壓工具解壓，找到Windows文件夾，直接雙擊GprMax2D.exe即可運行。但是為了方便批量模擬，建議把Win

2、dows文件夾下的文件（cygwin1.dll、GprMax2D.exe和GprMax3D.exe）全部復制到系統盤系統搜索路徑下，例如：C:Windows。如果只進行二維正演，只復制cygwin1.dll和GprMax2D.exe即可。當然也可放在任意路徑下，只要使用時包含所在路徑就OK了。筆者習慣放在系統要目錄下，即C:。另外，最好建立一個輸入文件和模擬結果存放的專用文件夾，并且把tools文件夾下的文件（gprmax.m、gprmax2g.m、gprmax3g.m、gprmaxde.m和gprmaxso.m）全部復制到該專用文件夾下，注意文件夾名最好使用英文。如果只進行二維正演，只復制g

3、prmax.m和gprmax2g.m即可。筆者習慣使用D:GPR。正演結果需要用MATLAB進行繪圖，因此需要安裝有MATLAB軟件。軟件下載地址和安裝方法此處不再說明，以后使用將假設讀者已經成功安裝MATLAB軟件。三、軟件使用一般進行數值模擬時通常都是若干個對比模型進行模擬，因此本文只介紹批量模擬的方法。如此一來，筆者將認為讀者的軟件和筆者一樣放在C:，輸入文件和gprmax.m、gprmax2g.m文件已經放在D:GPR。欲進行正演模擬需要先建立輸入文件，輸入文件的建立將在后面介紹，此處假設在D:GPR已經有若干輸入文件，文件名分別為fname1.in，fname2.in，fname3.

4、in。打開記事本，寫入以下兩行直線間的內容，另存為*.bat文件，例如：GPR.bat。其中括號里的內容為本條語句的說明。-path C:;%path% (把軟件所在路徑包含進系統搜索路徑)cd/d D:GPR (設置專用文件夾路徑為活動路徑)gprmax2d D:GPRfname1.in （按順序進行批量正演模擬）gprmax2d D:GPRfname2.ingprmax2d D:GPRfname3.in-完成后保存并關閉，然后雙擊運行*.bat文件即可進行批量正演模擬四、建立輸入文件輸入文件是純文本文件，可以用任何文字處理程序編輯，而且內容全部為英文。打開記事本，寫入輸入文件內容，另存為*

5、.in即可。輸入文件必需包括模擬模型的所有必要信息，每一個有效命令行都必須以符號（#）開始，否則該命令行無效。所有命令的一般語法為：#命令名: 參數1 參數2 參數3 為了方便介紹GprMax2D的命令及參數，做了以下約定：·f 代表實數，如1.5·i 代表整數，如15·c 代表單個字母，如y·str 代表字條串，如air·file 代表文件名，如test.in·所有的長度單位均為m·所有的時間單位均為s·所有的頻率單位均為Hz·模型空間的坐標原點（0,0）在左下角輸入文件中除了#analysis:和#e

6、nd_analysis:命令外，先后順序一般是沒有要求的。GprMax2D 2.0版一共有32個命令：-#title: str輸出文件的標題，包含在輸出文件中。#domain: f1 f2指定模型空間大小，單位是米。·f1 x方向的大小·f2 y方向的大小#dx_dy: f1 f2指定網格步長，單位是米。·f1 x方向的步長，即x·f2 y方向的步長，即y模型劃分的網格數目是由#domain:和#dx_dy:兩個命令共同決定的。網格劃分越小模型越精確，但模擬的計算量也越在，需要權衡確定。#time_step_stability_factor: f1時間步

7、長的穩定系數，取值范圍（0，1，實際使用時間為f1×t。#time_window: f1 or i1時窗大小，單位是秒。·f1 為時窗，如f1=16e-9，時窗為16ns·i1 為時間步長數，時窗為i1×t。#messages: c1是否開啟屏幕信息。y(es)開啟n(o)關閉，默認n。#number_of_media: i1介質數目，i1應大于10。系統分配了10種介質（自由空間即空氣free_space，金屬pec和8種用戶自定義介質）的內存空間，當模型中使用介質數量超過10種時需要使用該命令。#nips_number: i1指定分配的空間。存儲重要

8、信息模型的數組大小可以在內部計算，但通常比實際需要的大，可以使用較小的數字節省空間。如果分配過小，系統會發出錯誤并提醒分配更多空間#media_file: file1要包含的介質文件的文件名（如果有必要還包括文件路徑）。介質文件格式為：# f1 f2 f3 f4 f5 f6 str1參數的具體含義參考#medium:命令。#geometry_file: file1存儲模型文件（*.geo）。#medium: f1 f2 f3 f4 f5 f6 str1在輸入文件中定義介質屬性。·f1 靜態時介質的相對介電常數·f2 在理論上的無限頻率的相對介電常數·f3 介質的弛

9、豫時間（秒）·f4 靜態時介質的電導率（西門子/米）·f5 相對磁導率的介質·f6 介質的磁導率·str1 介質的名稱如果不想使用debye介質，可設置f3為0.0，此時系統將只使用f1和f4來描述介質介電性能。如果介質為非磁性，則設置f5為1.0，f6為0.0。如果使用f3的值，該值應始終高于模型中使用的時間步長t。-該組命令不建議使用，除非是熟悉ABC的人用來定制和優化吸收邊界條件。#abc_type: pml恢復默認的PML邊界。#abc_order: i1ABC次序。取值范圍1，2，3，默認為3。次序越低，性能越差。#abc_stability_

10、factors: f1 f2 f3指定ABC的穩定系數。#abc_optimization_angles: f1 f2 f3優化ABC為特定的入射角。#abc_mixing_parameters: f1 f2 f3 f4 f5 f6改變使用ABC的時間和偏導數的離散化。默認使用box。#pml_layers: i1指定PML使用的網格數量。默認為8。數量越多性能越佳。-#box: f1 f2 f3 f4 str1使用定義好的介質定義一個矩形模型。·f1 矩形左下角x坐標，單位米·f2 矩形左下角y坐標，單位米·f3 矩形右上角x坐標，單位米·f4 矩形右

11、上角y坐標，單位米·str1 使用的介質的名稱#cylinder: f1 f2 f3 str1使用定義好的介質定義一個圓形模型。·f1 圓心x坐標，單位米·f2 圓心y坐標，單位米·f3 圓半徑R，單位米·str1 使用的介質的名稱#x_segment: f1 f2 f3 f4 f5 str1使用定義好的介質定義一個圓沿x軸裁剪的模型。·f1 圓心x坐標，單位米·f2 圓心y坐標，單位米·f3 沿x軸開始的位置·f4 沿x軸結束的位置·f5 圓半徑R，單位米·str1 使用的介質的名稱

12、#y_segment: f1 f2 f3 f4 f5 str1使用定義好的介質定義一個圓沿y軸裁剪的模型。·f1 圓心x坐標，單位米·f2 圓心y坐標，單位米·f3 沿y軸開始的位置·f4 沿y軸結束的位置·f5 圓半徑R，單位米·str1 使用的介質的名稱#triangle: f1 f2 f3 f4 f5 f6 str1·f1 第一頂點x坐標，單位米·f2 第一頂點y坐標，單位米·f3 第二頂點x坐標，單位米·f4 第二頂點y坐標，單位米·f5 第三頂點x坐標，單位米·f6

13、 第三頂點y坐標，單位米·str1 使用的介質的名稱模型的建立是有先后順序的，后建立的模型會把先建立的模型覆蓋住。-該組命令是數值模擬的核心。一個輸入文件中可以出現多組數值模擬命令，但是數值模擬的控制命令必須放在#analysis:和#end_analysis:之間。#analysis: i1 file1 c1數值模擬開始的命令。·i1 總步數·file1 模擬結果輸出文件名（*.out）·c1 類型，取值范圍a（文本文件）或b（二進制文件）#end_analysis:數值模擬結束的命令，無參數。#tx: f1 f2 str1 f3 f4發射天線的屬性。

14、·f1 發射天線位置的x坐標，單位米·f2 發射天線位置的y坐標，單位米·str1 定義過的發射天線的名稱·f3 發射天線的延遲時間·f4 發射天線的清除時間#rx: f1 f2接收天線的屬性。·f1 接收天線位置的x坐標，單位米·f2 接收天線位置的y坐標，單位米#rx_box: f1 f2 f3 f4 f5 f6接收區域的屬性。·f1 接收區域左下角的x坐標，單位米·f2 接收區域左下角的y坐標，單位米·f3 接收區域右上角的x坐標，單位米·f4 接收區域右上角的y坐標，單位米&#

15、183;f5 接收區域在x軸的步長，即x，單位米·f6 接收區域在y軸的步長，即y，單位米#snapshot: i1 f1 f2 f3 f4 f5 f6 f7 file1 c1或#snapshot: i1 f1 f2 f3 f4 f5 f6 i2 file1 c1獲取模型一個區域內，某一時刻的電磁場的快照。·i1 第幾道的掃描，值介于1和#analysis:命令中的總步數。·f1 區域左下角的x坐標，單位米·f2 區域左下角的y坐標，單位米·f3 區域右上角的x坐標，單位米·f4 區域右上角的y坐標，單位米·f5 區域在x軸

16、方向上的采樣間隔，單位米·f6 區域在y軸方向上的采樣間隔，單位米·f7 快照的時刻，單位秒·i2 快照的迭代次數·file1 快照存儲的文件名（*.out）·c1 類型，取值范圍a（文本文件）或b（二進制文件）#tx_steps: f1 f2發射天線的移動步長。·f1 發射天線沿x軸的移動步長，單位米·f2 發射天線沿y軸的移動步長，單位米#rx_steps: f1 f2接收天線的移動步長。·f1 接收天線沿x軸的移動步長，單位米·f2 接收天線沿y軸的移動步長，單位米-#line_source: f1

17、 f2 str1 str2·f1 發射天線的振幅·f2 發射天線的頻率·str1 信號波形，有五種： ·cont_sine 連續正弦波 ·sine 脈沖正弦波 ·gaussian 高斯波 ·ricker 雷克子波 ·user 用戶自定義的激勵函數·str2 發射天線的名稱#excitation_file: str1選擇用戶自定義的激勵文件。自定義激勵文件是文本文件。-五、MATLAB繪圖筆者使用MATLAB7.0，以該版本為例進行說明，其他版本可參考使用。為方便說明，筆者假設讀者已對MATLAB具有最基本

18、的了解和使用能力，同時主要以命令來完成操作。讀者若對MATLAB一無所知，請另行學習。1.設定當前目錄和搜索路徑。-cd D:GPR %設定當前目錄 addpath D:GPR %設定本次搜索路徑-2.模型文件繪圖-geo = 'fname1.geo'mesh,header,media = gprmax2g(geo);M,N = size(mesh); %獲取模型大小imagesc(1:N)*header.dx,(1:M)*header.dy,mesh) %繪圖xlabel('Distance/m'); %設置x軸標簽ylabel('Depth/m

19、9;); %設置y軸標簽-3.輸出文件繪圖-out = 'fname1.out'Header,Fields = gprmax(out);NN = 1:Header.NSteps; %移動次數Position = Header.dx*Header.tx+(NN-1)*Header.dx*Header.TxStepX); %天線每次所在位置Data(:,:) = Fields.ez(:,1,:); %轉換數組格式imagesc(Position,Fields.t*1e9,Data) %繪圖colorbar; %顯示色彩條xlabel('Distance/m'); %

20、設置x軸標簽ylabel('Time/ns'); %設置y軸標簽-繪出的圖像請及時保存。例子：#medium: 6.0 0.0 0.0 0.01 1.0 0.0 concrete#medium: 20.0 0.0 0.0 0.1 1.0 0.0 wet_sand-#domain: 2.5 0.5#dx_dy: 0.0025 0.0025#time_window: 12e-9-#box: 0.0 0.0 2.5 0.45 wet_sand- Get a cylinder of free space and then put a slab of concrete to cut it in half -#cylinder: 1.05 0.3 0.2 free_space #box: 0.0 0.3 2.5 0.45 concrete#cylinder: 0.25 0.375 0.0125 pec#cylinder: 0.45 0.375 0.0125 pec#cylinder: 0.65