李展輝，曹學峰，朱自強

(1.自然資源部 航空地球物理與遙感地質(zhì)重點實驗室，北京 100083；2.中國自然資源航空物探遙感中心，北京 100083)

1 引 言

探地雷達(Ground Penetrating Radar, GPR)是一種在工程勘察、地質(zhì)調(diào)查、路基檢測等多領(lǐng)域有廣泛應用的淺層電磁探測方法[1-6]。探地雷達最主要的數(shù)據(jù)成像方法是偏移成像技術(shù)。隨著科技的不斷發(fā)展以及任務需求的不斷精細化，偏移成像已經(jīng)從最初的相移偏移成像、克希霍夫偏移成像、波動方程偏移成像，發(fā)展到了現(xiàn)今的三維逆時偏移成像，甚至是三維全波形反演。不論是逆時偏移成像還是全波形反演，對三維正演的需求都非常大[7-10]。此外，探地雷達的正演模擬亦能直觀反應典型目標體的回波響應特征，對于目標體的形態(tài)、位置等基本屬性分析具有重要意義[11]。為滿足繁重的三維正演計算任務，最常用的方法是采用時域有限差分(Finite Difference Time Domain, FDTD)結(jié)合MPI、OpenMP等加速技術(shù)，充分利用現(xiàn)今計算機集群的并行計算能力結(jié)合FDTD良好的并行性質(zhì)，開展快速的三維正演計算。尤其是基于圖形處理器(Graphics Processing Unit, GPU)的CUDA加速的算法已經(jīng)被廣泛采用[12-15]。其主要原因是在同等價格、功耗條件下，GPU的浮點計算性能已經(jīng)遠超過CPU的浮點計算性能。以消費級GPU和CPU對比為例，NVIDIA公司的RTX 2080Ti顯卡單精度計算能力約為13.5 TFOLPS，而當前同價位單精度計算能力最強的CPU，AMD Threadripper 3960X，其單精度計算能力僅為3.1 TFLOPS左右，不足RTX 2080Ti的1/4。此外，當今世界上許多頂級超級計算機，比如美國的Summit和Sierra，浮點計算力也主要由GPU提供。不過，GPU雖然能提供的算力遠高于CPU，但是其可用顯存是極為有限的，RTX 2080Ti的顯存僅為12 GB，而Threadripper 3960X最高可支持到256 GB內(nèi)存。因此在消費級的計算機上，一些較大的問題比如1 000×1 000×1 000網(wǎng)格維度的FDTD正演問題，使用單精度僅存儲電磁場的6個分量就需要24 GB的內(nèi)存容量，在消費級GPU上沒有運算的可能，只能利用CPU進行計算。

除了FDTD之外，在探地雷達正演中比較常用的還有時域有限元法(Finite Element Time Domain, FETD)[16, 17]、時域偽譜法(Pseudo-Spectral Time Domain, PSTD)[18-21]等方法。其中PSTD使用快速傅里葉變換(Fast Fourier Transform, FFT)開展電磁場空間一階導數(shù)的計算，具有遠高于FDTD的空間導數(shù)計算精度[18]，曾在探地雷達正演計算方面有著廣泛的應用。但是隨著并行計算要求的提高，PSTD因其在分布式計算機集群上并行性能較差等弱點，逐漸被FDTD等并行性能好的算法替代。不過在單個CPU計算節(jié)點或者單個GPU上，PSTD依然可以利用多個FFT計算的獨立性，開展多線程并行運算。由于PSTD的空間導數(shù)高精度特性，其空間網(wǎng)格采樣間隔可以接近奈奎斯特采樣定理，即波長的1/2，而FDTD一般要求空間采樣間隔不大于波長的1/10[22]，即使4階或更高階的FDTD方法，其采樣間隔一般不大于波長的1/4。那么使用PSTD進行正演計算，1 000×1 000×1 000網(wǎng)格維度問題理論上至少可以縮減為500×500×500網(wǎng)格維度問題，存儲6個電磁分量所需的內(nèi)存縮減到原來的1/8，僅為3 GB，完全可以在消費級的GPU上開展計算。而且由于網(wǎng)格數(shù)量的大幅減少，以及網(wǎng)格間距的增大，計算所需的時間也能夠大幅縮短。

