結(jié)合中尺度模式物理約束的雷達回波臨近外推預(yù)報方法研究*

孫泓川 吳海英 曾明劍 程叢蘭

1.江蘇省氣象臺，南京，210008

2.中國氣象局交通氣象重點開發(fā)實驗室，南京，210009

3.江蘇省氣象科學(xué)研究所，南京，210009

4.北京城市氣象研究院，北京，100089

1 引言

盡管數(shù)值天氣預(yù)報的最新進展使人們能夠以很高的分辨率預(yù)報、預(yù)測，但是對于頻繁更新的臨近預(yù)報應(yīng)用要求，計算成本過高，因此臨近預(yù)報應(yīng)用上，基于雷達回波的外推方法明顯勝過數(shù)值模式預(yù)報。方法主要有交叉相關(guān)法（陳明軒等，2007；黃旋旋等，2017）、光流法（曹春燕等，2015；Ayzel，et al，2019）、粒子濾波融合法（Chen，et al，2017）等，近年來對抗網(wǎng)絡(luò)（陳元昭等，2019）、循環(huán)神經(jīng)網(wǎng)絡(luò)（韓豐等，2019；郭瀚陽等，2019）等人工智能技術(shù)也在臨近外推預(yù)報上有所應(yīng)用。光流法通過從連續(xù)時次雷達回波圖像中計算光流場，用光流場代替雷達回波運動矢量場，并基于該運動矢量場對雷達回波進行外推。大量學(xué)者開展了光流法和交叉相關(guān)法的外推預(yù)報效果的定量對比評估。結(jié)果表明，對移動型局地生成及強度和形狀隨時間變化較快的雷暴及鋒面低槽降水，光流法優(yōu)于交叉相關(guān)法（吳劍坤等，2019）?？傮w而言，近20 年的時間里，光流算法一直都在全球最先進的臨近預(yù)報業(yè)務(wù)系統(tǒng)中起著非常重要的作用（Ayzel，et al，2019）。

目前的外推預(yù)報方法主要基于雷達反射率單個因子，基本上沒有考慮熱動力環(huán)境對回波強度變化的影響，因此這些方法大多對預(yù)報回波的強度變化趨勢預(yù)報能力欠佳（曹春燕等，2015）。這里需要指出：基于人工智能技術(shù)的外推為獲得目標(biāo)函數(shù)最優(yōu)在一定程度會調(diào)整回波強度（會形成諸如降低中心強度，擴大回波范圍等情況迎合目標(biāo)函數(shù)），但這也并非是嚴(yán)格意義上的預(yù)測回波強度變化。因此單外推方法的研究（陳元昭等，2019；韓豐等，2019）往往指出，臨近外推預(yù)報除了需要利用海量雷達數(shù)據(jù)給模型學(xué)習(xí)外，必須要結(jié)合高時空分辨率的動力、熱力信息提高模型的預(yù)測能力，例如利用數(shù)值模式的預(yù)報信息。其實，很多中外短時臨近研究（Wong，et al，2009；Atencia，et al，2010；Haiden，et al，2011；Liang，et al，2010；程叢蘭等，2013，2019）均明確了臨近外推預(yù)報和數(shù)值模式二者進行融合預(yù)報是目前提升預(yù)報能力的主要手段。但目前這些研究的融合方法主要是在一定預(yù)報時效后對外推和模式的結(jié)果進行權(quán)重融合，而從初始預(yù)報時效就開始融入了模式熱動力信息從而對外推預(yù)報進行約束的研究相對較少。

