基于非球形雨衰模型的微波鏈路雨強反演方法?

宋堃 高太長 劉西川 印敏 薛楊

(國防科技大學(xué)氣象海洋學(xué)院,南京 211101)

Song Kun Gao Tai-Chang?Liu Xi-Chuan Yin Min Xue Yang

(College of Meteorology and Oceanography,National University of Defense Technology,Nanjing 211101,China)

1 引 言

降雨是十分重要的物理現(xiàn)象,實時準(zhǔn)確地測量降雨對于人民生活、社會生產(chǎn)、部隊保障和國家安全具有十分重要的意義[1].目前,在氣象觀測業(yè)務(wù)中測量降雨的手段主要有雨量計、天氣雷達和氣象衛(wèi)星[2?4],地面雨量計能夠較為準(zhǔn)確地測量地表降雨強度,然而其空間代表性較差,無法監(jiān)測降雨的精細分布;天氣雷達雖能提供時空分辨率高的雨強數(shù)據(jù),但受雷達掃描方式和Z-R關(guān)系不確定性(Z為測雨雷達反射率,R為降雨強度)的制約,其測量結(jié)果與地表真實降雨存在一定差異;氣象衛(wèi)星可自上而下測量云頂或穿透云頂,但其探測目標(biāo)仍不是地表降雨,測雨精度仍有待提高.現(xiàn)階段,實時準(zhǔn)確地測量地表降水仍然存在一定的難度.

近些年,氣象學(xué)者提出利用目前廣泛存在的微波鏈路通信信號傳輸過程中的衰減信息反演地表降雨強度的設(shè)想,并已開展了大量的實驗研究,取得了一定的進展.2006年,Messer等[5]提出利用商用微波鏈路探測降雨的設(shè)想,并提出了基于雨強-微波衰減冪律關(guān)系的雨衰反演模型.2008年,Feng等[6]利用啟發(fā)式分割算法(brief l y BG algorithm)對中國南北方的降雨時序突變開展研究,在此基礎(chǔ)上分析了不同地區(qū)的降雨冪律衰減特性.2009年,Goldshtein等[7]開展了不同頻段商用微波鏈路降雨測量的實驗,驗證了多頻段微波雨強反演的可行性.2012年,Messer等[8]提出了利用微波鏈路探測雪、冰水混合物等其他相態(tài)降水的設(shè)想,并提出了相關(guān)的反演模型.2013年,David等[9]進行了不同頻段微波鏈路測雨實驗,通過與實際地表降雨進行比對,驗證了微波鏈路測雨的準(zhǔn)確性與可行性;同年,Overeem等[10]利用超過2400條微波鏈路對荷蘭全境進行了多頻段蜂窩無線網(wǎng)絡(luò)實時監(jiān)測降雨的試驗,驗證了基于組網(wǎng)微波鏈路進行區(qū)域雨強反演的可行性.目前,國內(nèi)在微波測雨方面研究較少,劉西川等[11,12]研究了降雨粒子對不同波段微波傳輸特性的影響特性,提高了復(fù)雜降雨條件下微波傳輸衰減的評估精度,并建立了考慮溫度的降水衰減模型;2013年,姜世泰等[13]研究了基于微波鏈路的區(qū)域降雨反演方法,并進行了仿真研究;2015年,印敏等[14]設(shè)計了視距微波通信鏈路系統(tǒng),并驗證了利用微波鏈路開展實時測量降雨的可行性;同年,高太長等[15]提出了基于實地雨滴譜分布的微波鏈路路徑雨強反演算法,并開展了15—20 GHz頻段的微波鏈路測量實驗,驗證了該方法的有效性.

目前,國內(nèi)外對微波鏈路雨強反演模型普遍采用由國際電信聯(lián)盟(ITU)推薦的雨衰冪律模型,該模型是基于數(shù)據(jù)統(tǒng)計和通信原理建立起的經(jīng)驗?zāi)Ｐ蚚16],其形式為

式中,A(dB/km)代表微波傳輸雨致衰減,R(mm/h)代表降雨強度,a和b為雨衰冪律系數(shù).對(1)式進行如下變形即可反解得到降雨強度:

