基于GPU的相位梯度變化質(zhì)量圖并行計算*

李 偉 劉殿剛 鐘何平

(1.91458部隊 三亞 572021)(2.海軍工程大學(xué)電子工程學(xué)院 武漢 430033)

?

基于GPU的相位梯度變化質(zhì)量圖并行計算*

李 偉1劉殿剛1鐘何平2

(1.91458部隊 三亞 572021)(2.海軍工程大學(xué)電子工程學(xué)院 武漢 430033)

相位梯度變化質(zhì)量圖是一種可靠的相位質(zhì)量圖,但當(dāng)計算窗口較大時,計算效率非常低下,嚴(yán)重影響信號處理的速度。論文提出了一種快速的相位梯度變化質(zhì)量圖的計算方法,采用GPU的多個計算核心同時計算分塊纏繞相位的相位質(zhì)量圖,然后將結(jié)果下載到主機(jī)內(nèi)存。最后通過試驗驗證了所提方法的高效性。

質(zhì)量圖; 相位梯度; 并行計算; GPU

Class Number P237

1 引言

干涉合成孔徑雷達(dá)(InSAR)和干涉合成孔徑聲納(InSAS)都是通過兩個信號接收器接收來自同一區(qū)域的回波,得到單視復(fù)數(shù)影像數(shù)據(jù),后經(jīng)成像處理和干涉運(yùn)算得到干涉圖像,進(jìn)一步反演出目標(biāo)區(qū)域高程信息[1～2]。兩幅精確配準(zhǔn)后的復(fù)數(shù)影像經(jīng)共軛相乘后形成干涉相位圖,干涉相位圖是一個二維矩陣,其中每一點的數(shù)值表示的是空間上一點與兩個傳感器之間因距離差引起的纏繞相位差。理想上,所有相鄰的干涉相位點之間滿足Nyquist采樣定理,但是實際中由于透視壓縮、陰影和層疊,從而引起兩幅圖像相干性下降,以及由其他原因引起的去相關(guān)現(xiàn)象和在信號處理中引入的相干噪聲等,使得部分區(qū)域不滿足采樣定理,干涉信號處理過程變得復(fù)雜化[3～4]。質(zhì)量圖[5]就是表示每一點干涉相位質(zhì)量好壞的數(shù)組,可以用于指導(dǎo)枝切法中枝切線的設(shè)置,質(zhì)量引導(dǎo)相位解纏算法[6～8]中路徑的選擇,全局優(yōu)化解纏算法[9～11]中權(quán)重的設(shè)置,在相位解纏過程中扮演著重要的作用。但是隨著計算窗口的增加,計算效率變得越來越低嚴(yán)重影響著相位解纏的速度。GPU的出現(xiàn)為高速進(jìn)行相位梯度質(zhì)量圖的計算提供了新的手段,它的優(yōu)點恰好是在處理高密度的數(shù)值運(yùn)算,并且已在天文計算、流體力學(xué)模擬、圖像處理、音視頻編碼解碼領(lǐng)域廣泛應(yīng)用,取得了顯著的效率提升。

本文借助GPU的強(qiáng)大計算能力,提出了一種GPU環(huán)境下的快速相位梯度變化質(zhì)量圖的計算方法。首先將干涉圖整塊上傳到GPU的顯存空間,然后將整個干涉圖分為小塊,每個線程塊分別計算對應(yīng)的小塊干涉相位的相位梯度變化質(zhì)量圖,每個線程塊完成計算后,將計算結(jié)果進(jìn)行拼接,然后下載到主機(jī)端的內(nèi)存。最后通過計算不同局部窗口大小下的相位梯度變化質(zhì)量圖,驗證了所提方法的高效性。

2 相位梯度變化質(zhì)量圖

InSAR和InSAS干涉測量系統(tǒng)是以信號的相干性為基礎(chǔ)的,相干性可以用相干系數(shù)來衡量。得到干涉圖后,需要對相位數(shù)據(jù)的一致性進(jìn)行分析,以確定哪些區(qū)域相位是一致的,哪些區(qū)域相位是不一致的。最常用的質(zhì)量圖就是相干圖,相干圖中相干值的高低表明了干涉圖不同區(qū)域的相干性大小,是最直觀的干涉質(zhì)量評價圖。相干系數(shù)定義如下:

(1)

式中u1和u2分別代表兩幅單視復(fù)圖像,*表示取共軛,E[·]表示期望值,γ的絕對值范圍是[0,1]。γ=0表示兩幅影像完全不相關(guān),γ=1表示兩幅影像完全一樣。