基于以上問題與分析，本文在探地雷達三維正演問題中引入基于GPU加速的三維PSTD正演方法，相比FDTD，在保持相同精度的同時，能大幅減小正演對顯存的需求，充分利用GPU的高性能浮點計算能力，提高探地雷達三維正演效率。

2 原理與方法

探地雷達電磁波的激發(fā)與傳播由Maxwell方程組控制。在給定合理的初始條件后，Maxwell方程組中的散度方程組自動滿足。因此，在正演計算的過程中，一般只考慮兩個旋度方程組：

其中，E為電場，H為磁場，ε為介電常數(shù)，σ為電導率，μ為磁導率，J為電流激勵源。介電常數(shù)又通常用相對介電常數(shù)表述ε≡εrε0，其中ε0為真空中介電常數(shù)，εr為相對介電常數(shù)。探地雷達正演一般采用吸收邊界條件，本文采用非分裂的時域卷積完全匹配層(Perfectly Matched Layer, PML)吸收邊界[23]，相關(guān)理論與實施方法已經(jīng)非常成熟，本文不再贅述。電流激勵源為高斯函數(shù)一階導數(shù)類:

J(t)=2πfct×e-(πfct)2

(3)

其中,fc為激勵源的中心頻率。

模型的網(wǎng)格采用最常用的Yee交錯網(wǎng)格。在交錯網(wǎng)格中，PSTD相對于FDTD的主要區(qū)別在于PSTD使用FFT求取交錯網(wǎng)格點上的空間一階導數(shù)。一般形式的FFT求導可表述如下：

設(shè)有一場f(x)，其FFT變換和反變換分別為

式(4)、式(5)中，i為虛數(shù)單位，k為波數(shù)。那么f(x)在x+Δx處的導數(shù)可表示為

(6)

即可以通過求ikeikΔxF(k)的FFT反變換來獲取f(x)在x+Δx處的導數(shù)。

PSTD在GPU中的加速通過NVIDIA公司的CUDA Toolkit來實現(xiàn)。CUDA Toolkit中的CUFFT庫函數(shù)在GPU上基于著名的FFT庫函數(shù)FFTW實現(xiàn)了FFT的并行計算，并進行了大量的優(yōu)化。在PSTD中，主要用到其中的一維FFT正變換函數(shù)cufftExecR2C、反變換函數(shù)cufftExecC2R，以及FFT計算空間開辟函數(shù)cufftPlan1D，操作流程如圖1所示。

圖1 使用CUDA Toolkit中的FFT庫函數(shù)求導流程圖

CUDA Toolkit中CUFFT庫函數(shù)有一個重要的特征，就是可以直接輸入三維數(shù)組開展一維FFT變換，但是要求開展FFT變換的維度處于第1個維度。比如有一三維電場Ex(i,j,k)，其中i代表x方向序號，j代表y方向序號，k代表z方向序號，三個方向維度大小分別為Nx、Ny、Nz。要計算其在y方向上的空間一階導數(shù)，那么需要將y方向的維度轉(zhuǎn)置到第1個維度，然后整體輸入到圖1所示的流程中，同時將圖1中的Howmany設(shè)置為Nx×Nz，表示開展Nx×Nz次沿y方向的空間一階導數(shù)計算。三維數(shù)據(jù)的直接輸入能夠大幅減少函數(shù)的調(diào)用次數(shù)，進而顯著提高計算效率。出于此目的，為了盡量減少計算過程中的電磁場轉(zhuǎn)置計算，將Ex、Ey、Ez的維度次序分別設(shè)置為Ex(j,k,i)，Ey(k,i,j),Ez(i,j,k)，磁場采用與電場同樣的維度次序。在該維度次序下，轉(zhuǎn)置計算的復雜度能夠顯著減少，比如Ex僅需在求對z方向空間導數(shù)時，開展第1、2維度的轉(zhuǎn)置計算，不涉及第3維度的轉(zhuǎn)置，其他分量類似。在電磁場時間迭代計算過程中，此維度次序設(shè)計亦能顯著降低在第三個維度上的頻繁尋址，進而提高正演效率。

在求出空間導數(shù)之后，時間步進迭代采用二階中心差分格式。二階中心差分格式已經(jīng)非常成熟，以Ex為例，第n+1步的迭代方式可表述為

(7)

3 數(shù)值對比與算例

3.1 基于GPU與CPU的FFT并行效率對比