冪律系數(shù)a,b與雨滴形狀、雨滴譜分布等因素有關(guān),為了得到具有實地譜分布的冪律系數(shù),本文利用Pruppacher-Beard(PB)降雨粒子模型[17]和T矩陣?yán)碚揫18]計算非球形粒子消光截面,結(jié)合OTT雨滴譜儀記錄的南京地區(qū)Gamma譜參數(shù)歷史數(shù)據(jù)計算得到不同頻率下微波的雨致衰減,對微波雨致衰減及其對應(yīng)的降雨強度進行(1)式形式的非線性擬合,擬合得到不同頻率下非球形雨衰模型冪律系數(shù).為驗證非球形雨衰模型反演降水強度的精度,開展15—20 GHz頻段微波反演降雨強度的實驗進行檢驗.

2 基于非球形雨衰模型的微波鏈路雨強反演模型

2.1 微波鏈路雨強反演原理

微波在空氣中傳播時會受到其中的氣體、降水等因素的影響而發(fā)生吸收、衰減、散射等能量衰減現(xiàn)象,其中降水是主要影響因素.微波雨致衰減主要取決于降雨強度、降雨粒子形狀、雨滴尺度分布等因素.基于此,根據(jù)微波雨致衰減的大小可以反演實際降雨強度.

本文基于T矩陣?yán)碚摵蚉B降雨粒子模型,采用與我國大部分地區(qū)降雨尺度分布更為相近的Gamma譜分布[19?21],利用OTT雨滴譜儀實測的歷史降雨數(shù)據(jù)進行非線性擬合,建立微波鏈路雨強反演模型;同時通過測量微波鏈路發(fā)射端與接收端的接收電平差值得到微波路徑傳輸總衰減,根據(jù)特征大氣衰減模型修正非降雨因素導(dǎo)致的微波傳輸衰減,從而計算得到鏈路上的微波雨致衰減Arain;通過反演模型和微波雨致衰減就可反演得到實際降雨強度.具體反演流程如圖1所示.

圖1 雨衰反演流程圖Fig.1.Flow chart of the rainfall inversion model.

2.2 微波鏈路實測雨致衰減值提取

微波路徑傳輸總衰減可表示為[22]

式中,A為微波路徑傳輸總衰減,f為微波發(fā)射頻率(GHz),D為微波傳輸總距離(km),Awater為水汽造成的微波傳輸衰減,Afog為霧或輕霧造成的微波衰減,AO2為氧氣造成的微波衰減,Aair為其他氣體吸收造成的衰減,Arain為降水造成的微波衰減.其中,降雨造成的衰減是微波傳輸?shù)淖钪饕p因素.

為提取出微波傳輸雨致衰減,本文在晴空條件下測得微波路徑傳輸衰減基值,設(shè)基值為降雨時實測得到微波鏈路路徑傳輸總衰減Atot.由于降雨現(xiàn)象發(fā)生時,大氣中的溫度、氣壓、濕度等條件和測量晴空基值時的條件會發(fā)生變化,因此要對傳輸衰減進行溫度、氣壓及濕度的修正,根據(jù)降雨前和降雨時的溫度、氣壓、濕度數(shù)據(jù)變化,結(jié)合特征大氣吸收衰減模型進行衰減修正,計算得到降雨過程中微波的衰減修正值A(chǔ)Modi,因此僅由降雨因素造成的微波傳輸衰減Arain可由(4)式計算:

2.3 非球形雨衰模型

降雨粒子的形狀取決于粒子的尺度,一般而言,降雨粒子的直徑在10μm—10 mm之間,大于7 mm的降雨粒子是不穩(wěn)定的,容易發(fā)生破裂.研究表明,等效半徑小于1 mm的降雨粒子基本為球狀粒子;等效半徑達到1.5 mm時,降雨粒子的形狀近似為一個扁旋轉(zhuǎn)橢球體;當(dāng)?shù)刃О霃竭_到2 mm時,降雨粒子的形狀變成一個底部扁平的扁旋轉(zhuǎn)橢球體;當(dāng)?shù)刃О霃竭_到2.75 mm時,橢球體形狀的降雨粒子底部形成一個凹槽,但是這種尺度的降雨粒子在我國很少出現(xiàn).

2.3.1 PB降雨粒子模型

降雨粒子的形狀通常用粒子的軸比來表示,軸比為橢球體短軸與長軸之比:

