鄭 和， 謝亞楠， 萬智龍， 劉文淵

(上海大學(xué)通信與信息工程學(xué)院，上海200072)

基于SAR測量的改進型VIE降雨反演算法

鄭 和， 謝亞楠， 萬智龍， 劉文淵

(上海大學(xué)通信與信息工程學(xué)院，上海200072)

提出一種新型降雨反演算法，該算法是對沃爾塔積分方程(Volterra integral equation，VIE)降雨反演算法的改進．首先，對雷達接收到的回波數(shù)據(jù)進行預(yù)處理，得到降雨起點和降雨寬度等信息;然后，僅利用沃爾塔積分反演算法處理少量特定區(qū)域，通過深度挖掘數(shù)據(jù)變化的規(guī)律來分析得出降雨分布的類型和降雨速率;最后，根據(jù)通常研究對象的分布對稱性，結(jié)合已經(jīng)獲得的信息反演得到整個區(qū)域的降雨分布情況．通過仿真結(jié)果可知，該算法高效、計算速度快，并且具有較高的精度．

反演算法;沃爾塔積分方程改進算法;面向統(tǒng)計反演模型;降雨;合成孔徑雷達

Abstract:This paper proposes a new type of rainfall retrieval algorithm，which improves the method using the Volterra integral equation(VIE)．The proposed algorithm makes pre-analysis on the data received by synthetic aperture radar to get information such as the rainfall start point and the scope of the rain．The VIE retrieval algorithm is then applied in a short distance to get information of the rain distribution．Further analysis is then made on the information obtained in the previous step to find parameters such as shape of the rain and the rain rate．According to the most commonly studied rainfall distributions that are symmetric，the rainfall distribution is obtained．Simulation results show that the proposed algorithm is effective and efficient，with high accuracy．

Key words:inversion algorithm;improved Volterra integral equation(VIE)method;statistic oriented inversion model;rainfall;synthetic aperture radar(SAR)

降雨不僅影響人們的日常生活，而且改變了全球的熱分布，因此對全球范圍內(nèi)的降雨情況進行監(jiān)控，對氣象和水循環(huán)研究都具有重要的意義．常用的降雨監(jiān)測方法主要有2種:一種是使用地基雷達，但是地基雷達的觀測范圍小，無法實現(xiàn)針對全球范圍的有效監(jiān)測;另外一種是采用星載雷達進行監(jiān)測，主要有紅外遙感和微波遙感．目前全球唯一在運行的是1997年美國和日本聯(lián)合研制的“熱帶降雨使命(tropical rainfall measuring mission，TRMM)”衛(wèi)星［1］，這顆衛(wèi)星已經(jīng)運行了十幾年，積累了大量的降雨數(shù)據(jù)．由于TRMM的“壽命”將盡，美國和日本聯(lián)合研制的“全球降雨使命(global precipitation mission，GPM)”衛(wèi)星也將發(fā)射．文獻［2-3］分別對因測雨雷達波束不均勻填充帶來的誤差和后向迭代降雨反演算法的適用性進行了分析．而且，文獻［4］提出了對星載測雨雷達相控陣天線實現(xiàn)低副瓣的粒子群算法．

合成孔徑雷達的回波數(shù)據(jù)含有豐富的降雨區(qū)域和地面背景散射信息，能夠準(zhǔn)確得到地面的降雨分布情況．1999年中國首顆SAR衛(wèi)星開始工程研制，但中國SAR衛(wèi)星都沒有氣象觀測任務(wù)．由于測雨雷達(precipitation radar，PR)造價昂貴、元器件采購困難，而且X波段和Ku波段很接近，因此探索如何將X-SAR用于全球降雨監(jiān)測，對于測雨理論的補充和滿足實際需要都有重要意義．Moore等［5-6］針對用于測雨的合成孔徑雷達，從雷達方程、最小可測雨量等方面進行了研究;Pichugin等［7］提出了VIE反演算法并進行了驗證，該算法數(shù)學(xué)推導(dǎo)嚴謹，但是計算繁雜，運算量很大;Weinamn和 Marzano等［8-9］分別詳細論述了面向模型的反演算法，以及將球面波簡化為平面波帶來的誤差，并把VIE反演算法和面向模型反演算法進行了比較;Marzano等［10］進一步論證了合成孔徑雷達測雨的能力，并且給出了根據(jù)雨區(qū)衰減補償回波圖像的方法;Xie等［11-12］提出了對于面向模型反演算法的改進并通過對雨滴形狀的分析來提高反演精度．本工作在上述研究的基礎(chǔ)上，為了克服VIE反演算法計算復(fù)雜、運算量大的缺點，結(jié)合回波數(shù)據(jù)預(yù)處理、VIE反演算法以及雨區(qū)分布的對稱性，提出了一種新型反演算法，達到了高效、準(zhǔn)確反演降雨量的目的．

