於敏佳，劉菡，王倩，胡瀟杰

（舟山市氣象局，浙江舟山316000）

1 引言

舟山地處中緯度地帶，是“絲綢之路經(jīng)濟(jì)帶”和“21世紀(jì)海上絲綢之路”的交匯點，區(qū)位優(yōu)勢明顯。習(xí)近平總書記2015年視察舟山時指出：“舟山港口優(yōu)勢、區(qū)位優(yōu)勢、資源優(yōu)勢獨特，其開發(fā)開放不僅具有區(qū)域性的戰(zhàn)略意義，而且具有國家層面的戰(zhàn)略意義”。舟山港航經(jīng)濟(jì)發(fā)展蓬勃，寧波-舟山港連續(xù)10 a貨物吞吐量全球第一，寶鋼馬跡山、岙山興中和新奧舟山LNG等港口也承接著大量的礦石中轉(zhuǎn)和油氣儲運等業(yè)務(wù)。然而，舟山災(zāi)害性天氣四季頻發(fā)，年平均8級以上大風(fēng)日數(shù)142.7 d，其中8～9級大風(fēng)日數(shù)129.2 d，10～11級大風(fēng)日數(shù)12.7 d，12級以上大風(fēng)日數(shù)0.8 d[1]。同時舟山群島地形復(fù)雜，舟山沿海各海域受島嶼和山系遮蔽影響，使得各海域風(fēng)力具有明顯的局地性特征[2]。針對舟山各港口風(fēng)力局地性特點，業(yè)務(wù)上已開展了各港口風(fēng)力專業(yè)專項預(yù)報服務(wù)。經(jīng)過多年的EC數(shù)值模式預(yù)報與站點實況風(fēng)力對比，已得出各港口風(fēng)力預(yù)報實況偏差的定性經(jīng)驗判斷，但并未獲得精確定量的訂正數(shù)值。因此，應(yīng)用EC細(xì)網(wǎng)格數(shù)值預(yù)報產(chǎn)品，做好舟山各港口風(fēng)力預(yù)報釋用，為舟山港航經(jīng)濟(jì)發(fā)展提供氣象保障，具有十分重要的意義。

關(guān)于數(shù)值預(yù)報產(chǎn)品的釋用，目前我國已有較多的科技人員進(jìn)行了深入研究。唐躍等[3]采用雙線性內(nèi)插法和滑動平均誤差法，進(jìn)行舟山港口大風(fēng)的本地化釋用，并建立了各類指標(biāo)和預(yù)報模型。張娟等[4]將BP（Back Propagation）人工神經(jīng)網(wǎng)絡(luò)引入到風(fēng)暴潮數(shù)值預(yù)報的解釋應(yīng)用中，對風(fēng)暴潮的數(shù)值預(yù)報結(jié)果進(jìn)行訂正；陳德花等[5]應(yīng)用BP人工神經(jīng)網(wǎng)絡(luò)方法對MM5中尺度數(shù)值模式進(jìn)行了福建沿海大風(fēng)的釋用；黨英娜等[6]應(yīng)用極差分析訂正的方法對T639和MM5數(shù)值模式產(chǎn)品進(jìn)行釋用；陳豫英等[7]利用MM5模式的站點輸出產(chǎn)品，采用多元線性和逐步回歸的模式輸出統(tǒng)計方法（Model Output Statistics，MOS）進(jìn)行風(fēng)的精細(xì)化釋用；金龍等[8]利用模塊化模糊神經(jīng)網(wǎng)絡(luò)方法進(jìn)行數(shù)值預(yù)報產(chǎn)品釋用；劉鴻升等[9]設(shè)計了最大值-閥值綜合集成法，利用歐洲中期天氣預(yù)報中心（European Centre for Mediumrange Weather Forecasts，ECMWF）中期數(shù)值預(yù)報，進(jìn)行了偏北大風(fēng)的數(shù)值預(yù)報釋用；李江萍等[10]初步探討了統(tǒng)計降尺度法在數(shù)值預(yù)報產(chǎn)品釋用中的應(yīng)用；錢燕珍等[11]將支持向量機(jī)（Support Vector Machine，SVM）回歸方法應(yīng)用于近海和登陸熱帶氣旋的強(qiáng)度預(yù)報，對全球預(yù)報系統(tǒng)（Global Forecast System，GFS）數(shù)值預(yù)報產(chǎn)品進(jìn)行了釋用。以上方法都提高了預(yù)報準(zhǔn)確率，取得了一定的釋用效果，但對于地形的風(fēng)力影響研究較少，且隨著模式預(yù)報水平的不斷改進(jìn)，釋用模型也應(yīng)隨之更新。本文在研究釋用方法的基礎(chǔ)上，進(jìn)一步探究釋用模型的滾動訂正方法。