式中,a為橢球體短軸,b為長軸.

Pruppacher和Beard[17]在風(fēng)洞實驗實測數(shù)據(jù)的基礎(chǔ)上進行分析,認為降雨粒子的軸比和粒子等效直徑之間進行存在線性關(guān)系,提出了PB模型:

式中,Deq為降雨粒子等效直徑.

2.3.2 基于 TT 矩陣?yán)碚撚嬎憬涤炅Ｗ酉饨孛?/h4>

T矩陣?yán)碚?extended boundary condition method)是對任意形狀粒子的Maxwell方程的嚴(yán)格解析解[23].假設(shè)電磁波入射到空間粒子上:

對入射場、散射場按矢量球面波展開如下:

式中,k1為周圍介質(zhì)的波數(shù);RgMmn,RgNmn,Mmn和Nmn是矢量球諧函數(shù).入射場的展開系數(shù)amn,bmn表示如下:

式中,C?mn,B?mn為球矢面變量.由于Maxwell方程及邊界條件是線性的,散射場的展開系數(shù)pmn和qmn與入射展開系數(shù)amn和bmn之間也是線性的.利用假設(shè)的轉(zhuǎn)換矩陣T給出如下關(guān)系:

對(12)和(13)式采用矩陣形式改寫:

(14)式表明散射場的展開系數(shù)由T矩陣和入射場展開系數(shù)相乘獲得,基于(14)式可計算得到粒子消光截面:

2.3.3 非球形雨衰模型

基于上述的T矩陣?yán)碚?計算得到PB降雨粒子模型的消光截面,根據(jù)雨致衰減理論公式計算理論雨致衰減[24],

式中,N(D)為Gamma雨滴譜分布(m?3·mm?1),其譜分布函數(shù)為

式中,N0為粒子數(shù)密度(m?3·mm?1+μ),Λ為尺度因子(mm?1),μ為形狀因子,為無量綱量.

利用(16)式,根據(jù)OTT雨滴譜儀記錄的Gamma譜參數(shù)歷史數(shù)據(jù)和基于T矩陣?yán)碚撚嬎愕慕涤炅Ｗ酉饨孛?計算得到微波理論雨致衰減,同時,OTT雨滴譜記錄了相對應(yīng)的降雨強度.對微波理論雨致衰減和實測雨強采用(1)式的雨衰冪律形式,利用LM最優(yōu)化算法[25](Levenberg-Marquardt algorithm)進行非線性擬合,得到適用于實地氣候雨滴譜特征的雨衰冪律模型.其擬合方法如下.

設(shè)定目標(biāo)函數(shù)為

式中,Ri為實測降雨強度(mm/h),Ai為理論雨致衰減(dB/km).(18)式向量形式如下:

此時,擬合問題轉(zhuǎn)化為求解最小化問題,

式中,x=[a,b],xopt=[aopt,bopt].

圖2 (網(wǎng)刊彩色)微波雨衰模型擬合結(jié)果Fig.2.(color online)The f i tting results of rain-attenuation relationship.

求解的具體步驟如下:

步驟1 設(shè)定初始值x0,誤差精度ε=10?5,計算目標(biāo)函數(shù)Fopt(x0);

步驟2 令μM=0.001,k=0;

步驟3 求解δxk,計算Fopt(xk+1),δxk的計算公式為

式中,Jk為雅可比(Jacobi)矩陣在xk點的值,I為單位矩陣,xk+1的計算公式為xk+1=xk+δxk;

步驟4 如果μM=0,迭代結(jié)束,否則進行下一步;

步驟5 如果Fopt(xk+1)＞Fopt(xk),則令xk+1=xk+δxk,μM=10μM,并轉(zhuǎn)入步驟3;否則,若δxk＜ε,則結(jié)束迭代,若δxk≥ε,令xk+1=xk+δxk,μM=0.1μM,轉(zhuǎn)入步驟3.

圖2為15—20 GHz基于非球形雨衰模型的擬合結(jié)果,其中OTT雨滴譜儀歷史數(shù)據(jù)為2013年6月的南京地區(qū)Gamma雨滴譜和降雨強度的實測數(shù)據(jù),對實測數(shù)據(jù)按分鐘平均,得到4620個數(shù)據(jù)樣本點.從圖2可以看出,微波頻率在15—20 GHz范圍時,雨衰模型的冪律系數(shù)b接近1,說明在該頻段范圍,微波雨致衰減與降雨強度呈近似線性關(guān)系.

