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

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

1 引言

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

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

2 資料

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

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

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

表1站點(diǎn)信息

3 釋用方法

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

由于風(fēng)速隨高度的升高而增大，各站點(diǎn)在對(duì)比EC細(xì)網(wǎng)格10 m風(fēng)資料時(shí)，首先需要進(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)速與實(shí)況風(fēng)速誤差明顯改善。

3.2 考慮地形的PP法建模

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

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

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

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

PP法回歸方程：

YS=XSB+ε （2）

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

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

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

按16方位的風(fēng)向進(jìn)行研究，分別計(jì)算嵊泗站在各個(gè)風(fēng)向條件下，PP法模型擬合與站點(diǎn)高度訂正的風(fēng)速誤差。如圖2所示，我們發(fā)現(xiàn)擬合風(fēng)速誤差和風(fēng)向角度呈現(xiàn)較明顯的三角函數(shù)分布特征，風(fēng)速誤差和風(fēng)向角度兩者的相關(guān)系數(shù)高達(dá)R2=0.82。對(duì)各風(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表示站點(diǎn)高度訂正風(fēng)力；Y?表示多元回歸預(yù)報(bào)風(fēng)力；E表示各風(fēng)向角度下的擬合風(fēng)速偏差；e表示隨機(jī)誤差；XS表示700 hPa位勢(shì)高度、850 hPa溫度、EC插值點(diǎn)風(fēng)力和1 000 hPa相對(duì)濕度等28個(gè)標(biāo)準(zhǔn)化預(yù)報(bào)因子；B表示待估參數(shù)；S表示Y的標(biāo)準(zhǔn)差；-Y表示Y的平均值。

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

4 效果檢驗(yàn)

4.1 訂正誤差評(píng)價(jià)方法

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

4.2 考慮地形的PP法模型檢驗(yàn)

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

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

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

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

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

量測(cè)方程：

狀態(tài)方程：

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

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

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

6 結(jié)論

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

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

（2）用PP法建立的站點(diǎn)高度訂正風(fēng)速和EC模式10 m風(fēng)場(chǎng)、700 hPa位勢(shì)高度和850 hPa溫度等預(yù)報(bào)因子間的多元線性回歸方程，方程線性關(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)行模型訂正釋用能力檢驗(yàn)。經(jīng)計(jì)算，模型預(yù)報(bào)風(fēng)速與站點(diǎn)高度訂正風(fēng)速的MAE為2.57 m/s，RMSE為3.27 m/s，模型具有較強(qiáng)的訂正釋用能力。

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