1 歸一化雷達反射截面方程

歸一化雷達截面來源于兩個方面:地面散射以及雷達波面遇到降雨粒子帶來的體散射．X-SAR測雨模型如圖1所示．本模型忽略地球曲率，不考慮波束寬度，將球面波簡化為平面波來處理．則歸一化雷達反射截面(normalized radar cross section，NRCS)為［8］

圖1 降雨反演算法原理圖Fig．1 Rainfall retrieval method schematic diagram

式中，

其中σ0∑為歸一化雷達反射截面，σsrf為面散射回波，σvol為體散射回波．如圖1所示，x為與觀測點之間的橫向距離，y為與觀測點之間的縱向距離，方向垂直指向紙張(向里)，Z軸表示高度．(ξ，μ)為體散射波回波路徑上任意一點的坐標(biāo)，(γ，τ)為波面處散射粒子的坐標(biāo)，h為降雨區(qū)的高度，θ為入射角．xL為雨區(qū)開始位置的橫坐標(biāo)，xR為雨區(qū)終止位置的橫坐標(biāo)，xZ為體散射波結(jié)束位置的橫坐標(biāo)．

2 改進的VIE降雨反演算法

由雷達測雨的近似公式［8］可知:

式中，系數(shù)a，b，c，d與雷達波的頻率、水滴的大小和分布、溫度等條件是密切相關(guān)的;α為降雨對雷達波的衰減;I為降雨強度;η為降雨對雷達波束的散射．對于降雨分布，Weinman等［8］提出了幾種典型的系數(shù)取值:a=2．6 ×10－3，b=1．11，c=300，d=1．35．如果把式(4)和(5)代入到式(1)中，可得到一個非線性的積分方程．應(yīng)用VIE反演算法對其進行處理，使得方程簡化為線性的第二類沃爾塔積分方程，從而實現(xiàn)降雨的反演，但是，使用上述方法使計算復(fù)雜度大幅增加．本工作提出的改進的VIE降雨反演算法有效解決了這個問題．

2．1 VIE反演算法推導(dǎo)

引入函數(shù)［7］

其中fx'(x)為f(x)對x的一階導(dǎo)數(shù)．這里沒有考慮y，是因為現(xiàn)在處理的是圖像的同一行數(shù)據(jù)．處理不同行數(shù)據(jù)就可以得到整個分布．考慮到

則將式(4)～(10)代入式(1)，并將觀測點從x移到x+h tan θ，可得［9］

式中，

由式(11)解出f(x)得

情況三：如圖6，作△ADB的外接圓⊙E，假設(shè)E在AB上，連接DE，在⊙E中，∠DEA=2∠DBA=60°，又因為DE=AE，所以△ADE為等邊三角形，所以AD=AE，因為AB=2AE，AC=AB，所以AC=2AE，因為AD=AE=CD，所以AC=AD+CD，在△ACD中AC

通過解第二類沃爾塔積分方程就可以解出式(14)．具體步驟如下:在區(qū)間x＞xZ，通過圖1可知，反射波不經(jīng)過降雨區(qū)，即α=0，代入式(6)易知［9］，

在區(qū)間 xZ－h(huán) tanθ≤x≤xZ，由式(12)和(14)得γmin＞xZ，由式(16)得 f(γ)=1，則 fγ'(γ)=0，則由式(14)易知［9］，