為定量評價擬合效果,采用擬合值與理論值間的均方根誤差(root mean square error,RMSE)、標(biāo)準(zhǔn)化平均誤差(normalized mean error,NME)和標(biāo)準(zhǔn)化平均偏差(normalized mean bias,NMB)對擬合結(jié)果進行評估.其中,RMSE反映的是擬合值與理論值之間的誤差大小,NME和NMB反映的是擬合值與理論值之間的相對偏差,若NME和NMB的值均小于20%,表示擬合效果較好.從表1中可以看出,各頻率下RMSE值均在0.04—0.05之間,NME值均在15%以下,NMB值在5.5%以下.評估結(jié)果表明,模型擬合值與理論值間相差較小,擬合誤差較小,擬合精度較高.

圖3分別給出了15—20 GHz頻率范圍內(nèi),微波在不同溫度條件下非球形雨衰模型的冪律系數(shù)a,b值.由表可知,同一頻率時雨衰模型參數(shù)a,b隨溫度的變化比較小,隨著溫度升高,a值逐漸變大.當(dāng)微波頻率為15 GHz時,b值隨溫度的變化不明顯,在16—19 GHz時,b值隨溫度的升高呈下降趨勢.同時可以發(fā)現(xiàn),微波頻率越高,a值呈上升趨勢,b值呈下降趨勢.

圖4是在溫度為25°C時,冪律關(guān)系中的b值隨頻率的變化情況.從圖4可以看出,頻率越高,b值越小.微波頻率在32 GHz左右時,b值接近1,說明在此頻率下微波雨致衰減與降雨強度接近呈線性關(guān)系.

表1 微波雨衰模型擬合結(jié)果Table 1.The f i tting results of rain-attenuation relationship.

圖3 (網(wǎng)刊彩色)非球形雨衰模型冪律系數(shù)Fig.3.(color online)Power-law relationship for different temperature and frequency.

圖4 (網(wǎng)刊彩色)非球形雨衰模型b值 (25°C)Fig.4.(color online)b in power-law relationship for different frequency(25 °C).

3 微波鏈路測量實驗及結(jié)果分析

不同類型降雨的雨滴形狀是類似的,其差別主要體現(xiàn)在雨滴譜分布上,因此非球形雨衰模型在理論上可適用于不同類型的降雨.為了檢驗非球形雨衰模型反演降雨強度的準(zhǔn)確性,本文針對南京地區(qū)一次層狀云降雨事件,開展了15—20 GHz頻段的視距微波鏈路反演路徑降雨實驗,進行個例分析,并與OTTPARSIVEL激光粒子譜儀進行同步對比.本次實驗時間段為2013年7月5日13:00—18:00及7月6日11:00—14:00,實驗區(qū)域介于南京市江寧區(qū)和玄武區(qū)之間,整個實驗覆蓋了從暴雨、大雨、中雨、小雨到毛毛雨的全部過程.實驗鏈路長度為6.57 km,微波發(fā)射機型號為Anritsu MG3694B信號發(fā)生器,信號接收機型號為Agilent E4440A頻譜分析儀,微波發(fā)射和接收天線均采用定向天線,信號傳輸為屏蔽電纜,以保證微波信號免受外界干擾,實現(xiàn)微波信號的無損接收.微波發(fā)射機的發(fā)射功率為14 dBm,接收機的采樣間隔為15 s,其分辨率為0.01 dB.

圖5和圖6分別給出了2013年7月5日和7月6日基于非球形雨衰模型反演結(jié)果與OTT雨滴譜儀測量值的對比.從實驗結(jié)果可以看出,基于非球形雨衰模型反演的降雨強度與OTT雨滴譜儀觀測的降雨強度之間均存在較好的相關(guān)性.對于降雨強度和累積降雨量,與OTT雨滴譜儀的實測雨強相比,非球形雨衰模型反演的雨強偏大.需要強調(diào)的是,降雨在空間分布是不均勻的,OTT雨滴譜儀只能測量單點降雨情況,而微波鏈路反演的是微波傳輸路徑的平均雨強.因此微波鏈路反演的雨強值與OTT雨滴譜儀觀測的降雨強度無法完全一致.