2 資料

模式資料采用EC細(xì)網(wǎng)格2017年1月1日—2019年6月30日每日20時（北京時，下同）再分析資料，其中2017年1月1日—12月31日為釋用模式建模樣本；2018年1月1日—12月31日為獨立釋用模式檢驗樣本；2018年1月1日—2019年6月30日為模型回歸參數(shù)滾動訂正樣本。

觀測資料為舟山嵊泗國家基準(zhǔn)站逐時10 m極大風(fēng)向風(fēng)速，該站站點觀測資料經(jīng)人工控制，數(shù)據(jù)完整，且具有較高的準(zhǔn)確性和代表性。

采用雙線性插值法將EC細(xì)網(wǎng)格10 m風(fēng)再分析產(chǎn)品插值到觀測站點上，站點信息見表1。

表1站點信息

3 釋用方法

3.1 站點風(fēng)速高度訂正

由于風(fēng)速隨高度的升高而增大，各站點在對比EC細(xì)網(wǎng)格10 m風(fēng)資料時，首先需要進(jìn)行風(fēng)速的高度訂正。根據(jù)之前的研究[12]，指數(shù)律模式能較好地反映風(fēng)速隨高度的變化。

指數(shù)律：

式中：u和u1分別為距地面高度z和z1處的風(fēng)速；m為風(fēng)速隨高度的切變指數(shù)，其值的大小即表明了風(fēng)速垂直切變的強(qiáng)度；嵊泗基準(zhǔn)站為草地下墊面，m取0.14。

如圖1所示，經(jīng)高度訂正后，EC再分析插值風(fēng)速與實況風(fēng)速誤差明顯改善。

3.2 考慮地形的PP法建模

在大尺度運動中，當(dāng)不考慮摩擦力和地形時，空氣質(zhì)點的運動遵從地轉(zhuǎn)風(fēng)規(guī)則，這是我們進(jìn)行風(fēng)預(yù)報的基礎(chǔ)。但在實際中，會有地形作用且摩擦力不為零的情況。

摩擦力方面：近地面層摩擦力來自于粗糙的地面及近地面層空氣的粘性力，而空氣粘性力主要由湍流活動引起，近地面層湍流與低層溫度層有關(guān)，即與大氣穩(wěn)定度相關(guān)。由于摩擦力的存在，定量表示地面風(fēng)有一定困難，難以精確估計摩擦力的大小。從實測風(fēng)的讀數(shù)可見，風(fēng)向和風(fēng)速均存在著明顯的脈動現(xiàn)象，只能以自記曲線的10 min平均作為平均風(fēng)，而以某一時段中最大平均風(fēng)作為該時段的平均風(fēng)[7]。挑選因子時應(yīng)從大氣穩(wěn)定度入手，尋求表示估計偏差風(fēng)的統(tǒng)計關(guān)系式?；谝陨侠碚摵腿斯そ?jīng)驗，挑選50個與大氣穩(wěn)定度有關(guān)的因子。這些因子包括4個層次（地面、850 hPa、700 hPa和500 hPa）的再分析場，即溫度、風(fēng)、比濕、相對濕度、溫度平流、水汽通量散度、渦度、垂直速度、散度、位渦和位勢高度等熱力動力因子。

圖1 2017年1月1日—12月31日每日20時EC插值點風(fēng)速、站點小時極大風(fēng)和高度訂正后站點小時極大風(fēng)對比