實際計算時,假設(shè)在k×k窗口內(nèi)的估值區(qū)域內(nèi)隨機(jī)過程是平穩(wěn)和各態(tài)歷經(jīng)的,即相干系數(shù)可以用樣本的空間平均值按照式(2)計算。

(2)

逐點計算以(m,n)為中心的k×k窗口內(nèi)的γ值,其結(jié)果是相關(guān)系數(shù)文件,每一點相關(guān)系數(shù)的絕對值表示該點配準(zhǔn)結(jié)果的好壞或干涉相位可信度的大小。

在無法獲得相關(guān)圖時,常采用相位梯度變化質(zhì)量圖作為相位質(zhì)量的評價指標(biāo)。相位梯度變化質(zhì)量圖定義如下:

(3)

3 質(zhì)量圖的并行計算

將干涉相位數(shù)組表示為φ(m,n),1≤m≤M,1≤n≤N,其中M、N分別表示干涉相位數(shù)組的行數(shù)和列數(shù)。對于數(shù)組中的任一點(i,j),其相位質(zhì)量值q(i,j)是根據(jù)以點(i,j)為中心,2Wx+1和2Wy+1為窗寬的局部干涉相位計算所得。圖1為干涉相位質(zhì)量圖計算示意圖,待計算的干涉相位點坐標(biāo)為(6,6),局部窗口大小設(shè)置為5×5,取得當(dāng)前局部窗口內(nèi)的相位質(zhì)量值后,立刻調(diào)用相位梯度變化質(zhì)量函數(shù)計算對應(yīng)的相位質(zhì)量值。完成質(zhì)量值得計算后,將其寫入到對應(yīng)的結(jié)果數(shù)組中。相位質(zhì)量值的計算過程非常簡單,但是其計算密度非常高,特別是局部窗口加大時,效率顯著降低,嚴(yán)重影響著干涉信號的后續(xù)處理。

圖1 質(zhì)量值計算示意圖

圖2 計算流程

GPU為實現(xiàn)快速計算相位梯度變化質(zhì)量圖提供了可能,采用GPU計算相位梯度變化質(zhì)量圖的流程如圖2所示。具體流程如下:

1) 從原始數(shù)據(jù)中讀入干涉相位圖;

2) 根據(jù)干涉相位圖的行、列數(shù),調(diào)用cudaMalloc函數(shù),在GPU上分配顯存;

3) 在主機(jī)端調(diào)用cudaMemcpy函數(shù),將主機(jī)內(nèi)存中存放的干涉相位圖上傳至新開普的顯存空間;

4) 通過在主機(jī)端調(diào)用核函數(shù),啟動相位梯度變化質(zhì)量圖的并行計算函數(shù),計算后的結(jié)果存儲在臨時顯存中;

5) 調(diào)用cudaMemcpy函數(shù),將保存在顯存中的計算結(jié)果拷貝到主機(jī)內(nèi)存;

6) 釋放GPU上分配的顯存空間。

4 實驗結(jié)果分析

為了驗證基于GPU的干涉相位質(zhì)量圖計算效率,在如下環(huán)境進(jìn)行了相位梯度變化質(zhì)量圖的計算試驗:處理器Intel(R) Xeon(R) CPU X5650;顯卡Tesla C2050;操作系統(tǒng)Windows 7專業(yè)版;軟件環(huán)境VS2008+CUDA5.0。

試驗數(shù)據(jù)選用了一幅TerraSAR干涉數(shù)據(jù),如圖3(a)所示,其大小為1024×1024像素。從干涉圖中可以看出成像區(qū)域存在大量的噪聲區(qū)域,嚴(yán)重影響著后續(xù)的相位解纏。為了比較不同大小干涉圖質(zhì)量圖的計算效率,將圖3(a)從左上角開始,依次截取大小為512×512、256×256的兩幅子圖,分別如圖3(b)和3(c)所示。

表1給出了CPU和GPU計算方式下相位梯度變化質(zhì)量圖的計算時間,其中質(zhì)量圖的計算時間取10次計算時間的平均值,數(shù)據(jù)上傳和數(shù)據(jù)下載時間取的是不同計算窗口大小下的三次數(shù)據(jù)傳輸時間平均值。從表1中可以看出采用CPU計算不需要進(jìn)行數(shù)據(jù)傳輸,采用GPU計算需要進(jìn)行數(shù)據(jù)的上傳和下載操作,但由于傳輸帶寬大,可以在非常短的時間完成。在CPU計算環(huán)境下,干涉圖大小為1024×1024,計算窗口為7×7時,最大相位梯度質(zhì)量圖的計算時間為1291ms。但是在GPU計算環(huán)境下,相位梯度變化質(zhì)量圖的計算時間與CPU中計算時間相比顯著降低。維數(shù)為1024×1024的干涉圖,當(dāng)計算窗口為7×7時,計算時間為20.8ms。

