雷達回波外推方法在臨近降雨預報中的應用

張衛(wèi)國，范仲麗，鐘 偉，江雨田，孫飛飛，陳 娟

(1.寧波市水利水電規(guī)劃設計研究院，浙江 寧波 315192；2.河海大學水文水資源學院，南京 210098)

數(shù)值預報對強對流天氣的臨近預報技術還不夠成熟，但強對流天氣引發(fā)的氣象災害對民眾的生命財產(chǎn)安全造成極大的危害。因此，對強對流天氣的預警和臨近預報意義重大。而雷達估測降雨在城市監(jiān)測預警方面的應用將是未來水利防洪領域的前沿[1]。

天氣雷達具有高時空分辨率的優(yōu)勢，可以快速、大范圍捕捉強對流天氣情況，也因此越來越多地被應用到強對流天氣短時臨近降雨預報中。雷達回波最優(yōu)空間相關方法最早被用于Austin和Bellon[2]研究中，利用平均移動矢量對降雨進行了預報；在此基礎上，Rinehart和Gravey[3]發(fā)展了TREC(Tracking Radar Echo By Correction)方法，將雷達回波分為多個區(qū)域，得到不同的移動矢量；Li[4]等在TREC基礎上發(fā)展了COTREC方法，并進行了降水預報。而多數(shù)利用TREC方法的研究中，只作雷達回波位置和大概形狀的外推，對回波內(nèi)部的量變不作修正。

雨強的臨近估計主要根據(jù)雷達回波反射率因子Z和雨強R的關系式進行。在多數(shù)研究中，該關系式中參數(shù)都采用固定參數(shù)或經(jīng)驗參數(shù)。但實際上，不同地區(qū)或者不同降雨過程中，參數(shù)值存在較大差別，使得雷達降雨估測與實際降雨偏差較大。因此，利用實測雨量數(shù)據(jù)對雨強關系式中參數(shù)的實時校準，以減少更多的估算誤差是很有必要的。

本文以臺風多發(fā)的寧波市為例，利用給定研究區(qū)的雷達回波最優(yōu)空間相關方法，得到前一時次每6 min的雷達回波移動矢量，對該矢量的2個分量分別利用最小二乘回歸得到其與時間的關系，由此可得下一時次的移動矢量；疊加上一時次雷達回波內(nèi)部變化量后，得到雷達回波預報圖像，實現(xiàn)對雷達回波位置、形狀的外推及內(nèi)部回波量的修正。對Z～R關系式進行變形，采用SCE-UA算法對雷達回波強度N與降雨強度R關系中參數(shù)進行估算與實時校準，嘗試預測下一時次的降水分布情況。

1 理論方法

本文采用給定研究區(qū)的雷達回波最優(yōu)空間相關方法對強對流天氣下雷達回波的位置和形狀進行外推，并把上一時次的雷達回波變化網(wǎng)格疊加到下一時次的雷達回波預報中，得到最終的雷達預報圖像。根據(jù)雷達反射率因子與雨強的關系式的變形公式對雨強進行估算，利用SCE-UA算法對雨強關系式中參數(shù)進行實時校準，得到時效性較好的A，b參數(shù)，從而得到相對準確的下一時次預報降雨分布。

1.1 雷達回波最優(yōu)空間相關方法

從雷達回波圖像上可以看出，雷達探測回波存在明顯的、較為明確的移動軌跡(見圖1)。根據(jù)已發(fā)生的短時間內(nèi)的雷達回波移動軌跡，推測下一時次雷達回波的位置及形狀，即是雷達回波最優(yōu)空間相關方法的核心思想。

雷達回波最優(yōu)空間相關方法，是在假設大氣中水汽的移動受環(huán)境風場的引導，而環(huán)境風場在短時間內(nèi)不存在嚴重突變的前提下成立的。將已探測到的雷達回波的移動矢量作為當前環(huán)境風的移動矢量，來預測未來短時間內(nèi)雷達回波的移動矢量。該方法簡便，運算量小，運行時間短，在有大面積層狀云時跟蹤效果較好[5]。

圖1 雷達回波移動軌跡示意圖(2018年4月4日)Fig.1 The schematic of radar echo movement track (April 4, 2018)