由于不同的物理量有不同量綱，需要先對各預(yù)報因子作標(biāo)準(zhǔn)化處理，本文采用正規(guī)化方法處理。應(yīng)用多元線性回歸方法對2017年整年觀測和再分析標(biāo)準(zhǔn)化資料建立回歸方程。經(jīng)計算，方程線性關(guān)系顯著性F=2.66，查F分布表F（50，314）的99%分位數(shù)為2，可見線性關(guān)系為顯著。計算出單個自變量顯著性t值（見表2），按0.01顯著性（t 0.99（314）=2.59）T檢驗標(biāo)準(zhǔn)對因子進(jìn)行排序篩選，保留了28個相關(guān)性較大的因子，建立PP法模型。

PP法回歸方程：

YS=XSB+ε （2）

式中：YS表示站點高度訂正標(biāo)準(zhǔn)化風(fēng)力；XS表示700 hPa位勢高度、850 hPa溫度、EC插值點風(fēng)力和1 000 hPa相對濕度等28個標(biāo)準(zhǔn)化預(yù)報因子；B表示待估參數(shù)；ε表示隨機(jī)誤差。

地形方面：風(fēng)速的大小除與水平氣壓梯度有直接關(guān)系外，與地形也有較大關(guān)系。舟山群島地形復(fù)雜，各觀測站點多有地形影響，主要包括山體阻擋使風(fēng)力減小、狹管效應(yīng)和爬坡使風(fēng)力增大等。地形不管以哪種形式對風(fēng)力產(chǎn)生影響，具體表現(xiàn)在站點風(fēng)力測量都在相同氣壓梯度下，不同風(fēng)向的站點風(fēng)速大小不同，因此，本文采用分風(fēng)向的方法，具體研究地形對站點風(fēng)速的影響。

圖2 嵊泗站各風(fēng)向多元回歸預(yù)報與高度訂正后風(fēng)速誤差箱線圖

按16方位的風(fēng)向進(jìn)行研究，分別計算嵊泗站在各個風(fēng)向條件下，PP法模型擬合與站點高度訂正的風(fēng)速誤差。如圖2所示，我們發(fā)現(xiàn)擬合風(fēng)速誤差和風(fēng)向角度呈現(xiàn)較明顯的三角函數(shù)分布特征，風(fēng)速誤差和風(fēng)向角度兩者的相關(guān)系數(shù)高達(dá)R2=0.82。對各風(fēng)向風(fēng)速進(jìn)行擬合，均方根誤差（Root Mean Squared Error，RMSE）為0.63，擬合效果較好（見圖3）。

故可得風(fēng)速在不同風(fēng)向角度的誤差公式如下：

式中：E表示各風(fēng)向角度下的擬合風(fēng)速誤差；x為風(fēng)向角度。據(jù)此建立考慮地形的PP法模型。

考慮地形的pp法回歸方程：

式中：Y表示站點高度訂正風(fēng)力；Y?表示多元回歸預(yù)報風(fēng)力；E表示各風(fēng)向角度下的擬合風(fēng)速偏差；e表示隨機(jī)誤差；XS表示700 hPa位勢高度、850 hPa溫度、EC插值點風(fēng)力和1 000 hPa相對濕度等28個標(biāo)準(zhǔn)化預(yù)報因子；B表示待估參數(shù)；S表示Y的標(biāo)準(zhǔn)差；-Y表示Y的平均值。

圖3 嵊泗站各風(fēng)向風(fēng)速誤差擬合圖

4 效果檢驗

4.1 訂正誤差評價方法

本文采用統(tǒng)計指標(biāo)參數(shù)平均絕對誤差（Mean Absolute Error，MAE）和RMSE評價模型訂正前后的風(fēng)速誤差。具體公式如下：

4.2 考慮地形的PP法模型檢驗

本文通過高度訂正、PP法多元線性回歸擬合和考慮地形的風(fēng)力誤差訂正，對2017年嵊泗基準(zhǔn)站高度訂正風(fēng)速和EC再分析資料進(jìn)行了預(yù)報建模。每經(jīng)一次訂正，誤差都有改進(jìn)，其中高度訂正后誤差改進(jìn)最明顯，最終各項評價指標(biāo)的訂正結(jié)果都超過50%，具體見表2。