圖4 加速比

相位梯度變化質(zhì)量圖的計算加速比如圖4所示,在計算加速比時,GPU計算時間中加入了數(shù)據(jù)上傳和下載時間。對于256×256數(shù)據(jù)量小的干涉圖,在計算窗口為3×3時,塊之間的重疊比例較高,此外由于數(shù)據(jù)傳輸時間的影響,加速比不是非常理想。相位梯度變化質(zhì)量圖的加速比隨著計算窗口和干涉圖維數(shù)的增加而增加,當(dāng)干涉圖維數(shù)為1024×1024,計算窗口為7×7,其加速比高達(dá)50。

表1 相位梯度變化質(zhì)量圖計算效率比較

5 結(jié)語

本文提出了一種基于GPU的快速相位梯度變化質(zhì)量圖計算方法。首先將干涉圖上傳至顯存,然后采用多個線程塊同時計算每個分塊干涉相位的質(zhì)量值,并將每個線程塊計算的質(zhì)量值進(jìn)行合并。最后將計算的相位質(zhì)量圖下載到主機(jī)內(nèi)存。試驗結(jié)果表明,對于大小為1024×1024的干涉圖,采用GPU計算相位梯度變化質(zhì)量圖,可以獲得高達(dá)50倍的加速比。

[1] 李平湘,楊杰.雷達(dá)干涉測量原理與應(yīng)用 [M].北京:測繪出版社,2006.

[2] Ghiglia D C, Pritt M D. Two-dimensional phase unwrapping: theory, algorithm, and software[M]. New York: John Wiley & Sons. Inc,1998.

[3] 彭石寶,袁俊泉,向家彬.一種基于加權(quán)迭代貪婪算法的InSAR相位解纏的新方法[J].電子與信息學(xué)報,2008,30(6):1326-1330.

[4] Costantini M. A novel phase unwrapping method based on network programming[J]. IEEE Transactions on Geoscience and Remote Sensing,1998,36(3):813-821.

[5] 楊夏,于起峰,伏思華.由新質(zhì)量圖引導(dǎo)的InSAR快速解纏方法[J].電子與信息學(xué)報,2007,29(10):2367-2370.

[6] 鐘何平,唐勁松,陳鳴.一種新的InSAR相位質(zhì)量圖生成方法[J].遙感學(xué)報,2011,15(4):766-777.

[7] 陳家鳳,李春芳.基于自適應(yīng)加權(quán)窗口的相位質(zhì)量圖提取[J].光學(xué)技術(shù),2012,38(6):734-739.

[8] Cho B. L., Kong Y. K., Kim Y. S. Quality map extraction for radar interferometry using weighted window[J]. Electronics Letters,2004,40(8):472-473.

[9] Costantini M. A novel phase unwrapping method based on network programming[J]. IEEE Transactions on Geoscience and Remote Sensing,1998,36(3):813-821.

[10] Yang Fengtao, Lu Xiaoxu, Wang Dianyuan, et al. Weighted minimum cost flow phase unwrapping algorithm based on second difference[J]. Laser Technology,2006,30(6):667-672.

[11] Ghiglia, D. C., Romero, L. A. Robust two-dimensional weighted and unweighted phase unwrapping that uses fast transforms and iterative methods[J]. Journal of the Optical Society of America A: Optics and Image Science, and Vision,1994,11(1):107-117.

Parallel Computing of Phase Derivative Variance Quality Map Based on GPU

LI Wei1LIU Diangang1ZHONG Heping2

(1. No. 91458 Troops of PLA, Sanya 572021) (2. Electronics Engineering College, Naval University of Engineering, Wuhan 430033)

The phase quality map of phase derivative variance is reliable. When the size of computing window is larger, the computing efficiency is too low, which affects the speed of signal processing. A fast computing method for phase derivative variance is proposed in this paper, which uses multi-computing cores to compute the phase quality map for every blocked wrapped phase. Then the result is downloaded to the memory of host. At last, the high efficiency of the proposed method is verified through experiment.

quality map, phase derivative, parallel computing, GPU

2015年4月5日,

2015年5月17日

李偉,男,工程師,研究方向:信號處理,預(yù)警探測。劉殿剛,男,工程師,研究方向:信號處理,預(yù)警探測。鐘何平,男,博士,講師,研究方向:信號處理,并行計算。

P237

10.3969/j.issn.1672-9730.2015.10.024