為進一步檢驗?zāi)Ｐ头囱菪Ч?本文定量分析了基于非球形雨衰模型的微波鏈路反演的降雨數(shù)據(jù),采用雨強時序相關(guān)系數(shù)、均方根誤差、累計降雨強度絕對偏差和累計降雨相對偏差進行評估,結(jié)果列于表2.

圖5 (網(wǎng)刊彩色)微波鏈路反演降雨結(jié)果(2013年7月5日)Fig.5.(color online)Inversion results of microwave link(2013.07.05).

圖6 (網(wǎng)刊彩色)微波鏈路反演降雨結(jié)果(2013年7月6日)Fig.6.(color online)Inversion results of microwave link(2013.07.06).

表2 微波鏈路降雨反演結(jié)果Table 2.Inversion results of microwave link.

由表2可知,基于非球形雨衰模型的微波鏈路降雨強度時序反演值與OTT雨滴譜儀實測的降雨強度值存在一定的相關(guān)性,相關(guān)系數(shù)都在0.6以上,最高達到了0.96,說明基于微波鏈路雨致衰減反演降雨強度具有一定的可行性;對于鏈路反演雨強的RMSE,其結(jié)果在0.79—3.67 mm/h之間.從累積降雨的絕對誤差來看,結(jié)果在0.28—2.47 mm之間,說明微波鏈路反演的累積降雨量略高于OTT實測降雨強度;從相對偏差來看,反演結(jié)果在0.44%—1.84%之間,說明累積降雨量的反演結(jié)果較為理想.同時,從表2還可以看出,第一組實驗的相關(guān)系數(shù)為0.65,低于其他三組,然而其對應(yīng)雨強RMSE、累積降雨絕對偏差卻低于其他三組試驗的結(jié)果,原因主要為第一組實驗降雨量較小,所以其反演值與實測值的絕對差值要小于其他三組,因而該組雨強RMSE、累積降雨絕對偏差小于另外三組,但從累積降雨相對偏差看,第一組實驗大于其他三組,進而導(dǎo)致其相關(guān)系數(shù)較其他三組而言相對較低.

從表2還可以看出,非球形雨衰模型對小雨的反演精度低于中雨和大雨,其原因主要為降雨強度較小時,微波的雨致衰減值也較低,由于無法完全消除非雨致因素造成的微波衰減,因此降雨較弱時,非雨致微波衰減對反演精度的影響要大于強降雨時對反演精度的影響.

綜上所述,本文提出的基于非球形雨衰模型的微波鏈路反演雨強方法的具有一定的可行性,基于T矩陣?yán)碚摵蚉B降雨粒子模型,對OTT雨滴譜儀的Gamma譜參數(shù)歷史資料與降雨強度值進行非線性擬合得到具有實地譜分布的冪律系數(shù),建立適合于本地區(qū)的雨衰模型,便于進一步開展微波鏈路雨強反演研究.

4 結(jié) 論

本文以降雨粒子對電磁波的衰減效應(yīng)為物理原理,以T矩陣?yán)碚?、Gamma譜分布為理論基礎(chǔ),基于PB降雨粒子模型,提出了基于非球形雨衰模型的微波鏈路雨強反演方法,并對南京地區(qū)一次層狀云降雨開展了15—20 GHz頻段的視距微波鏈路反演路徑降雨實驗,對非球形雨衰模型進行驗證.實驗結(jié)果表明,反演雨強的相關(guān)系數(shù)全部高于0.6,最高達到了0.96,RMSE最小值為0.79,累積降雨量的絕對偏差在2.47 mm以內(nèi),最小偏差僅為0.28 mm,相對誤差低于1.84%.本次實驗結(jié)果驗證了本文提出的基于非球形雨衰模型的微波鏈路雨強反演方法的準(zhǔn)確性和可行性,對于進一步提高微波鏈路反演降雨精度、改善降水監(jiān)測效果具有重要意義.

同時,實驗結(jié)果也表明,非球形雨衰模型對小雨的反演精度低于中雨和大雨,主要是受非雨致因素造成的微波衰減的影響,降雨強度較小時,微波雨致衰減值較低,非雨致微波衰減對反演精度的影響要大于強降雨時對反演精度的影響,從而影響了雨強反演精度.因此,需要進一步開展消除非雨致微波衰減的工作.