如上所述類推，即可得到f(x)的分布情況;然后采用數(shù)值積分的方法，利用插值得到所有分布;最后代入式(9)即得反演結(jié)果．對于上述反演算法，每隔h tanθ距離，反演的表達式都不相同．而且在計算式(14)時，需要使用數(shù)值積分，比如使用Newton-Cotes公式進行數(shù)值積分計算．由于計算科特斯系數(shù)的限制，進行數(shù)值積分的點數(shù)不宜太多，這樣缺少的信息只能通過插值來獲得，因此帶來了較大的誤差．而且很顯然，隨著降雨區(qū)域長度的增加，VIE反演算法的計算量和計算的復(fù)雜程度一直在大幅增加，這大大降低了其反演的處理能力．

2．2 VIE反演算法的改進

為了解決上述VIE算法存在的問題，本工作提出了改進算法，以達到不僅計算量小而且高效準(zhǔn)確的目的．本算法的核心思想是通過對回波數(shù)據(jù)的預(yù)處理，得到降雨起點、回波最小值、降雨區(qū)域?qū)挾鹊刃畔?然后，通過在一個特定的區(qū)域使用VIE反演算法，獲得降雨分布的形狀和降雨速率等信息．由于通常將降雨區(qū)域視為對稱分布，這樣根據(jù)上面的兩個步驟就可以對稱得出整個降雨區(qū)域的分布．這種方法隨著降雨區(qū)域的增大，計算量增加不多且保證了反演精度．具體步驟如下:

(1)反演降雨起始點x^0和回波最小值x^min，σSAR(x)＜{An(σSAR)－3Rn［σSAR－An(σSAR)］}，(18)式中，An和Rn值分別表示在x點前5個回波數(shù)值的平均值和均方根［4］．當(dāng)x滿足上述條件時，即為降雨起始點x^0．

通過比較回波數(shù)據(jù)就可以得到回波最小值x^min．但是為了避免隨機噪聲帶來的影響導(dǎo)致誤判，需要多考慮兩個參數(shù)，D1和D2．D1為對回波數(shù)據(jù)進行一階求導(dǎo)，D2為對回波數(shù)據(jù)進行二階求導(dǎo)．根據(jù)拐點的特點容易得到，在拐點處D1=0，D2值最大，

(2)反演降雨寬度w^．根據(jù)統(tǒng)計公式，對于矩形分布，

對于三角形分布，對于梯形分布，采用矩形分布和三角形分布的加權(quán)平均來表示，即

(3)反演降雨水平分布的形狀和降雨速率．

如式(16)和(17)所示，本算法在處理降雨區(qū)域時，僅僅在區(qū)間長度為h tanθ這個較小的區(qū)域使用VIE反演算法，然后對其他數(shù)值進行填充．通過分析數(shù)據(jù)變化的規(guī)律，反演降雨水平分布的形狀和降雨速率．接著，代入式(9)中得到區(qū)間(0，xZ)降雨的分布類型，存儲于數(shù)組I(·)中并設(shè)Δx為計算步長．結(jié)合矩形分布、梯形分布和三角形分布的特點易知，若k=ΔI(x)/Δx＞30，則可以認為斜邊斜率足夠大而視為矩形分布．對于矩形分布，斜邊的斜率理論上是無窮大，但是由于數(shù)值運算求導(dǎo)的限制，實際運算出的斜率與網(wǎng)格的劃分有關(guān)，這里選取一個典型值．如果不符合上述條件，對于I(·)中所有的點計算以下兩個參數(shù):

式中，k1和k2意味著I(x)數(shù)值內(nèi)每個點的斜率，I(n)表示數(shù)組I(·)中的任一數(shù)值．根據(jù)梯形分布和三角分布的特點易知，如果出現(xiàn)k1≥2k2且k2≥0，說明是梯形分布;如果k1≥2k2且k2＜0，說明是三角形分布．對數(shù)組I(·)進行處理，記數(shù)組中第一個拐點為I(M)，第一個不為零的點為I(N)，則對于梯形分布，I(M)為反演的降雨量;對于矩形分布，查找I(·)的最大值即為反演的降雨量;對于三角形分布，如果 I(M)和 I(N)之間的距離大于等于w^/2，則I(M)是最大降雨量．否則，根據(jù)斜率相等的原則，填充數(shù)組I(·)，使得I(M)和 I(N)之間的距離等于w^/2，此時I(M)是三角形分布的最大降雨量．