為對比GPU并行相比CPU并行的FFT計算效率，本文以常規(guī)消費級筆記本為基礎(chǔ)，對比計算效率。計算機平臺為微星GT72，其中CPU模塊為Intel I7 5770HQ，全核心單精度計算能力為0.45 TFLOPS，GPU模塊為NVIDIA GeForce 970 M，單精度計算能力為2.7 TFLOPS。測試內(nèi)容為對一三維數(shù)組沿其中第1個維度開展1 000次FFT正變換與反變換，三維數(shù)組的大小如表1所示。CPU并行計算采用當前最高效的Intel Fortran編譯器以及FFTW庫函數(shù)，GPU并行采用PGI編譯器中CUDA Fortran編譯模塊，并調(diào)用CUDA Toolkit 10.1中的CUFFT庫函數(shù)進行編譯與運行。運行結(jié)果如表1所示。與GPU和CPU的單精度計算能力差距相似，不管哪種情形下，基于GPU的FFT運算速度都顯著高于基于CPU的，雖然未達到單精度計算能力的倍率，但是基本上達到了4倍以上，具有非常高的加速效果。

表1 分別基于GPU與CPU并行的三維數(shù)組的一維FFT計算效率對比.

3.2 PSTD與4階FDTD同網(wǎng)格精度對比

盡管PSTD相對FDTD的高精度特性已有大量的對比工作，本文從存儲節(jié)約方面再次說明PSTD相對FDTD的重要意義。由于二階FDTD精度過低，因此本文與常用的空間四階FDTD進行對比。對比所用模型如圖2所示。該模型實際上是一個二維模型，電性參數(shù)沿z方向不變，但是采用三維網(wǎng)格進行離散，三維源激發(fā)。背景介質(zhì)為空氣，即相對介電常數(shù)為1、電導率為0，左上角為一相對介電常數(shù)為9、電導率為0的低速異常體，內(nèi)部電磁波速為背景介質(zhì)中的1/3，右下角為一相對介電常數(shù)為4、電導率為0的低速異常體，內(nèi)部電磁波速為背景介質(zhì)中的1/2。激勵源位于模型的正中間，中心頻率fc為100 MHz，x、y、z三方向的網(wǎng)格采樣間距均為0.2 m，網(wǎng)格大小為400×200×200，時間步長設(shè)定為0.24 ns，同時滿足PSTD和四階FDTD的穩(wěn)定性條件。

圖2 PSTD與FDTD正演效果對比所用異常體模型

不同時間點PSTD與四階FDTD的正演結(jié)果在z軸中間的x-y切面上的不同時刻波場快照如圖3所示。在43 ns時，電磁波尚未傳播到異常體位置，二者的結(jié)果是一致的。在91 ns時，電磁波已經(jīng)初步進入兩個異常體內(nèi)部，可以看出，在兩個異常體內(nèi)部PSTD依然體現(xiàn)出很高的精度，沒有出現(xiàn)明顯的空間頻散現(xiàn)象，而對于四階FDTD，左上角低速異常體中已經(jīng)出現(xiàn)了較大的空間頻散現(xiàn)象，右下角的異常體內(nèi)部尚未出現(xiàn)頻散問題。隨著時間的進一步推進，在139 ns時，PSTD依然在所有區(qū)域完全正常，而FDTD在左上角的低速區(qū)已經(jīng)有非常嚴重的頻散問題，在右下角的異常體中也隨著電磁波的推進出現(xiàn)了頻散問題。圖4(a)進一步對比了在t= 43 ns時y=0 m剖面線上Ez剖面，PSTD與FDTD基本一致，也說明在網(wǎng)格間距相對波長較小的時候，二者都可以獲取到好的正演結(jié)果；圖4(b)和圖4(c)分別對比了在t= 139 ns時，y= 5 m(穿過右下角異常體)和y= 35 m(穿過左上角異常體)剖面線上Ez剖面，由于在異常體內(nèi)部波長分別壓縮為空氣中的1/2(左下角異常體)和1/3(右上角異常體)，F(xiàn)DTD由于采樣不足而出現(xiàn)了嚴重的空間頻散問題，PSTD則始終維持著非常好的波形形態(tài)?？梢?，PSTD相比4階FDTD表現(xiàn)出了其超高的空間精度特征，4階FDTD至少需要加密一倍以上的網(wǎng)格才能接近PSTD的空間精度，但是如此一來，存儲占用要增加為原來的8倍，而且會超出消費級GPU顯存的容量，無法再利用普通GPU加速，同時時間步長也要減半，計算效率大幅降低?？梢?，PSTD在節(jié)約存儲方面有非常大的優(yōu)勢，非常適合于GPU顯存較小的特性。

