張衛(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ù)測下一時次的降水分布情況。

1 理論方法

本文采用給定研究區(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ù)報降雨分布。

1.1 雷達(dá)回波最優(yōu)空間相關(guān)方法

從雷達(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á)回波圖像。

1.2 SCE-UA算法

復(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

2 應(yīng)用實例

本文選取臺風(fēng)多發(fā)的寧波市為研究區(qū)，區(qū)內(nèi)汛期期間(4-10月)，受臺風(fēng)和熱帶風(fēng)暴登陸侵襲，造成全流域或部分流域的強(qiáng)降雨，降水量較大。

2.1 數(shù)據(jù)準(zhǔn)備

寧波天氣雷達(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

2.2 考慮雷達(dá)回波形變與量變的外推預(yù)報

本文采用的外推方法是一種線性外推算法，不考慮雷達(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

2.3 臨近降雨預(yù)測

雨量站每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ù)的率定是很有必要的。

3 結(jié) 語

給定研究區(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