[1]Gao T C 2012Meteorol.Hydrol.Eq.23 1(in Chinese)[高太長2012氣象水文裝備23 1]

[2]Lü D R,Wang P C,Qiu J H,Tao S Y 2003Chin.J.Atmos.Sci.27 552(in Chinese)[呂達仁,王普才,邱金恒,陶詩言2003大氣科學(xué)27 552]

[3]Liang H H,Xu B X,Liu L P,Ge R S 2005Adv.Earth Sci.20 541(in Chinese)[梁海河,徐寶祥,劉黎平,葛潤生2005地球科學(xué)進展20 541]

[4]Qie X S,Lü D R,Chen H B,Wang P C,Duan S,Zhang W X 2008Chin.J.Atmos.Sci.32 867(in Chinese)[郄秀書,呂達仁,陳洪濱,王普才,段樹,章文星 2008大氣科學(xué)32 867]

[5]Messer H,Zinevich A,Alpert P 2006Science312 713

[6]Feng G L,Gong Z Q,Zhi R,Zhang D Q 2008Chin.Phys.B17 2745

[7]Goldshtein O,Messer H,Zinevich A 2009IEEE Trans.Signal Proces.57 1616

[8]Messer H,Zinevich A,Alpert P 2012IEEE Trans.Instrum.Meas.15 32

[9]David N,Alpert P,Messer H 2013Atmos.Res.131 13

[10]Overeem A,Leijnse H,Uijlenhoet R 2013Proc.Natl.Acad.Sci.USA110 2741

[11]Liu X C,Liu L,Gao T C,Ren J P 2013J.Infrared Millim.Waves32 379(in Chinese)[劉西川,劉磊,高太長,任景鵬2013紅外與毫米波學(xué)報32 379]

[12]Liu X C,Gao T C,Liu L,Zhai D L 2014Acta Phys.Sin.63 199201(in Chinese)[劉西川,高太長,劉磊,翟東力2014物理學(xué)報63 199201]

[13]Jiang S T,Gao T C,Liu X C,Liu L,Liu Z T 2013Acta Phys.Sin.62 154303(in Chinese)[姜世泰,高太長,劉西川,劉磊,劉志田2013物理學(xué)報62 154303]

[14]Yin M,Jiang S T,Gao T C,Liu X C,Liang M Y,Ge S R,Cao C K 2015Meteorol.Sci.Technol.43 1(in Chinese)[印敏,姜世泰,高太長,劉西川,梁妙元,戈書睿,曹承堃2015氣象科技43 1]

[15]Gao T C,Song K,Liu X C,Yin M,Liu L,Jiang S T 2015Acta Phys.Sin.64 174301(in Chinese)[高太長,宋堃,劉西川,印敏,劉磊,姜世泰2015物理學(xué)報64 174301]

[16]International Telecommunication Union 2005Rec.ITURP.838-3

[17]Pruppacher H R,Beard K V 1970J.Quart.J.R.Met.Soc.96 247

[18]Michael I M,Larry D T 1998J.Quant.Spectrosc.Radiat.Transfer60 309

[19]Chen B J,Li Z H,Liu J C,Gong F J 1998Acta Meteorol.Sin.56 123(in Chinese)[陳寶君,李子華,劉吉成,宮福久1998氣象學(xué)報56 123]

[20]Yuan C,Fan L,Li Y B 2001J.Nanjing I.Meteorol.24 250(in Chinese)[袁成,樊玲,李亞濱 2001南京氣象學(xué)院學(xué)報24 250]

[21]Zheng J H,Chen B J 2007J.Meteorol.Sci.27 17(in Chinese)[鄭嬌恒,陳寶君 2007氣象科學(xué) 27 17]

[22]FreemanR1991TelecommunicationsTransmission Handbook(3rd Ed.)(Canda:John Wiley&Sons Inc.)p279

[23]Somerville W R C,Auguie B,Ru E C L 2013J.Quant.Spectrosc.Radiat.Transfer123 153

[24]David C H,Chu T S 1975Proc.IEEE63 1308

[25]Lampton M 1997Comput.Phys.11 110