利用2018年全年每日20時EC再分析資料，應(yīng)用所建模型進(jìn)行預(yù)報，對比嵊泗站點高度訂正后風(fēng)速。如圖4所示，模型預(yù)報誤差較原始EC插值點預(yù)報誤差有較大改進(jìn)，計算得出MAE為2.57，RMSE為3.27，可見模型具有較好的訂正能力。

5 卡爾曼濾波回歸系數(shù)訂正與檢驗

由于EC模式的預(yù)報性能和站點的氣候特征都在發(fā)生變化，因此模型中的回歸系數(shù)也應(yīng)隨時間發(fā)生變化。本文采用卡爾曼濾波法對模型回歸系數(shù)進(jìn)行滾動訂正。卡爾曼濾波是一種數(shù)學(xué)統(tǒng)計方法，該方法采用量測方程和狀態(tài)方程組成的線性隨機(jī)系統(tǒng)狀態(tài)空間模型來描述濾波器，并且利用狀態(tài)方程的遞推算法對濾波器的狀態(tài)變量做出最佳估計，從而求得濾掉噪聲的有用信號的最佳估計[13]。卡爾曼濾波一般用來處理濾波、預(yù)測和平滑3方面的問題。

天氣預(yù)報領(lǐng)域中的卡爾曼濾波的基本公式如下：

量測方程：

狀態(tài)方程：

表2嵊泗站點2017年1月1日—12月31日20時高度訂正前后、多元回歸擬合后和風(fēng)向修正后風(fēng)速誤差對比

圖4 2018年1月1日—12月31日每日20時原始誤差與模型預(yù)報誤差對比

圖5 2018年1月1日—2019年6月30日每日20時卡爾曼濾波滾動訂正前后風(fēng)速絕對誤差對比

6 結(jié)論

本文根據(jù)EC細(xì)網(wǎng)格資料和嵊泗站點實況資料，建立了考慮地形的PP法模型，用卡爾曼濾波法滾動更新回歸參數(shù)，以MAE和RMSE作為誤差訂正評價指標(biāo)，得出以下結(jié)論：

（1）對風(fēng)速的高度訂正，能有效降低EC細(xì)網(wǎng)格模式的預(yù)報誤差，MAE和RMSE由高度訂正前的4.47 m/s和5.51 m/s降至高度訂正后的2.70 m/s和3.53 m/s，誤差減幅分別為40%和36%。

（2）用PP法建立的站點高度訂正風(fēng)速和EC模式10 m風(fēng)場、700 hPa位勢高度和850 hPa溫度等預(yù)報因子間的多元線性回歸方程，方程線性關(guān)系顯著性F=2.66，查F分布表，F(xiàn)（50,314）的99%分位數(shù)為2，具有顯著的線性關(guān)系。

（3）分析各風(fēng)向的擬合風(fēng)速誤差，發(fā)現(xiàn)風(fēng)速誤差與風(fēng)向角度存在明顯的正余弦函數(shù)關(guān)系，擬合后相關(guān)系數(shù)R2=0.82，RMSE=0.63，擬合效果較好。

（4）經(jīng)高度訂正、多元線性回歸和分風(fēng)向誤差風(fēng)速訂正后，MAE和RMSE由高度訂正前的4.47 m/s和5.51 m/s降至模型訂正后的2.09 m/s和2.71 m/s，誤差降幅達(dá)到53%和51%，訂正效果明顯。

（5）利用2018年全年資料進(jìn)行模型訂正釋用能力檢驗。經(jīng)計算，模型預(yù)報風(fēng)速與站點高度訂正風(fēng)速的MAE為2.57 m/s，RMSE為3.27 m/s，模型具有較強(qiáng)的訂正釋用能力。

（6）利用卡爾曼濾波法，實時吸收最新EC細(xì)網(wǎng)格資料和站點實況風(fēng)速資料，滾動訂正所建模型的回歸系數(shù)，計算得出MAE和RMSE由訂正前的2.91 m/s和3.59 m/s滾動降至訂正后的2.66 m/s和3.34 m/s，可見該方法具有滾動更新能力，且能進(jìn)一步提高模型訂正釋用能力。