(4)反演整個區(qū)域的降雨水平分布．

如上所述，降雨分布的起點、終點、形狀、大小都已經(jīng)確定，下面分不同情況將降雨分布信息反演出來．對于矩形分布，給定起點、終點和大小即可確定;對于梯形分布，左側(cè)斜邊的長度為L1=(M－N)Δx，則梯形的上底寬度為w1=w^－2L1;對于三角形分布，如果L1＜w^/2，說明左側(cè)三角形并沒有達到降雨速率的最大值．由此依據(jù)三角形左側(cè)圖形斜率相等，將左側(cè)三角形補充完整，然后左右對稱即可．反之，則直接在L1=w^/2處左右對稱．

(5)補償降雨水平分布區(qū)域的位置．

從對I(x)的處理過程可知，反演出的降雨起點并不準(zhǔn)確，但式(18)已經(jīng)得到了降雨的起始點．這樣將整個區(qū)域?qū)?yīng)進行降雨位置調(diào)整即可得到最終的反演結(jié)果．

3 仿真結(jié)果

下面分別對常見的矩形降雨分布、梯形降雨分布和三角形降雨分布進行仿真，驗證算法的可靠性．重點探討對VIE反演算法的改進效果，并對此進行仿真驗證，相關(guān)應(yīng)用會在今后進行進一步研究．

降雨分布類型為矩形分布，降雨起點為8．34 km，降雨寬度為10 km，降雨速率為20 mm/h．反演結(jié)果如下:降雨分布類型為矩形分布，降雨起點為7．902 km，降雨寬度為 10．024 km，降雨速率為19．29 mm/h，降雨速率反演誤差為 3．55%(見圖2)．

降雨分布類型為梯形分布，降雨起點為8．34 km，降雨寬度為10 km，左側(cè)邊斜率為6．7，降雨速率為20 mm/h．反演結(jié)果如下:降雨分布類型為梯形分布，降雨起點為 8．510 km，降雨寬度為9．456 km，降雨速率為17．67 mm/h，左側(cè)邊斜率為5．8，降雨速率反演誤差為11．65%(見圖3)．

降雨分布類型為三角形分布，降雨起點為8．34 km，降雨寬度為10 km，降雨速率為20 mm/h．反演結(jié)果如下:降雨分布類型為三角形分布，降雨起點為8．206 km，降雨寬度為10．614 km，降雨速率為18．73 mm/h，降雨速率反演誤差為6．35%(見圖4)．

仿真結(jié)果分析:上述仿真設(shè)定的雷達入射角為40°，地面最大面散射為－5 dB，是常見的典型取值．仿真實際的降雨起點和降雨寬度以及降雨速率的選取不影響仿真結(jié)果．上述結(jié)果表明:①梯形分布和三角形分布的降雨反演誤差均在10%以內(nèi)．梯形分布由于在數(shù)據(jù)預(yù)處理時無法獲得準(zhǔn)確的統(tǒng)計公式，導(dǎo)致反演誤差稍大，達到11．65%，但是仍能滿足反演算法要求，擬合出梯形分布公式則可以大幅提高梯形分布反演精度．本算法的反演精度與VIE反演算法的精度11%左右［7］相當(dāng);②本算法僅僅計算長度為h tanθ的數(shù)據(jù)，然后通過數(shù)據(jù)預(yù)處理和對稱來得到整體降雨水平分布;而VIE反演算法需要計算(0，xZ)內(nèi)的降雨分布，且每隔h tanθ的長度，都需要重新計算式(14)．所以，本算法要比VIE算法的計算量小得多，且隨著雨區(qū)寬度的增大，計算效率高的優(yōu)勢將更加明顯．

圖2 矩形降雨分布反演結(jié)果Fig．2 Retrieved results of rectangular rainfall distribution

圖3 梯形降雨分布反演結(jié)果Fig．3 Retrieved results of trapezoidal rainfall distribution

圖4 三角形降雨分布反演結(jié)果Fig．4 Retrieved results of triangle rainfall distribution