給定研究區(qū)的雷達回波最優(yōu)空間相關方法具體步驟是，將整個回波區(qū)域按照研究區(qū)的網(wǎng)格大小劃分成若干個相同大小的矩形網(wǎng)格，研究區(qū)雷達回波矩形網(wǎng)格與上一時刻的若干矩形網(wǎng)格之間計算相關系數(shù)，尋得的最大相關系數(shù)的網(wǎng)格到研究區(qū)網(wǎng)格對應的矢量即為該天氣系統(tǒng)下相鄰時刻雷達回波的位移矢量。

相關系數(shù)R的計算公式如下：

(1)

尋找與研究區(qū)網(wǎng)格最大相關系數(shù)的網(wǎng)格時，為縮減計算矩形網(wǎng)格的數(shù)量，可根據(jù)雨帶移動的最大期望速度與分析時間間隔得到最大網(wǎng)格搜索半徑。為減小網(wǎng)格距離計算量，可采用最大網(wǎng)格搜索半徑為邊長的正方形作為搜索區(qū)域。按網(wǎng)格大小為一個移動單位進行移動，求得在最大搜索半徑范圍內(nèi)的所有網(wǎng)格(計算區(qū))與下一相鄰時次研究區(qū)網(wǎng)格的相關系數(shù)。其中相關系數(shù)最大的計算區(qū)網(wǎng)格中心到研究區(qū)網(wǎng)格中心的矢量即為t1到t2時刻的位移矢量(見圖2)，并保存該計算區(qū)網(wǎng)格與研究區(qū)網(wǎng)格的變化量網(wǎng)格E。重復尋找當前時刻之前一定時間段內(nèi)相鄰時次的雷達回波圖像的位移矢量(Δx1,Δy1)，(Δx2,Δy2)，…,(Δxn,Δyn)，記錄一系列變化量網(wǎng)格(E1,E2,…,En)。假設環(huán)境風場方向及風速在一定時間段內(nèi)穩(wěn)定，可以分別建立X方向與Y方向上的位移分量Δx、Δy與絕對時間t之間的線性關系式：

Δx=a1t+b1

(2)

Δy=a2t+b2

(3)

根據(jù)至少2組數(shù)據(jù)[t1,(Δx1,Δy1)]，[t2,(Δx2,Δy2)]即可確定a1、b1、a2、b24個參數(shù)，理論上可用來計算未來任意時間與其相鄰前一時次間的雷達回波位移矢量。在實際運用時，需要考慮方法的時效性。

圖2 給定研究區(qū)的雷達回波最優(yōu)空間相關方法Fig.2 Radar echo optimal spatial correlation method for a given study area

在確定雷達回波外推位置與形狀后，疊加上一時段的變化量網(wǎng)格，可以得到最終的外推雷達回波圖像。

1.2 SCE-UA算法

復合形交叉進化算法(SCE-UA)是一種全局優(yōu)化算法，集成了隨機搜索算法、單純形法、聚類分析法等方法的優(yōu)點，是一種有效解決非線性約束最優(yōu)化問題的混合算法[6,7]。該算法引入了種群的概念。復合形點在可行域內(nèi)隨機生成、競爭演化，在多個吸引域內(nèi)獲得全局收斂點，能有效地表達參數(shù)敏感性與參數(shù)之間的相關性。SCE-UA算法基本思路是首先在參數(shù)可行域內(nèi)引入隨機分布的點群，將其分為多個復合形，每個復合形中有2n+1個點(其中n是需要優(yōu)化的參數(shù)個數(shù))；每個復合形根據(jù)下降單純形算法進行進化；定時地將整個點群混合進而形成新的復合形，也就將之前的復合形所包含的信息引入新的復合形中；進化與混合不斷進行直至滿足收斂準則[8,9]。

SCE-UA算法中大部分參數(shù)取值可以采用默認值：

(4)

式中：n為參數(shù)個數(shù)；m為每個復合形的頂點個數(shù)；q為每個子復合形頂點個數(shù)；s為種群大?。沪?、β為父輩產(chǎn)生的子輩個數(shù)與代數(shù)；p為復合形個數(shù)，需根據(jù)具體情況確定。

該算法流程見圖3。

圖3 SCE-UA算法流程Fig.3 SCE-UA algorithm flowchart

2 應用實例

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

2.1 數(shù)據(jù)準備