圖3 PSTD(a、c、e)與四階FDTD(b、d、f)正演結(jié)果在z軸中間x-y切面上不同時刻的波場快照對比

圖4 PSTD與四階FDTD正演結(jié)果在y=0 m, t=43 ns(a);y=5 m t=139 ns(b);y=35 m, t=139 ns(c)剖面線上的Ez剖面對比

4 地下三維空洞模型GPR探測的PSTD模擬

為進一步體現(xiàn)基于GPU加速的PSTD正演方法的實用性，設(shè)計了如圖5(a)所示的地下空洞探測模型。模型左側(cè)和覆蓋層是一個相對介電常數(shù)為4、電導率為0.001 S/m的基巖體，右下側(cè)為一個相對介電常數(shù)為7、電導率為0.002 S/m的弱破碎帶，左右兩側(cè)各有一個埋深為10 m、直徑為4 m的球形空洞。在空洞上方的地表，從左至右以0.8 m為間隔，布置有85個自激自收GPR探測點，如圖5(a)中地表的黑色小點所示。網(wǎng)格剖分在x、y、z方向的數(shù)量分別為400、150、150，網(wǎng)格間距為0.2 m，天線中心頻率為100 MHz，時間步長為0.24 ns，總時間步數(shù)為1 600，數(shù)據(jù)采集總時間為384 ns。采用NVIDIA GeForce 970M進行GPU加速運算，經(jīng)統(tǒng)計，平均每個測點的正演時間為6.5 min，總運行時間為9.2 h。如采用RTX 2080Ti等更高級的GPU加速卡，預計還能獲取到4倍以上的加速性能，整個模擬過程在2 h左右便能完成，能為GPR數(shù)據(jù)的三維逆時偏移成像與全波形反演的快速開展提供新的并行加速方案。

GPR探測模擬結(jié)果的測點-時間數(shù)據(jù)成像如圖5(b)所示。為清楚顯示地下空洞的回波，每一道數(shù)據(jù)都進行了一定程度的道內(nèi)均衡增益，這也是GPR在數(shù)據(jù)展示過程中常用的手段?？梢钥闯?，即使采用了0.2 m的大尺寸網(wǎng)格，模擬結(jié)果也沒有頻散現(xiàn)象的出現(xiàn)。以激勵源中心頻率的2倍作為最高頻率，即200 MHz，那么在右側(cè)低速區(qū)內(nèi)的最小波長為0.57 m，0.2 m的網(wǎng)格尺寸意味著每個波長采樣不足3個。PSTD在該條件下依然能夠維持長時間內(nèi)無明顯空間頻散現(xiàn)象，再次證明了其在精度方面的巨大優(yōu)勢。兩個地下空洞的反射信號組成了典型的雙曲線型回波，不過由于所處背景介質(zhì)的不同，信號強度、回波時間均有不同，而且摻雜有接觸角點引起的回波和接觸界面引起的多次波，更準確的目標識別需要通過信號處理、偏移成像等手段來實現(xiàn)。

圖5 地下空洞模型與GPR探測模擬結(jié)果

5 結(jié) 論

本文將基于GPU加速技術(shù)的PSTD正演方法應用到探地雷達三維模擬中，充分利用PSTD高空間精度的優(yōu)勢壓縮顯存需求，并開發(fā)非常規(guī)的三維數(shù)組維度設(shè)置方案減少轉(zhuǎn)置計算復雜度與第三維度尋址次數(shù)，高效利用了GPU的超高浮點計算性能，實現(xiàn)了探地雷達三維正演的GPU大幅度加速。本文研究表明，在消費級的計算機上，GPU加速相比CPU并行計算能夠獲取接近于單精度浮點計算倍率的加速效率；PSTD相比FDTD能夠更充分地綜合利用GPU高浮點計算能力和有限顯存的特征，采用較粗的網(wǎng)格、較少的內(nèi)存高效獲取高精度的正演結(jié)果，對于開展三維逆時偏移成像和全波形反演具有重要意義。