4 結(jié)束語

本工作分析了VIE降雨反演算法的過程和特點，對算法存在的問題進行了討論，解決了原算法存在的數(shù)學(xué)計算過程復(fù)雜、計算量大等問題;提出了依據(jù)回波預(yù)處理、降雨區(qū)域的對稱性并結(jié)合VIE反演算法來反演降雨的改進反演算法．從仿真結(jié)果可以看出，該算法的反演精度在10%左右，與VIE反演算法精度相當(dāng)．而且從上述理論分析可知，本算法能大大簡化數(shù)學(xué)運算過程，大幅降低運算量．針對本算法的下一步工作是對相關(guān)應(yīng)用情況進行研究．

［1］ MENEGHINI R，JONES J A． Standard deviation of spatially averaged surface cross section data from the TRMM precipitation radar［J］．IEEE Geoscience and Remote Sensing Letters，2011，8(2):293-297．

［2］ DURDEN SL，TANELLI S．Predicted effects of nouniform beam filling on GPM radar data［J］．IEEE Geoscience and Remote Sensing Letters，2008，5(2):308-310．

［3］ SETO S， IGUCHI T． Applicability of the iterative backward retrieval method for the GPM dual-frequency precipitation radar ［J］． IEEE Transactions on Geoscience and Remote Sensing，2011，49(6):1827-1838．

［4］ 徐鋒明，孟令琴，謝亞楠．基于改進的粒子群算法實現(xiàn)階梯幅度量化相控陣天線的低副瓣［J］．上海大學(xué)學(xué)報:自然科學(xué)版，2010，16(4):361-366．

［5］ ATLAS D， MOORE R K． The measurement of precipitation with synthetic aperture radar［J］．Journal of Atmospheric and Oceanic Technology，1987，4(3):368-376．

［6］ MOORE R K，MOGILI A， FANG Y， et al． Rain measurement with SIR-C/X-SAR ［J］．Remote Sensing of Environment，1997，59(2):280-293．

［7］ PICHUGIN A P，SPIRIDONOV Y G．Spatial distributions of rainfall intensity recovery from space radar images［J］．Sov J Remote Sens，1991，8:917-932．

［8］ WEINMAN JA，MARZANO F S．An exploratory study to derive precipitation over land from X-band synthetic aperture radar measurements［J］．Journal of Applied Meteorology and Climatology，2008，47:562-575．

［9］ MARZANO F S，WEINMAN J A．Inversion of spaceborne X-band synthetic aperture radar measurements for precipitation remote sensing over land ［J］．IEEE Transactions on Geoscience and Remote Sensing，2008，46(11):3472-3487．

［10］ MARZANO F S，MORI S，WEINMAN J A．Evidence of rainfall signatures on X-band synthetic aperture radar imagery overland［J］．IEEE Transactions on Geoscience and Remote Sensing，2010，48(2):950-964．

［11］ XIE Y N，HUAN J P，TAO Y．Study on the algorithm to retrieve precipitation with X-band synthetic aperture radar［J］．Acta Meteorologica Sinica，2010，24(5):614-621．

［12］ XIE Y N，ZHOU X L，YANG Z D．Error analysis of nonspherical raindrops on precipitation measurement［J］．Journal of Shanghai University:English Edition，2011，15(2):92-95．

Improved VIE Rainfall Retrieval Algorithm Based on SAR Measurements

ZHENG He， XIE Ya-nan， WAN Zhi-long， LIU Wen-yuan
(School of Communication and Information Engineering，Shanghai University，Shanghai 200072，China)

P 435

A

1007-2861(2012)05-0464-06

10．3969/j．issn．1007-2861．2012．05．005

2011-11-09

國家自然科學(xué)基金資助項目(61071185);上海市重點學(xué)科建設(shè)資助項目(S30108);上海市科委重點實驗室資助項目(08DZ2231100);上海市科委科技攻關(guān)計劃資助項目(10511501702)

謝亞楠(1962～)，男，研究員，博士生導(dǎo)師，博士，研究方向為微波遙感、電磁散射、射頻電路等．E-mail:yxie@shu．edu．cn