寧波天氣雷達基站位于慈溪市達蓬山，站點高程454.4 m。基站為CINRAD/SA系列雷達，空間分辨率為1 500 m，時間分辨率為5～6 min。本文在進行去噪、地物回波等預處理后，插值成1 km柵格以供雷達回波圖像外推使用。

本文中選用于實驗的雷達回波圖像為2015年21號臺風“杜鵑”期間2015-09-30，11∶00-12∶00，1 h內(nèi)12張雷達回波圖像中大于25 dBZ的區(qū)域。雷達回波區(qū)域為寧波雷達站附近350 km×380 km范圍區(qū)域，研究區(qū)選取寧波市附近250 km×230 km范圍區(qū)域(見圖4)。

從寧波市水文站收集所需的寧波全市63個國家雨量站點5 min降雨數(shù)據(jù)，以供對預測雨強的精度分析使用。

圖4 雷達回波區(qū)域及研究區(qū)示意圖Fig.4 The sketch map of Radar echo area and study area

2.2 考慮雷達回波形變與量變的外推預報

本文采用的外推方法是一種線性外推算法，不考慮雷達回波在移動過程中的非線性變化、垂直運動、強度演變等。因此，若外推時間太長，則此方法外推就沒有太大意義。在強對流天氣條件下，60 min內(nèi)雷達回波位置和形狀的外推預報具有一定的物理意義和可靠性，其結果具有較好的指示意義。本文考慮60 min的雷達外推時間。同時，雷達回波強度大于25 dBZ時，具有良好聚攏效果，有明顯的外廓邊界，更有利于雷達回波外推。因此，本文采用前1 h內(nèi)12張雷達回波圖像中大于25 dBZ的區(qū)域，分析并外推下1 h的雷達回波形狀及位置，利用Python語言實現(xiàn)圖像批處理及相關計算。

雷達外推的具體過程如下：首先裁剪出寧波雷達基站附近350 km×380 km的雷達回波區(qū)域作為雷達外推的搜索區(qū)域及寧波市附近250 km×230 km范圍的研究區(qū)；以圖像每一個柵格大小為移動單元依次平移，用250 km×230 km裁剪窗口依次裁剪搜索區(qū)域內(nèi)雷達回波圖像，每張雷達回波圖像搜索區(qū)域內(nèi)共有100×130張裁剪出的雷達回波圖像以待計算(稱為計算窗口圖像)；計算裁剪出的所有計算窗口雷達回波圖像與下一張圖像中研究區(qū)雷達回波圖像的相關系數(shù)(見表1、表2)；保存與下一張研究區(qū)雷達回波圖像相關系數(shù)最大的計算窗口雷達回波圖像(見圖5、圖6)，記錄該計算窗口網(wǎng)格中心到研究區(qū)網(wǎng)格中心的X方向與Y方向2個矢量(北為Y的正方向，東為X的正方向，反之則為負)；計算并記錄研究區(qū)雷達回波圖像與相關系數(shù)最大的計算窗口雷達回波圖像的回波強度變化量網(wǎng)格E；得到1 h內(nèi)12張雷達回波圖像中前11張圖像的X方向與Y方向2個移動矢量序列[(Δx1,Δy1)，(Δx2,Δy2)，…，(Δx11,Δy11)](見表3)；疊加1 h內(nèi)前11張回波圖像的回波強度變化量網(wǎng)格∑E，作為下1 h預報雷達回波圖像形狀的依據(jù)。

圖5 與11∶06研究區(qū)雷達回波(彩色圖像)相關系數(shù) 最大的11∶00雷達回波(黑色輪廓)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ū)雷達回波(彩色圖像)相關系數(shù) 最大的11∶06雷達回波(黑色輪廓)Fig.6 The radar echo of 11∶06 having largest correlation coefficient with the echo of 11∶12 in study area