鑒于上述方法存在的不足，本研究嘗試使用智能外推方法來改進對流的強度變化預(yù)報。將人工智能技術(shù)與氣象的物理模型相結(jié)合，能夠從海量的氣象數(shù)據(jù)中提取并解讀有用的信息，發(fā)現(xiàn)回波強度的變化，實現(xiàn)更精準(zhǔn)的短時臨近外推預(yù)報。文中選取江蘇地區(qū)多部雷達組網(wǎng)數(shù)據(jù)的長時序歷史資料作為試驗對象，在光流法外推基礎(chǔ)上在初始預(yù)報時效即利用中尺度高分辨率模式的物理約束改進整個外推預(yù)報時效的質(zhì)量。具體方法分為3 步：（1）融合模式風(fēng)場修正光流風(fēng)場；（2）利用針對強對流的模式釋用產(chǎn)品和雷達歷史資料構(gòu)建訓(xùn)練集，使用投票回歸器集成多種深度學(xué)習(xí)算法構(gòu)建智能回波強度預(yù)測模型；（3）利用智能回波強度預(yù)測模型實時預(yù)測預(yù)報時效回波強度頻率分布，并使用降水頻率匹配訂正技術(shù)（李俊等，2014；Zhu，et al，2015；蘇翔等，2021）對外推預(yù)測的雷達回波強度進行訂正。最后基于這種方法檢驗了兩個強對流天氣過程的預(yù)報個例，也給出長時間序列的檢驗評分結(jié)果。

2 資 料

2.1 雷達資料

文中使用的雷達資料來自江蘇及其周邊主要的SA 雷達，具體分布見圖1，紅色曲線范圍內(nèi)是雷達的覆蓋區(qū)，其中藍色框線是研究范圍（29°—37°N，114°—122°E），研究區(qū)域基本都在雷達覆蓋區(qū)內(nèi)。在研究范圍內(nèi)，利用組網(wǎng)技術(shù)計算多部雷達的組網(wǎng)產(chǎn)品，利用組網(wǎng)產(chǎn)品中的組合反射率產(chǎn)品進行外推預(yù)報。本研究利用2018 年汛期雷達數(shù)據(jù)建模，研究時段為2019—2020 年主汛期（6—9 月）。

圖1 選取的試驗范圍（藍色方框為選取的試驗區(qū)域（29°—37°N，114°—122°E）；紅色圓點為雷達站點位置，數(shù)字為雷達站號，紅色曲線范圍為雷達站點探測覆蓋范圍）Fig.1 Distribution of radar stations（red points）and radar detection range（red line）（The blue rectangle denotes the selected area（29°—37°N，114°—122°E））

2.2 其他資料

被用于對外推結(jié)果進行物理約束試驗的資料包括江蘇本地精細化天氣分析預(yù)報系統(tǒng)（Li，et al，2016）（簡稱為PWAFS）的預(yù)報及釋用產(chǎn)品。PWAFS 由WRF 模式和WRFDA 資料同化系統(tǒng)構(gòu)成，模式采用雙層嵌套，水平分辨率分別為15 km和3 km，垂直分為45 層，每3 h 起報一次，每日運行8 次。文中主要使用：

（1）模式業(yè)務(wù)預(yù)報的500 hPa 至850 hPa 多層風(fēng)場；（2）基于模式產(chǎn)品診斷的分類強對流概率預(yù)報產(chǎn)品，時間分辨率為1 h；（3）基于模式產(chǎn)品診斷的強對流潛勢指數(shù)，時間分辨率為1 h。

3 研究方法

本研究基于光流法（Optical Flow method，OF）的雷達回波外推預(yù)報技術(shù)（Ayzel，et al，2019），該方法采用了金字塔分層技術(shù)（Liu，et al，2015）對不同尺度雷達回波進行追蹤從而減小由于回波移速較快造成的反演誤差，提高反演風(fēng)場的計算精度和運算效率；采用半拉格朗日外推方案改進外推，保持回波的旋轉(zhuǎn)性，提高回波預(yù)報效果。

3.1 光流風(fēng)場填充

光流法通過對雷達回波圖像追蹤來獲取風(fēng)場，在無回波區(qū)域基于圖像追蹤無法獲得風(fēng)場。因此，會導(dǎo)致雷達回波無法走出原始回波區(qū)域，造成回波在原始風(fēng)場邊緣出現(xiàn)堆積扭曲的現(xiàn)象。文中參考INCA系統(tǒng)（Haiden，et al，2011）引入對應(yīng)時效PWAFS 風(fēng)場來補充光流法的風(fēng)場。下面將介紹光流風(fēng)場填充的步驟。

