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

(國防科技大學氣象海洋學院,南京 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 引 言

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

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

目前,國內外對微波鏈路雨強反演模型普遍采用由國際電信聯盟(ITU)推薦的雨衰冪律模型,該模型是基于數據統計和通信原理建立起的經驗模型[16],其形式為

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

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

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

2.1 微波鏈路雨強反演原理

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

本文基于T矩陣理論和PB降雨粒子模型,采用與我國大部分地區(qū)降雨尺度分布更為相近的Gamma譜分布[19?21],利用OTT雨滴譜儀實測的歷史降雨數據進行非線性擬合,建立微波鏈路雨強反演模型;同時通過測量微波鏈路發(fā)射端與接收端的接收電平差值得到微波路徑傳輸總衰減,根據特征大氣衰減模型修正非降雨因素導致的微波傳輸衰減,從而計算得到鏈路上的微波雨致衰減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為降水造成的微波衰減.其中,降雨造成的衰減是微波傳輸的最主要衰減因素.

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

2.3 非球形雨衰模型

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

2.3.1 PB降雨粒子模型

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

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

Pruppacher和Beard[17]在風洞實驗實測數據的基礎上進行分析,認為降雨粒子的軸比和粒子等效直徑之間進行存在線性關系,提出了PB模型:

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

2.3.2 基于 TT 矩陣理論計算降雨粒子消光截面

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

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

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

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

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

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

2.3.3 非球形雨衰模型

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

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

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

利用(16)式,根據OTT雨滴譜儀記錄的Gamma譜參數歷史數據和基于T矩陣理論計算的降雨粒子消光截面,計算得到微波理論雨致衰減,同時,OTT雨滴譜記錄了相對應的降雨強度.對微波理論雨致衰減和實測雨強采用(1)式的雨衰冪律形式,利用LM最優(yōu)化算法[25](Levenberg-Marquardt algorithm)進行非線性擬合,得到適用于實地氣候雨滴譜特征的雨衰冪律模型.其擬合方法如下.

設定目標函數為

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

此時,擬合問題轉化為求解最小化問題,

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

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

求解的具體步驟如下:

步驟1 設定初始值x0,誤差精度ε=10?5,計算目標函數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,迭代結束,否則進行下一步;

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

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

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

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

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

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

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

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

3 微波鏈路測量實驗及結果分析

不同類型降雨的雨滴形狀是類似的,其差別主要體現在雨滴譜分布上,因此非球形雨衰模型在理論上可適用于不同類型的降雨.為了檢驗非球形雨衰模型反演降雨強度的準確性,本文針對南京地區(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ā)射和接收天線均采用定向天線,信號傳輸為屏蔽電纜,以保證微波信號免受外界干擾,實現微波信號的無損接收.微波發(fā)射機的發(fā)射功率為14 dBm,接收機的采樣間隔為15 s,其分辨率為0.01 dB.

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

為進一步檢驗模型反演效果,本文定量分析了基于非球形雨衰模型的微波鏈路反演的降雨數據,采用雨強時序相關系數、均方根誤差、累計降雨強度絕對偏差和累計降雨相對偏差進行評估,結果列于表2.

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

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

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

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

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

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

4 結 論

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

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

[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大氣科學27 552]

[3]Liang H H,Xu B X,Liu L P,Ge R S 2005Adv.Earth Sci.20 541(in Chinese)[梁海河,徐寶祥,劉黎平,葛潤生2005地球科學進展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大氣科學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紅外與毫米波學報32 379]

[12]Liu X C,Gao T C,Liu L,Zhai D L 2014Acta Phys.Sin.63 199201(in Chinese)[劉西川,高太長,劉磊,翟東力2014物理學報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物理學報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物理學報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氣象學報56 123]

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

[21]Zheng J H,Chen B J 2007J.Meteorol.Sci.27 17(in Chinese)[鄭嬌恒,陳寶君 2007氣象科學 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