將時間換算成時間戳形式，利用最小二乘算法建立X、Y分量累積量與時間的關系式(即以11∶00雷達回波圖像為準，各時間點的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轉換為時間戳形式，代入以上2個關系式中，得到13∶00雷達回波相對11∶00時的X、Y偏移量，該偏移量減去12∶00雷達回波相對于11∶00的X、Y偏移量，可得到13∶00雷達回波相對于12∶00的X、Y偏移量(114，26)；在12∶00雷達回波區(qū)域尋找(114，26)偏移量的相反方向窗口，即為13∶00研究區(qū)窗口的雷達回波位置；在該雷達回波窗口基礎上，疊加11∶00-12∶00的回波強度變化量網(wǎng)格E，即研究區(qū)內(nèi)13∶00的雷達外推預報圖像。

表1 11∶06研究區(qū)雷達回波與11∶00各計算窗口 的相關系數(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ū)雷達回波與11∶06各計算窗口 的相關系數(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張雷達 回波圖像移動矢量序列(Δx,Δy)Tab.3 The first 11 radar echo image motion vector sequences of 2015-09-30 11∶00 to 12∶00

由圖7可以看出，考慮了雷達回波的移動過程及內(nèi)部量變的雷達外推預報圖像與回波實況接近，且該方法運算簡便、迅速，具有較高的實用價值和可操作性。

圖7 2015-09-30 13∶00的60 min雷達外推 預報圖像范圍(黑色輪廓)及實況(彩色圖像)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 臨近降雨預測

雨量站每5 min一次實測降雨資料，按雷達探測時間整理成雨強數(shù)據(jù)。根據(jù)雷達反射率因子Z和降雨強度R的關系式Z=ARb及雷達反射因子Z與雷達回波強度N的關系式N=10 lgZ可知：

(7)

式中：A、b為待確定參數(shù)。

由這2個參數(shù)來確定雷達回波強度N和降雨強度R(單位：mm/h)的關系。本文對A，b2個參數(shù)用以下2種方法進行估算。

(1)Lemon[10]提出考慮2種對流類型給定參數(shù)：①大陸強對流型，A=300，b=1.4；②熱帶型，A=230，b=1.25。本文采用大陸強對流型參數(shù)，利用雷達回波預報圖像，計算出預報降雨強度，用2015-09-30 13∶00降雨預報數(shù)據(jù)中有效雨量站點實測數(shù)據(jù)及預報數(shù)據(jù)求得相關系數(shù)為0.56，均方根誤差δ為9.35。

(2)SCE-UA算法可以有效、快速地搜索到參數(shù)全局的最優(yōu)解。為提高收斂速度，避免陷入局部最優(yōu)解，需給優(yōu)化參數(shù)設置上下限，建立其尋優(yōu)區(qū)間。式(7)中2個參數(shù)取值范圍為：A∈[200，350]，b∈[0.7，1.8]，作為參數(shù)尋優(yōu)的上下邊界。

參數(shù)優(yōu)化時樣本點需按目標函數(shù)值進行排序。本文構造目標函數(shù)，為降雨實測值與參數(shù)優(yōu)化后的關系式計算所得估計值的殘差平方和最?。?/p>

(8)

式中：Ri為降雨強度的實測值，mm/h；n為樣本數(shù)量，i=1,2,…,n。

將樣本數(shù)據(jù)代入式(7)中，在滿足目標函數(shù)的情況下，采用SCE-UA算法優(yōu)化參數(shù)A，b。算法中復合形個數(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算法對雨強與回波強度關系式進行參數(shù)估計，對60 min雷達外推預報成果進行臨近降雨預報；利用雨量站實測數(shù)據(jù)進行誤差分析計算，以均方根誤差作為衡量指標。δ的計算公式如下：

(9)

本文采用11∶00-12∶00實測降雨資料及雷達回波強度數(shù)據(jù)，得到參數(shù)A=318，b=1.12。用13∶00降雨預報數(shù)據(jù)中有效雨量站點實測數(shù)據(jù)及預報數(shù)據(jù)(見圖8)求得相關系數(shù)達到0.85，均方根誤差δ=2.48(見表4)。

事實上，A，b2參數(shù)對降雨預報精度的參數(shù)敏感性很高，尤其b參數(shù)的選取，對預報降雨精度影響很大。因此，對不同場降雨事件的A，b參數(shù)的率定是很有必要的。

3 結 語

給定研究區(qū)的雷達回波最優(yōu)空間相關方法對雷達回波位置和形狀進行30～60 min的外推預報，雷達外推預報圖像與回波實況接近，具有一定的可行性和可靠性。在降雨預報方面，采用SCE-UA算法估算參數(shù)的預報精度明顯比采用經(jīng)驗參數(shù)的預報精度高，均方根誤差δ=2.48(見表4)，不同場降雨事件應率定相應的參數(shù)。

□

圖8 2015-09-30 13∶00 實測雨強與雷達外推預測雨強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