首先采用500、700 和850 hPa 三層模式風(fēng)場的矢量平均，用調(diào)和余弦譜展開方法（Chen，et al，1992；鄧滌菲等，2011）進行無輻散處理（文中光流法為了防止回波外推時結(jié)構(gòu)發(fā)散對光流風(fēng)場做了無輻散處理，所以對需要融合的模式風(fēng)場也做同樣處理）。具體做法為：在引入流函數(shù)（ψ）和速度勢（χ）后，水平風(fēng)矢量可以分解為無輻散風(fēng)分量和無旋轉(zhuǎn)風(fēng)分量，相應(yīng)的流函數(shù)（ψ）和速度勢（χ）滿足以下兩個Possion方程（式（1）、（2）），其中 ?是垂直渦度，D是水平散度，在滿足有限區(qū)域給定的邊界條件（式（3）、（4））下求解方程。

模式風(fēng)場經(jīng)無輻散處理后，存在光流風(fēng)場的區(qū)域以光流風(fēng)場為基準(zhǔn)，而無光流風(fēng)場區(qū)域以模式風(fēng)場為基準(zhǔn)，二者過渡區(qū)域通過對二者加權(quán)的方式生成。具體權(quán)重計算方式如下：首先建立初始權(quán)重矩陣w0，其中光流風(fēng)場大于5 m/s 區(qū)域權(quán)重為1，其余為0。然后對該矩陣進行σ=1，半徑為5 像素的高斯模糊。根據(jù)高斯分布公式（式（5））計算高斯掩膜，其中x為像素點坐標(biāo)，u為中心經(jīng)緯度坐標(biāo)。

將建立好的11×11 大小的高斯掩膜G對初始權(quán)重矩陣進行卷積處理（式（6）），生成加權(quán)系數(shù)（w）的矩陣。最后每個格點將模式風(fēng)場和光流風(fēng)場進行加權(quán)融合（式（7）、（8））得到融合后的風(fēng)場。

圖2 是外推預(yù)報的原始光流風(fēng)場、處理后的模式風(fēng)場和動力約束后二者融合的風(fēng)場示例?？梢钥闯觯朔椒ń鉀Q了無回波區(qū)域風(fēng)場估計的問題。文中把填充了PWAFS 風(fēng)場的光流外推方法稱為動力 約 束 外 推（Dynamic Constraint Optical Flow method，DCOF）。

圖2 動力約束風(fēng)場融合效果（a.光流風(fēng)場（大于 5 m/s 區(qū)域用紫色實線圈出），b.模式無輻散風(fēng)場，c.融合風(fēng)場）Fig.2 Comparison between radar the optical flow model produced wind field（a）（the purple line denotes the area where the wind speed is greater than 5 m/s），PWAFS harmonic wind field（b）and the combined wind field（c）

3.2 回波強度智能調(diào)整

3.2.1 物理約束因子

曾明劍等（2015a）基于PWAFS 模式預(yù)報的對流參數(shù)，采用相對偏差模糊矩陣評價技術(shù)構(gòu)建了分類強對流天氣預(yù)報概率。同時，曾明劍等（2015b）基于模式預(yù)測和實況信息利用接近度概念開發(fā)了強對流潛勢預(yù)報指數(shù)，這些方法都已經(jīng)用于實際預(yù)報取得較好的業(yè)務(wù)效果。不同于一般的強對流熱動力指標(biāo)和模式對流指標(biāo)產(chǎn)品，這些產(chǎn)品不僅綜合了熱動力信息，還同時考慮到氣候情況和PWAFS模式本身的偏差。研究發(fā)現(xiàn)若逐格點考察其和回波強度的關(guān)系，之間并不存在特別好的相關(guān)關(guān)系，但如果考慮研究區(qū)域整體情況，回波強度和不同模式強對流產(chǎn)品則呈現(xiàn)出明顯的相關(guān)。圖3 是研究區(qū)域2019 年7 月回波大于不同閾值（10、20、30 和40 dBz）的格點數(shù)時間序列以及和這些產(chǎn)品（多種分類強對流概率和強對流潛勢預(yù)報指數(shù)）統(tǒng)計之后其中相關(guān)最優(yōu)的一種產(chǎn)品的歸一化時間序列，圖中最優(yōu)的產(chǎn)品分別是大于0.3、0.4 閾值的短時強降水概率格點數(shù)（圖3a、b）和大于0.1 閾值的強對流潛勢格點數(shù)（圖3c、d）?？梢钥闯霾煌考壍幕夭ㄇ闆r和這些強對流產(chǎn)品存在明顯的相關(guān)，其中最優(yōu)的相關(guān)系數(shù)分別達到0.40、0.41、0.54、0.61，且通過p=0.05 顯著性t檢驗，都達到了中等相關(guān)和強相關(guān)的程度。這說明可以利用這些相關(guān)較強的模式強對流預(yù)報產(chǎn)品信息通過機器學(xué)習(xí)技術(shù)來預(yù)測外推預(yù)報時效內(nèi)研究區(qū)域整體回波強度的分布情況，從而改善外推的回波強度預(yù)報。文中選擇了相關(guān)最好的短時強降水概率、強對流潛勢指數(shù)兩種產(chǎn)品，統(tǒng)計兩種產(chǎn)品大于不同閾值的格點比例作為預(yù)測模型的物理因子。

