李 偉 劉殿剛 鐘何平

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

?

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

李 偉1劉殿剛1鐘何平2

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

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

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

Class Number P237

1 引言

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

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

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

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

(1)

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

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

(2)

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

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

(3)

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

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

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

圖2 計(jì)算流程

GPU為實(shí)現(xiàn)快速計(jì)算相位梯度變化質(zhì)量圖提供了可能,采用GPU計(jì)算相位梯度變化質(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ù),啟動(dòng)相位梯度變化質(zhì)量圖的并行計(jì)算函數(shù),計(jì)算后的結(jié)果存儲(chǔ)在臨時(shí)顯存中;

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

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

4 實(shí)驗(yàn)結(jié)果分析

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

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

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

圖4 加速比

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

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

5 結(jié)語

本文提出了一種基于GPU的快速相位梯度變化質(zhì)量圖計(jì)算方法。首先將干涉圖上傳至顯存,然后采用多個(gè)線程塊同時(shí)計(jì)算每個(gè)分塊干涉相位的質(zhì)量值,并將每個(gè)線程塊計(jì)算的質(zhì)量值進(jìn)行合并。最后將計(jì)算的相位質(zhì)量圖下載到主機(jī)內(nèi)存。試驗(yàn)結(jié)果表明,對(duì)于大小為1024×1024的干涉圖,采用GPU計(jì)算相位梯度變化質(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é)報(bào),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é)報(bào),2007,29(10):2367-2370.

[6] 鐘何平,唐勁松,陳鳴.一種新的InSAR相位質(zhì)量圖生成方法[J].遙感學(xué)報(bào),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日

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

P237

10.3969/j.issn.1672-9730.2015.10.024