張衛(wèi)國,范仲麗,鐘 偉,江雨田,孫飛飛,陳 娟
(1.寧波市水利水電規(guī)劃設(shè)計研究院,浙江 寧波 315192;2.河海大學(xué)水文水資源學(xué)院,南京 210098)
數(shù)值預(yù)報對強(qiáng)對流天氣的臨近預(yù)報技術(shù)還不夠成熟,但強(qiáng)對流天氣引發(fā)的氣象災(zāi)害對民眾的生命財產(chǎn)安全造成極大的危害。因此,對強(qiáng)對流天氣的預(yù)警和臨近預(yù)報意義重大。而雷達(dá)估測降雨在城市監(jiān)測預(yù)警方面的應(yīng)用將是未來水利防洪領(lǐng)域的前沿[1]。
天氣雷達(dá)具有高時空分辨率的優(yōu)勢,可以快速、大范圍捕捉強(qiáng)對流天氣情況,也因此越來越多地被應(yīng)用到強(qiáng)對流天氣短時臨近降雨預(yù)報中。雷達(dá)回波最優(yōu)空間相關(guān)方法最早被用于Austin和Bellon[2]研究中,利用平均移動矢量對降雨進(jìn)行了預(yù)報;在此基礎(chǔ)上,Rinehart和Gravey[3]發(fā)展了TREC(Tracking Radar Echo By Correction)方法,將雷達(dá)回波分為多個區(qū)域,得到不同的移動矢量;Li[4]等在TREC基礎(chǔ)上發(fā)展了COTREC方法,并進(jìn)行了降水預(yù)報。而多數(shù)利用TREC方法的研究中,只作雷達(dá)回波位置和大概形狀的外推,對回波內(nèi)部的量變不作修正。
雨強(qiáng)的臨近估計主要根據(jù)雷達(dá)回波反射率因子Z和雨強(qiáng)R的關(guān)系式進(jìn)行。在多數(shù)研究中,該關(guān)系式中參數(shù)都采用固定參數(shù)或經(jīng)驗參數(shù)。但實際上,不同地區(qū)或者不同降雨過程中,參數(shù)值存在較大差別,使得雷達(dá)降雨估測與實際降雨偏差較大。因此,利用實測雨量數(shù)據(jù)對雨強(qiáng)關(guān)系式中參數(shù)的實時校準(zhǔn),以減少更多的估算誤差是很有必要的。
本文以臺風(fēng)多發(fā)的寧波市為例,利用給定研究區(qū)的雷達(dá)回波最優(yōu)空間相關(guān)方法,得到前一時次每6 min的雷達(dá)回波移動矢量,對該矢量的2個分量分別利用最小二乘回歸得到其與時間的關(guān)系,由此可得下一時次的移動矢量;疊加上一時次雷達(dá)回波內(nèi)部變化量后,得到雷達(dá)回波預(yù)報圖像,實現(xiàn)對雷達(dá)回波位置、形狀的外推及內(nèi)部回波量的修正。對Z~R關(guān)系式進(jìn)行變形,采用SCE-UA算法對雷達(dá)回波強(qiáng)度N與降雨強(qiáng)度R關(guān)系中參數(shù)進(jìn)行估算與實時校準(zhǔn),嘗試預(yù)測下一時次的降水分布情況。
本文采用給定研究區(qū)的雷達(dá)回波最優(yōu)空間相關(guān)方法對強(qiáng)對流天氣下雷達(dá)回波的位置和形狀進(jìn)行外推,并把上一時次的雷達(dá)回波變化網(wǎng)格疊加到下一時次的雷達(dá)回波預(yù)報中,得到最終的雷達(dá)預(yù)報圖像。根據(jù)雷達(dá)反射率因子與雨強(qiáng)的關(guān)系式的變形公式對雨強(qiáng)進(jìn)行估算,利用SCE-UA算法對雨強(qiáng)關(guān)系式中參數(shù)進(jìn)行實時校準(zhǔn),得到時效性較好的A,b參數(shù),從而得到相對準(zhǔn)確的下一時次預(yù)報降雨分布。
從雷達(dá)回波圖像上可以看出,雷達(dá)探測回波存在明顯的、較為明確的移動軌跡(見圖1)。根據(jù)已發(fā)生的短時間內(nèi)的雷達(dá)回波移動軌跡,推測下一時次雷達(dá)回波的位置及形狀,即是雷達(dá)回波最優(yōu)空間相關(guān)方法的核心思想。
雷達(dá)回波最優(yōu)空間相關(guān)方法,是在假設(shè)大氣中水汽的移動受環(huán)境風(fēng)場的引導(dǎo),而環(huán)境風(fēng)場在短時間內(nèi)不存在嚴(yán)重突變的前提下成立的。將已探測到的雷達(dá)回波的移動矢量作為當(dāng)前環(huán)境風(fēng)的移動矢量,來預(yù)測未來短時間內(nèi)雷達(dá)回波的移動矢量。該方法簡便,運算量小,運行時間短,在有大面積層狀云時跟蹤效果較好[5]。
圖1 雷達(dá)回波移動軌跡示意圖(2018年4月4日)Fig.1 The schematic of radar echo movement track (April 4, 2018)
給定研究區(qū)的雷達(dá)回波最優(yōu)空間相關(guān)方法具體步驟是,將整個回波區(qū)域按照研究區(qū)的網(wǎng)格大小劃分成若干個相同大小的矩形網(wǎng)格,研究區(qū)雷達(dá)回波矩形網(wǎng)格與上一時刻的若干矩形網(wǎng)格之間計算相關(guān)系數(shù),尋得的最大相關(guān)系數(shù)的網(wǎng)格到研究區(qū)網(wǎng)格對應(yīng)的矢量即為該天氣系統(tǒng)下相鄰時刻雷達(dá)回波的位移矢量。
相關(guān)系數(shù)R的計算公式如下:
(1)
尋找與研究區(qū)網(wǎng)格最大相關(guān)系數(shù)的網(wǎng)格時,為縮減計算矩形網(wǎng)格的數(shù)量,可根據(jù)雨帶移動的最大期望速度與分析時間間隔得到最大網(wǎng)格搜索半徑。為減小網(wǎng)格距離計算量,可采用最大網(wǎng)格搜索半徑為邊長的正方形作為搜索區(qū)域。按網(wǎng)格大小為一個移動單位進(jìn)行移動,求得在最大搜索半徑范圍內(nèi)的所有網(wǎng)格(計算區(qū))與下一相鄰時次研究區(qū)網(wǎng)格的相關(guān)系數(shù)。其中相關(guān)系數(shù)最大的計算區(qū)網(wǎng)格中心到研究區(qū)網(wǎng)格中心的矢量即為t1到t2時刻的位移矢量(見圖2),并保存該計算區(qū)網(wǎng)格與研究區(qū)網(wǎng)格的變化量網(wǎng)格E。重復(fù)尋找當(dāng)前時刻之前一定時間段內(nèi)相鄰時次的雷達(dá)回波圖像的位移矢量(Δx1,Δy1),(Δx2,Δy2),…,(Δxn,Δyn),記錄一系列變化量網(wǎng)格(E1,E2,…,En)。假設(shè)環(huán)境風(fēng)場方向及風(fēng)速在一定時間段內(nèi)穩(wěn)定,可以分別建立X方向與Y方向上的位移分量Δx、Δy與絕對時間t之間的線性關(guān)系式:
Δx=a1t+b1
(2)
Δy=a2t+b2
(3)
根據(jù)至少2組數(shù)據(jù)[t1,(Δx1,Δy1)],[t2,(Δx2,Δy2)]即可確定a1、b1、a2、b24個參數(shù),理論上可用來計算未來任意時間與其相鄰前一時次間的雷達(dá)回波位移矢量。在實際運用時,需要考慮方法的時效性。
圖2 給定研究區(qū)的雷達(dá)回波最優(yōu)空間相關(guān)方法Fig.2 Radar echo optimal spatial correlation method for a given study area
在確定雷達(dá)回波外推位置與形狀后,疊加上一時段的變化量網(wǎng)格,可以得到最終的外推雷達(dá)回波圖像。
復(fù)合形交叉進(jìn)化算法(SCE-UA)是一種全局優(yōu)化算法,集成了隨機(jī)搜索算法、單純形法、聚類分析法等方法的優(yōu)點,是一種有效解決非線性約束最優(yōu)化問題的混合算法[6,7]。該算法引入了種群的概念。復(fù)合形點在可行域內(nèi)隨機(jī)生成、競爭演化,在多個吸引域內(nèi)獲得全局收斂點,能有效地表達(dá)參數(shù)敏感性與參數(shù)之間的相關(guān)性。SCE-UA算法基本思路是首先在參數(shù)可行域內(nèi)引入隨機(jī)分布的點群,將其分為多個復(fù)合形,每個復(fù)合形中有2n+1個點(其中n是需要優(yōu)化的參數(shù)個數(shù));每個復(fù)合形根據(jù)下降單純形算法進(jìn)行進(jìn)化;定時地將整個點群混合進(jìn)而形成新的復(fù)合形,也就將之前的復(fù)合形所包含的信息引入新的復(fù)合形中;進(jìn)化與混合不斷進(jìn)行直至滿足收斂準(zhǔn)則[8,9]。
SCE-UA算法中大部分參數(shù)取值可以采用默認(rèn)值:
(4)
式中:n為參數(shù)個數(shù);m為每個復(fù)合形的頂點個數(shù);q為每個子復(fù)合形頂點個數(shù);s為種群大??;α、β為父輩產(chǎn)生的子輩個數(shù)與代數(shù);p為復(fù)合形個數(shù),需根據(jù)具體情況確定。
該算法流程見圖3。
圖3 SCE-UA算法流程Fig.3 SCE-UA algorithm flowchart
本文選取臺風(fēng)多發(fā)的寧波市為研究區(qū),區(qū)內(nèi)汛期期間(4-10月),受臺風(fēng)和熱帶風(fēng)暴登陸侵襲,造成全流域或部分流域的強(qiáng)降雨,降水量較大。
寧波天氣雷達(dá)基站位于慈溪市達(dá)蓬山,站點高程454.4 m?;緸镃INRAD/SA系列雷達(dá),空間分辨率為1 500 m,時間分辨率為5~6 min。本文在進(jìn)行去噪、地物回波等預(yù)處理后,插值成1 km柵格以供雷達(dá)回波圖像外推使用。
本文中選用于實驗的雷達(dá)回波圖像為2015年21號臺風(fēng)“杜鵑”期間2015-09-30,11∶00-12∶00,1 h內(nèi)12張雷達(dá)回波圖像中大于25 dBZ的區(qū)域。雷達(dá)回波區(qū)域為寧波雷達(dá)站附近350 km×380 km范圍區(qū)域,研究區(qū)選取寧波市附近250 km×230 km范圍區(qū)域(見圖4)。
從寧波市水文站收集所需的寧波全市63個國家雨量站點5 min降雨數(shù)據(jù),以供對預(yù)測雨強(qiáng)的精度分析使用。
圖4 雷達(dá)回波區(qū)域及研究區(qū)示意圖Fig.4 The sketch map of Radar echo area and study area
本文采用的外推方法是一種線性外推算法,不考慮雷達(dá)回波在移動過程中的非線性變化、垂直運動、強(qiáng)度演變等。因此,若外推時間太長,則此方法外推就沒有太大意義。在強(qiáng)對流天氣條件下,60 min內(nèi)雷達(dá)回波位置和形狀的外推預(yù)報具有一定的物理意義和可靠性,其結(jié)果具有較好的指示意義。本文考慮60 min的雷達(dá)外推時間。同時,雷達(dá)回波強(qiáng)度大于25 dBZ時,具有良好聚攏效果,有明顯的外廓邊界,更有利于雷達(dá)回波外推。因此,本文采用前1 h內(nèi)12張雷達(dá)回波圖像中大于25 dBZ的區(qū)域,分析并外推下1 h的雷達(dá)回波形狀及位置,利用Python語言實現(xiàn)圖像批處理及相關(guān)計算。
雷達(dá)外推的具體過程如下:首先裁剪出寧波雷達(dá)基站附近350 km×380 km的雷達(dá)回波區(qū)域作為雷達(dá)外推的搜索區(qū)域及寧波市附近250 km×230 km范圍的研究區(qū);以圖像每一個柵格大小為移動單元依次平移,用250 km×230 km裁剪窗口依次裁剪搜索區(qū)域內(nèi)雷達(dá)回波圖像,每張雷達(dá)回波圖像搜索區(qū)域內(nèi)共有100×130張裁剪出的雷達(dá)回波圖像以待計算(稱為計算窗口圖像);計算裁剪出的所有計算窗口雷達(dá)回波圖像與下一張圖像中研究區(qū)雷達(dá)回波圖像的相關(guān)系數(shù)(見表1、表2);保存與下一張研究區(qū)雷達(dá)回波圖像相關(guān)系數(shù)最大的計算窗口雷達(dá)回波圖像(見圖5、圖6),記錄該計算窗口網(wǎng)格中心到研究區(qū)網(wǎng)格中心的X方向與Y方向2個矢量(北為Y的正方向,東為X的正方向,反之則為負(fù));計算并記錄研究區(qū)雷達(dá)回波圖像與相關(guān)系數(shù)最大的計算窗口雷達(dá)回波圖像的回波強(qiáng)度變化量網(wǎng)格E;得到1 h內(nèi)12張雷達(dá)回波圖像中前11張圖像的X方向與Y方向2個移動矢量序列[(Δx1,Δy1),(Δx2,Δy2),…,(Δx11,Δy11)](見表3);疊加1 h內(nèi)前11張回波圖像的回波強(qiáng)度變化量網(wǎng)格∑E,作為下1 h預(yù)報雷達(dá)回波圖像形狀的依據(jù)。
圖5 與11∶06研究區(qū)雷達(dá)回波(彩色圖像)相關(guān)系數(shù) 最大的11∶00雷達(dá)回波(黑色輪廓)Fig.5 The radar echo of 11∶00 having largest correlation coefficient with the echo of 11∶06 in study area
圖6 與11∶12研究區(qū)雷達(dá)回波(彩色圖像)相關(guān)系數(shù) 最大的11∶06雷達(dá)回波(黑色輪廓)Fig.6 The radar echo of 11∶06 having largest correlation coefficient with the echo of 11∶12 in study area
將時間換算成時間戳形式,利用最小二乘算法建立X、Y分量累積量與時間的關(guān)系式(即以11∶00雷達(dá)回波圖像為準(zhǔn),各時間點的X、Y方向偏移量Δx、Δy):
Δx=273.5t-1 247.3R2=0.999 6
(5)
Δy=628.43t-286.21R2=0.993 5
(6)
將13∶00轉(zhuǎn)換為時間戳形式,代入以上2個關(guān)系式中,得到13∶00雷達(dá)回波相對11∶00時的X、Y偏移量,該偏移量減去12∶00雷達(dá)回波相對于11∶00的X、Y偏移量,可得到13∶00雷達(dá)回波相對于12∶00的X、Y偏移量(114,26);在12∶00雷達(dá)回波區(qū)域?qū)ふ?114,26)偏移量的相反方向窗口,即為13∶00研究區(qū)窗口的雷達(dá)回波位置;在該雷達(dá)回波窗口基礎(chǔ)上,疊加11∶00-12∶00的回波強(qiáng)度變化量網(wǎng)格E,即研究區(qū)內(nèi)13∶00的雷達(dá)外推預(yù)報圖像。
表1 11∶06研究區(qū)雷達(dá)回波與11∶00各計算窗口 的相關(guān)系數(shù)(由高到低)Tab.1 Correlation coefficient of the radar echo of 11∶00 in study area and each calculation window of 11∶06
表2 11∶12研究區(qū)雷達(dá)回波與11∶06各計算窗口 的相關(guān)系數(shù)(由高到低)Tab.2 Correlation coefficient of the radar echo of 11∶06 in study area and each calculation window of 11∶12
表3 2015-09-30 11∶00-12∶00 前11張雷達(dá) 回波圖像移動矢量序列(Δx,Δy)Tab.3 The first 11 radar echo image motion vector sequences of 2015-09-30 11∶00 to 12∶00
由圖7可以看出,考慮了雷達(dá)回波的移動過程及內(nèi)部量變的雷達(dá)外推預(yù)報圖像與回波實況接近,且該方法運算簡便、迅速,具有較高的實用價值和可操作性。
圖7 2015-09-30 13∶00的60 min雷達(dá)外推 預(yù)報圖像范圍(黑色輪廓)及實況(彩色圖像)Fig.7 60 min Radar extrapolation forecast image range (black outline) and actual image (color image) of 13∶00 on Sept 30,2015
雨量站每5 min一次實測降雨資料,按雷達(dá)探測時間整理成雨強(qiáng)數(shù)據(jù)。根據(jù)雷達(dá)反射率因子Z和降雨強(qiáng)度R的關(guān)系式Z=ARb及雷達(dá)反射因子Z與雷達(dá)回波強(qiáng)度N的關(guān)系式N=10 lgZ可知:
(7)
式中:A、b為待確定參數(shù)。
由這2個參數(shù)來確定雷達(dá)回波強(qiáng)度N和降雨強(qiáng)度R(單位:mm/h)的關(guān)系。本文對A,b2個參數(shù)用以下2種方法進(jìn)行估算。
(1)Lemon[10]提出考慮2種對流類型給定參數(shù):①大陸強(qiáng)對流型,A=300,b=1.4;②熱帶型,A=230,b=1.25。本文采用大陸強(qiáng)對流型參數(shù),利用雷達(dá)回波預(yù)報圖像,計算出預(yù)報降雨強(qiáng)度,用2015-09-30 13∶00降雨預(yù)報數(shù)據(jù)中有效雨量站點實測數(shù)據(jù)及預(yù)報數(shù)據(jù)求得相關(guān)系數(shù)為0.56,均方根誤差δ為9.35。
(2)SCE-UA算法可以有效、快速地搜索到參數(shù)全局的最優(yōu)解。為提高收斂速度,避免陷入局部最優(yōu)解,需給優(yōu)化參數(shù)設(shè)置上下限,建立其尋優(yōu)區(qū)間。式(7)中2個參數(shù)取值范圍為:A∈[200,350],b∈[0.7,1.8],作為參數(shù)尋優(yōu)的上下邊界。
參數(shù)優(yōu)化時樣本點需按目標(biāo)函數(shù)值進(jìn)行排序。本文構(gòu)造目標(biāo)函數(shù),為降雨實測值與參數(shù)優(yōu)化后的關(guān)系式計算所得估計值的殘差平方和最?。?/p>
(8)
式中:Ri為降雨強(qiáng)度的實測值,mm/h;n為樣本數(shù)量,i=1,2,…,n。
將樣本數(shù)據(jù)代入式(7)中,在滿足目標(biāo)函數(shù)的情況下,采用SCE-UA算法優(yōu)化參數(shù)A,b。算法中復(fù)合形個數(shù)p是唯一需要確定的參數(shù),本文p取值2,則m取值5,q取值3,α為1、β為5。算法終止條件主要從3方面控制:①迭代次數(shù)10次,容許值為0.1%;②最大循環(huán)次數(shù)10 000 次;③參數(shù)收斂值為0.001。
利用SCE-UA算法對雨強(qiáng)與回波強(qiáng)度關(guān)系式進(jìn)行參數(shù)估計,對60 min雷達(dá)外推預(yù)報成果進(jìn)行臨近降雨預(yù)報;利用雨量站實測數(shù)據(jù)進(jìn)行誤差分析計算,以均方根誤差作為衡量指標(biāo)。δ的計算公式如下:
(9)
本文采用11∶00-12∶00實測降雨資料及雷達(dá)回波強(qiáng)度數(shù)據(jù),得到參數(shù)A=318,b=1.12。用13∶00降雨預(yù)報數(shù)據(jù)中有效雨量站點實測數(shù)據(jù)及預(yù)報數(shù)據(jù)(見圖8)求得相關(guān)系數(shù)達(dá)到0.85,均方根誤差δ=2.48(見表4)。
事實上,A,b2參數(shù)對降雨預(yù)報精度的參數(shù)敏感性很高,尤其b參數(shù)的選取,對預(yù)報降雨精度影響很大。因此,對不同場降雨事件的A,b參數(shù)的率定是很有必要的。
給定研究區(qū)的雷達(dá)回波最優(yōu)空間相關(guān)方法對雷達(dá)回波位置和形狀進(jìn)行30~60 min的外推預(yù)報,雷達(dá)外推預(yù)報圖像與回波實況接近,具有一定的可行性和可靠性。在降雨預(yù)報方面,采用SCE-UA算法估算參數(shù)的預(yù)報精度明顯比采用經(jīng)驗參數(shù)的預(yù)報精度高,均方根誤差δ=2.48(見表4),不同場降雨事件應(yīng)率定相應(yīng)的參數(shù)。
□
圖8 2015-09-30 13∶00 實測雨強(qiáng)與雷達(dá)外推預(yù)測雨強(qiáng)Fig.8 Measured and radar extrapolation forecast rain intensity at 13∶00 on Sept 30,2015
方法參數(shù)A參數(shù)b均方根誤差δ經(jīng)驗參數(shù)法3001.409.35SCE-UA算法3181.122.48