圖3 回波強度和其最優(yōu)相關(guān)的強對流產(chǎn)品歸一化時間序列（a.回波強度大于10 dBz，b.回波強度大于20 dBz，c.回波強度大于30 dBz，d.回波強度大于40 dBz）Fig.3 Normalized time series of precipitation echo and the most relevant convective products（a.precipitation echo above 10 dBz，b.precipitation echo above 20 dBz，c.precipitation echo above 30 dBz，d.precipitation echo above 40 dBz）

3.2.2 構(gòu)建訓(xùn)練集

訓(xùn)練集時間是2018 年主汛期（6—9 月），按照外推預(yù)報更新頻率每6 min 取一次樣本，在2018年6—9 月內(nèi)去除缺失數(shù)據(jù)時次后樣本共計20000 多次。首先針對短時強降水概率、強對流潛勢指數(shù)兩種產(chǎn)品，設(shè)定0.1—0.6 間隔0.1 為閾值（兩種產(chǎn)品值較少出現(xiàn)大于0.6 情況），分別統(tǒng)計其在研究區(qū)域大于不同閾值的格點比例作為建模的物理因子產(chǎn)品（每次2 h 的預(yù)報時效內(nèi)共有2 個時次模式物理因子產(chǎn)品對應(yīng)（因為產(chǎn)品為1 h 預(yù)報間隔），以其構(gòu)建預(yù)報因子。同時把每個起報時刻前1 h 研究區(qū)域內(nèi)雷達回波在不同間隔內(nèi)的格點比例也構(gòu)建為預(yù)報因子。然后，利用預(yù)報時效的雷達回波在不同強度間隔內(nèi)的格點比例作為預(yù)報量。文中對不同數(shù)據(jù)樣本進行預(yù)處理之后按照回波強度區(qū)間并分預(yù)報時效分別構(gòu)建訓(xùn)練集，詳見表1。訓(xùn)練集模型利用3.2.1 節(jié)選出來的物理約束因子結(jié)合起報時刻雷達回波和較前雷達回波強度變化預(yù)測共同預(yù)報未來預(yù)報時段的回波在不同強度區(qū)間（[5 dBz，10 dBz）、[10 dBz，15 dBz）、[15 dBz，20 dBz）、[20 dBz，25 dBz）、[25 dBz，30 dBz）、[30 dBz，35 dBz）、[35 dBz，∞））格點比例。

表1 訓(xùn)練集構(gòu)建方式Table 1 Construction of training set

3.2.3 回波強度智能預(yù)測模型

文中采用機器學(xué)習(xí)中的一種投票回歸器算法（Voting Regressor）（Pedregosa，et al，2011）將概念上不同的機器學(xué)習(xí)回歸器組合起來，以便平衡它們各自的弱點，并返回平均預(yù)測值。文中的投票回歸器集成了SVR（Support Vector Regression）、隨機森林（Breiman，2001）、Adaboost（Adaptive Boosting）、多元線性回歸4 種回歸算法。其中SVR 為支持向量回歸，實際就是支持向量機在做分類問題上的推廣；隨機森林算法是一種分類樹算法，通過組合多個分類器解決單一預(yù)測問題。Adaboost 是自適應(yīng)增強算法，根據(jù)回歸誤差率等參數(shù)來更新每個樣本的權(quán)重迭代計算形成強回歸器。具體流程見圖4。雖然訓(xùn)練過程耗時較長，但訓(xùn)練好的機器學(xué)習(xí)模型則可以直接部署在業(yè)務(wù)平臺上進行實時的預(yù)測產(chǎn)品推送，做到時間尺度為秒級的預(yù)測。

圖4 回波強度頻率智能預(yù)測模型流程Fig.4 Flow chart of intelligent echo intensity frequency prediction model

基于訓(xùn)練好的模型實時預(yù)測每次外推預(yù)報0—2 h 不同預(yù)報時效的回波在不同強度間隔內(nèi)的格點比例，即可得出回波強度頻率分布曲線。之后使用頻率匹配訂正方法對外推預(yù)報回波進行調(diào)整，具體做法是基于原始頻率曲線和預(yù)測頻率曲線，把原始預(yù)報曲線中的回波值調(diào)整為在預(yù)測頻率曲線中具有同樣頻率的點所對應(yīng)的回波值（李俊等，2014），從而實現(xiàn)高分辨率中尺度模式物理信息對外推預(yù)報的回波強度調(diào)整，即物理約束外推方法（Physical Constraint Optical Flow method，PCOF）。

4 檢驗評估

4.1 個例檢驗

2019 年7 月6 日華東地區(qū)在東北冷渦背景下出現(xiàn)罕見大范圍冰雹等強對流天氣，圖5a 給出了7 月6 日11 時（北京時，下同）的組合反射率因子實況，圖5b、c、d 分別為原始光流外推、動力約束外推、物理約束外推09 時起報的120 min臨近預(yù)報。

圖5 2019 年7 月6 日09 時雷達外推預(yù)報120 min 效果對比（a.實況，b.光流外推，c.動力約束外推，d.物理約束外推）Fig.5 Comparison between observations and forecasts on 6 July 2019（a.observations of composite reflectivity at 11：00 BT，b.120 min forecast at 09：00 BT using OF，c.120 min forecast at 09：00 BT using DCOF，d.120 min forecast at 09：00 BT using PCOF）

對比原始光流（圖5b）和動力約束外推方法（圖5c）可以看出：原始光流外推法的120 min 預(yù)報在回波的邊緣出現(xiàn)了明顯的扭曲堆積（該現(xiàn)象在60 min 預(yù)報已經(jīng)出現(xiàn)，120 min 更為明顯，圖略）與真實回波形態(tài)出現(xiàn)偏差；而加入動力約束外推以后，系統(tǒng)整體移動趨勢更合理，特別是在系統(tǒng)邊緣區(qū)域通過模式風(fēng)場的調(diào)整，預(yù)報與實況更接近，這種改進隨時效延長（如圖所示在120 min）預(yù)報效果更明顯。這次預(yù)報個例實況從起報時次開始有加強的趨勢（圖6），相較于其他兩種方法，物理約束在20—30 dBz 和超過30 dBz 兩個區(qū)間格點數(shù)明顯和觀測更為接近（原始光流法和動力約束方法并不改變回波強度，強回波格點數(shù)減少的原因是部分回波移出研究區(qū)域）?？梢钥闯霰M管動力約束外推可以有效提高雷達回波預(yù)測準(zhǔn)確率，但是使用動力約束無法實現(xiàn)回波強度調(diào)整的目標(biāo)，所以最終整體回波強度比實況弱。對比動力約束試驗方法（圖5c）和物理約束試驗方法（圖5d）可以看出：物理約束通過之前建立的回波強度頻率預(yù)測模型，預(yù)測的回波強度整體要比動力約束方法強，所以通過對動力約束外推方法預(yù)報進行回波強度頻率匹配訂正，加強了整體回波強度，在120 min這個預(yù)報時效的強度和觀測更為接近。

圖6 2019 年7 月6 日09 時雷達外推預(yù)報和實況在研究范圍內(nèi)20—30 dBz（a）、超過30 dBz（b）的格點數(shù)時間序列Fig.6 Time series of number of grid points with radar echoes within 20—30 dBz（a）and above 30 dBz（b）from radar extrapolation forecast and observations in the study area at 09：00 BT 6 July 2019

圖7a、b 給出了2019 年7 月6 日09 時起報的0—120 min 預(yù)報時效內(nèi)幾種外推方法的客觀檢驗結(jié)果（不同閾值的CSI 評分），分析可知：根據(jù)20 dBz的回波強度CSI 表現(xiàn)，前期動力約束外推和原始光流外推相差不大，40 min 預(yù)報時效之后逐漸體現(xiàn)出改進效果，改進隨著預(yù)報時效延長而增大，這對應(yīng)著原始光流外推在回波邊緣扭曲堆積現(xiàn)象的逐漸出現(xiàn)，而動力約束外推能明顯改進這種不真實的現(xiàn)象；30 dBz 的回波強度CSI 表現(xiàn)說明動力約束外推一直優(yōu)于原始光流外推，可以看出原始光流外推的扭曲堆積現(xiàn)象對強回波強度的CSI 影響更為明顯，動力約束的改進也更為明顯。本次預(yù)報個例中回波強度頻率預(yù)測模型預(yù)測的回波強度頻率分布整體高于動力約束外推結(jié)果，通過頻率匹配訂正，物理約束外推方法的回波強度整體有一定加強，更符合實況。對不同強度的CSI 評分，物理約束都明顯優(yōu)于原始光流和動力約束，尤其是在相對強回波（30 dBz）下物理約束外推的CSI 評分改進更為明顯。物理約束相比于其他方法的訂正效果隨預(yù)報時效的延長而增大。聯(lián)合BIAS 評分進行評判，整體上物理約束外推的BIAS 評分更接近于1，只是初始時效由于物理約束加強回波強度，一定程度增加了空報，但是在后面都明顯較之另外兩種方法更好的平衡了空報和漏報。

圖7e 給出同一個例的0—120 min 預(yù)報時效的均方根誤差（RMSE）結(jié)果，同樣可以看出本次預(yù)報個例中物理約束外推和動力約束外推也明顯優(yōu)于原始光流法，其中物理約束外推由于增加了整體回波強度，均方根誤差比動力約束外推略有增大，但基本上不存在大的差異。圖7f 給出同一個例的0—120 min 預(yù)報時效的相關(guān)系數(shù)檢驗結(jié)果（通過p=0.05 顯著性t檢驗），呈現(xiàn)物理約束外推優(yōu)于動力約束外推，動力約束外推優(yōu)于原始光流外推的情況。

圖7 2019 年7 月6 日09 時3 種外推預(yù)報方案（光流外推：黑線，動力約束外推：藍線，物理約束外推：紅線）20 dBz 的CSI評分（a）和BIAS 評分（b），30 dBz 的CSI 評分（c）和BIAS 評分（d），以及均方根誤差（e）和相關(guān)系數(shù)（f）Fig.7 Experimental results in terms of CSI（a.20 dBz，c.30 dBz），BIAS（b.20 dBz，d.30 dBz），RMSE（e）and CC（f）by OF（black line），DCOF（blue line）and PCOF（red line）at 09：00 BT 6 July 2019

2020 年5 月17—18 日受高空冷渦東移南下的影響，沿淮和淮北地區(qū)出現(xiàn)雷雨大風(fēng)等強對流天氣，局部地區(qū)出現(xiàn)冰雹。強對流天氣主要出現(xiàn)在17 日21 時—18 日03 時。圖8a 給 出 了5 月18 日03 時的組合反射率因子實況，圖8b、c、d 分別為原始光流法、動力約束外推、物理約束外推02 時起報的60 min 臨近預(yù)報。

對比原始光流外推（圖8b）和動力約束外推方法（圖8c）可以看出：原始光流外推方法的60 min預(yù)報在回波的邊緣（主要位置在山東省黃海海上區(qū)域）已經(jīng)出現(xiàn)了明顯的扭曲堆積，而加入動力約束外推以后，系統(tǒng)整體移動趨勢明顯更為合理，與實況更接近。但也可以看出約束外推無法實現(xiàn)實況中回波強度逐漸減弱的情況，在60 min 預(yù)報中回波強度整體比實況要強。對比動力約束試驗方法（圖8c）和物理約束試驗方法（圖8d）可以看出：本次個例選取的時段是強對流時段的相對后期，整體回波強度有減弱趨勢（圖9）。其中物理約束較之其他兩種方法更為接近觀測，通過對回波強度的頻率匹配訂正，使得回波強度有所減弱，可以看出在這個預(yù)報時效內(nèi)的回波強度和觀測更加接近。

圖8 2020 年5 月18 日02 時雷達外推預(yù)報60 min 效果對比（a.實況，b.光流外推，c.動力約束外推，d.物理約束外推）Fig.8 Comparison between observations and forecasts on 18 May 2020（a.observations of composite reflectivity at 03：00 BT，b.60 min forecast at 02：00 BT using OF，c.60 min forecast at 02：00 BT using DCOF，d.60 min forecast at 02：00 BT using PCOF）

圖9 2020 年5 月18 日02 時雷達外推預(yù)報和實況在研究范圍內(nèi)20—30 dBz（a）、超過30 dBz（b）的格點數(shù)時間序列Fig.9 Time series of number of grid points with radar echoes within 20—30 dBz（a）and above 30 dBz（b）from radar extrapolation forecast and observations in the study area at 02：00 BT 18 May 2020

圖10a、b 給出了2020 年5 月18 日02 時起報的0—120 min 預(yù)報時效內(nèi)幾種外推方法的客觀檢驗結(jié)果（不同閾值的CSI 評分），分析可知：20 dBz的回波強度CSI 顯示，前期動力約束外推和熱力約束外推對比原始光流法相差不大，但在20 min 預(yù)報時效之后有所改進，且改進隨著預(yù)報時效延長而增大，這對應(yīng)著光流外推風(fēng)場缺失的問題，而物理約束外推對回波強度的調(diào)整相對動力約束又有進一步改善；30 dBz 的回波強度CSI 表現(xiàn)說明動力約束和物理約束外推都優(yōu)于原始光流外推，而物理約束外推相對更好一些，可以看出在強回波強度的CSI表現(xiàn)上動力約束和物理約束外推改進影響更為明顯。結(jié)合BIAS 評分進行聯(lián)合評判，此個例整體上物理約束外推的BIAS 評分更接近于1。初始時效由于物理約束加強回波強度，一定程度增加了漏報，但在中期大部分時效都更好地平衡了空報和漏報，只是在較長的幾個預(yù)報時效這種優(yōu)勢有所下降。

圖10e 給出同一個例的0—120 min 預(yù)報時效的均方根誤差，同樣可以看出物理約束外推和動力約束外推也明顯優(yōu)于原始光流外推，在此個例中，由于整體回波強度調(diào)整減弱，物理約束外推在均方根誤差的優(yōu)勢上比CSI 略微明顯。圖10f 給出同一個例的0—120 min 預(yù)報時效的相關(guān)系數(shù)檢驗結(jié)果（通過p=0.05 顯著性t檢驗），呈現(xiàn)物理約束外推略優(yōu)于動力約束外推優(yōu)于原始光流外推的情況。

圖10 2020 年5 月18 日02 時3 種外推預(yù)報方案（光流外推：黑線，動力約束外推：藍線，物理約束外推：紅線）20 dBz 的CSI評分（a）和BIAS 評分（b），30 dBz 的CSI 評分（c）和BIAS 評分（d），以及圴方根誤差（e）和相關(guān)系數(shù)（f）Fig.10 Experimental results in terms of CSI（a.20 dBz，c.30 dBz），BIAS（b.20 dBz，d.30 dBz），RMSE（e）and CC（f）by OF（black line），DCOF（blue line）and PCOF（red line）at 02：00 BT 18 May 2020

4.2 統(tǒng)計檢驗

為了系統(tǒng)地檢驗幾種方法的效果，圖11 給出了2019、2020 年主汛期（6—9 月）整體的預(yù)報檢驗結(jié)果,包括20、30 dBz 的CSI 評分和BIAS 評分、均方根誤差、相關(guān)系數(shù)（通過p=0.05 顯著性t檢驗）。

CSI 的結(jié)果顯示，整體上，物理約束外推優(yōu)于動力約束外推，動力約束外推優(yōu)于原始光流外推。訂正效果隨預(yù)報時效有變化，在24 min 前，動力約束外推和原始光流外推差不多，物理約束外推調(diào)整過強度之后還略差一些。24 min 之后，物理約束外推明顯優(yōu)于動力約束外推和原始光流外推，改進效果隨著預(yù)報時效延長而增大。分析表明：在初始預(yù)報時效，原始光流法外推的堆積扭曲現(xiàn)象并不明顯，所以相應(yīng)的動力約束外推對其并沒有明顯的改進，當(dāng)預(yù)報時效增長，動力約束外推能更好改進原始光流外推存在的問題；同樣，預(yù)報初期，在短時間內(nèi)回

波強度變化較小，基于原始光流外推不對回波強度調(diào)整也可以反映真實的降水實況，而基于模式物理因子預(yù)測出的回波頻率也可能存在誤差，所以在前幾個預(yù)報時次，物理約束外推改進并不明顯。但隨著預(yù)報時效的延長，光流法外推不調(diào)整回波的強度則會和實況有很大的偏差，所以物理約束外推的改進效果隨著預(yù)報時效越來越明顯，尤其是在更強一些的回波量級CSI 改進更為明顯。同時結(jié)合BIAS評分進行聯(lián)合評判，整體上3 種外推方法都是漏報大于空報，其中物理約束外推的BIAS 評分相對更接近于1。

從均方根誤差的角度，物理約束外推和動力約束外推都優(yōu)于原始光流外推，但同時也注意到物理約束外推的均方根誤差要略高于動力約束外推。這表明：總體而言，基于中尺度模式信息對回波頻率進行訂正時，明顯改善了CSI 評分的同時均方根誤差也略有增大，這種情況在降水頻率匹配訂正時也偶有出現(xiàn)，這是因為頻率匹配訂正的目標(biāo)并不是均方根誤差而是分間隔的頻率分布。在相關(guān)系數(shù)檢驗中呈現(xiàn)物理約束外推優(yōu)于動力約束外推，動力約束外推優(yōu)于原始光流外推。

5 結(jié)論和討論

詳細介紹了基于中尺度模式物理約束的智能外推預(yù)報方法，通過個例分析和統(tǒng)計檢驗對比了原始光流外推、動力約束外推、物理約束外推3 種方案的0—120 min 預(yù)報性能，得到以下主要結(jié)論。

（1）利用基于中尺度模式物理約束的智能外推預(yù)報方法對2019、2020 年的兩個強對流個例進行檢驗。結(jié)果表明：動力約束外推通過改善光流法回波在邊緣的堆積扭曲改進了預(yù)報性能，物理約束外推通過基于模式信息預(yù)測的回波頻率分布調(diào)整回波強度實現(xiàn)回波的增強或減弱從而改善預(yù)報性能，隨著預(yù)報時效延長改善越來越明顯。

（2）2019、2020 年主汛期整體的預(yù)報檢驗結(jié)果表明：在初始預(yù)報時效，動力約束對原始光流外推結(jié)果并沒有明顯的改進，當(dāng)預(yù)報時效延長，動力約束能更好改進原始光流外推結(jié)果；同樣，在預(yù)報初期，物理約束外推甚至還不如其他方法，但隨著預(yù)報時效的延長，改進越來越明顯，聯(lián)合BIAS 綜合評判，物理約束外推整體上是最優(yōu)的方案。

雖然動力約束外推和物理約束外推方法通過修正光流風(fēng)場，回波強度調(diào)整改進了光流法降水回波預(yù)報能力，但是同時也應(yīng)注意到只能以預(yù)報中有回波區(qū)的回波增強或減弱從而在一定程度上反映降水系統(tǒng)生消的變化，對于新生回波的預(yù)報欠佳。以后研究將進一步利用諸如衛(wèi)星等多源觀測資料改善該方法中初生對流的預(yù)報問題。