林逸凡， 董 勝

(中國(guó)海洋大學(xué)工程學(xué)院，山東 青島 266100)

風(fēng)作為一種自然現(xiàn)象，其本身蘊(yùn)含著巨大的能量，作為清潔可再生能源，風(fēng)能的開發(fā)利用日益受到關(guān)注和投入[1-2]。沿海地區(qū)受陸風(fēng)和海風(fēng)的共同影響，風(fēng)能資源相對(duì)內(nèi)陸地區(qū)更加豐富[3]。觀測(cè)、分析日常、極端天氣條件下的風(fēng)況在風(fēng)能、航運(yùn)、旅游業(yè)、漁業(yè)、城市空氣污染控制規(guī)劃等眾多行業(yè)皆有直接的需求[4]。

風(fēng)玫瑰圖用于描述風(fēng)速和風(fēng)向的分布情況，是統(tǒng)計(jì)分析風(fēng)氣象要素最常用的資料。而實(shí)際工程項(xiàng)目對(duì)于數(shù)據(jù)挖掘分析日益充分、可靠性日益嚴(yán)格，在此要求下，采用概率聯(lián)合分布模型分析風(fēng)速風(fēng)向分布更為準(zhǔn)確，建立的聯(lián)合概率模型可應(yīng)用與和風(fēng)氣象要素密切相關(guān)的污染擴(kuò)散[5]、風(fēng)能資源評(píng)價(jià)[6-7]，風(fēng)電場(chǎng)布置[8-9]等領(lǐng)域，相比風(fēng)玫瑰圖用途更為廣泛。

芝罘島位于山東省煙臺(tái)市西北約9 000 m海上，是我國(guó)最大的陸連島，受季風(fēng)環(huán)流及其他天氣因素影響明顯。本文基于角度-線性聯(lián)合分布模型，提出利用風(fēng)玫瑰圖建立風(fēng)速風(fēng)向聯(lián)合分布模型的方法，對(duì)芝罘島各季節(jié)風(fēng)玫瑰圖建立聯(lián)合概率分布模型，驗(yàn)證其可靠性。同時(shí)，基于條件概率理論，給出該聯(lián)合分布模型隨機(jī)數(shù)生成的相關(guān)算法，利用各季節(jié)風(fēng)向風(fēng)速聯(lián)合分布模型采用蒙特卡洛模擬對(duì)各方位可能出現(xiàn)的極值風(fēng)速進(jìn)行了預(yù)測(cè)分析。

1 數(shù)據(jù)和方法

1.1 數(shù)據(jù)來源及組成

芝罘島海洋環(huán)境氣象監(jiān)測(cè)站位置為37°36′N，121°26′E。收集得到該監(jiān)測(cè)站根據(jù)1981—1986年風(fēng)速風(fēng)向觀測(cè)數(shù)據(jù)繪制的春、夏、秋、冬四季節(jié)風(fēng)玫瑰圖(如圖1)及其對(duì)應(yīng)的風(fēng)速風(fēng)向頻率統(tǒng)計(jì)表。統(tǒng)計(jì)表記錄了16方位0～4級(jí)(0～7 m/s)、5級(jí)(8～10 m/s)、6級(jí)(11～13 m/s)、7級(jí)(14～17 m/s)以及8級(jí)以上(≥18 m/s)風(fēng)力的頻率，無風(fēng)頻率另記為C。另收集得到該測(cè)站統(tǒng)計(jì)的1981—1992年各月份、方位極大風(fēng)速數(shù)據(jù)。

1.2 研究方法及模型

1.2.1 風(fēng)速邊緣分布模型 在海洋工程及風(fēng)能領(lǐng)域，威布爾分布作為最常用的風(fēng)速分布模型，適用性較強(qiáng)?？紤]靜風(fēng)頻率C不為零，常采用由Diracδ函數(shù)和Weibull分布組成的混合模型描述風(fēng)速分布[10]，其密度函數(shù)如下：

fV(v)=F0δ(v)+(1-F0)fW(v)。

(1)

式中：v≥0；F0=C為靜風(fēng)頻率；δ(x)為Dirac函數(shù)，當(dāng)x=0時(shí)δ(0)=1，當(dāng)x≠0時(shí)δ(0)=0；fW(v)為威布爾分布密度函數(shù)：

(2)

式中，α,β>0分別為形狀和尺度參數(shù)。

當(dāng)v>0時(shí)，混合分布的累積頻率函數(shù)為

FV(v)=F0+(1-F0)FW(v)，

(3)

對(duì)式(3)進(jìn)行變換可得

FW(v)=(FV(v)-F0)/(1-F0)。

(4)

圖1 芝罘島四季風(fēng)玫瑰圖Fig.1 Seasonal wind roses of Zhifudao

令H表示可觀測(cè)風(fēng)(v>0),根據(jù)式(4)，威布爾分布FW(v)可作為其分布函數(shù)HV(v)。設(shè)各級(jí)風(fēng)力的風(fēng)速上限集合為[Vi]，i=1,…,5，其中V1=7，V2=10，…，V5設(shè)為23，由統(tǒng)計(jì)表可計(jì)算得到H的風(fēng)速累積頻率Pv，采用非線性最小二乘法(NLS)以累計(jì)頻率均方根誤差最小為目標(biāo)擬合威布爾分布。

令Xi=ln(Vi),Yi=ln[-ln(1-Pvi)]，α，β兩參數(shù)的理論最小二乘估計(jì)值可根據(jù)下式(5)計(jì)算。由于V5為大致估計(jì)值，故式(5)計(jì)算的分布參數(shù)只能作為參數(shù)非線性最小二乘擬合的初始迭代值。

(5.a)

(5.b)

(5.c)

均方根誤差指標(biāo)顯示各季度風(fēng)速邊緣分布擬合均良好。圖2為各季度風(fēng)速分布擬合圖像。

1.2.2 風(fēng)向邊緣分布模型 von Mise分布為圓周上的正態(tài)分布，其混合分布fvMs(θ)可描述有多個(gè)主風(fēng)向地區(qū)的風(fēng)向分布，密度函數(shù)為

(6)

式中：θ∈[0,2π)為角度變量；N為von Mises分布的混合組數(shù)，Carta等[11]研究表明6組von Mise分布即可提供足夠的擬合精度；wi>0為von Mise分布的權(quán)比例參數(shù)；0≤μi≤2π為位置參數(shù)；Кi≥0反映數(shù)據(jù)關(guān)于μi的離散程度，為集中度參數(shù)。I0(Кi)為零階修正第一類貝塞爾(Bessel)方程，即

圖2 四季風(fēng)速邊緣分布擬合圖Fig.2 Fitting plots of the wind speed distribution for every season

(7)

令正北偏西11.5°方向?yàn)榻嵌攘泓c(diǎn)，風(fēng)向值按順時(shí)針計(jì)，由統(tǒng)計(jì)表計(jì)算H的各向風(fēng)頻pd，可觀測(cè)風(fēng)von Mises分布hΘ(θ)的參數(shù)采用邊了解，邊近似的兩步擬合估計(jì)，具體步驟如圖3所示流程圖。圖4為擬合得到的各季度可觀測(cè)風(fēng)(H)風(fēng)向分布擬合圖。

1.2.3 二維風(fēng)速風(fēng)向聯(lián)合分布模型 Johnson和Wehrly[12]定義的角度-線性聯(lián)合分布模型用于可測(cè)量風(fēng)H，其密度函數(shù)為

hv,Θ=2πs(ψ)hV(v)hΘ(θ)。

(8)

式中：hv(v)和hΘ(θ)分別為H的風(fēng)速及風(fēng)向分布密度函數(shù)；ψ為風(fēng)速風(fēng)向相關(guān)系數(shù)。作為角度變量，通常采用兩組von Mises混合分布近似其密度函數(shù)s(ψ)[11]。

(9)

對(duì)于總體風(fēng)，當(dāng)v>0時(shí)，聯(lián)合分布的密度函數(shù)如下：

fv,Θ=2π(1-F0)g(ψ)hV(v)hΘ(θ)。

(10)

假設(shè)有一組完整風(fēng)速風(fēng)向數(shù)據(jù)集{(v,θ)}，篩選出風(fēng)向介于任意方位(θj,θj+1)之間，風(fēng)速介于某級(jí)風(fēng)力區(qū)間(vi,vi)的子數(shù)據(jù)集Cij。根據(jù)式(9)，該組數(shù)據(jù)的相關(guān)系數(shù)ψ將在兩個(gè)特征值ψ1到ψ2之間變化。而ψ1和ψ2的值又由HV(v1)，HV(v2)，HΘ(θ1)，HΘ(θ2)共同影響。對(duì)于任意方位某級(jí)風(fēng)力的可觀測(cè)風(fēng)頻pij(i=1,2,…,16;j=1,2,…,5)，由于HV(v)和HΘ(θ)的變化范圍都不大，其ψ值的波動(dòng)也不會(huì)很大，假設(shè)存在一個(gè)特征相關(guān)系數(shù)值ψij，則ψij和其對(duì)應(yīng)的風(fēng)頻pij將共同對(duì)相關(guān)系數(shù)的分布產(chǎn)生影響。基于該假設(shè)，相關(guān)系數(shù)分布參數(shù)的估計(jì)采用以下方法：

圖3 風(fēng)速von Mises混合分布擬合流程圖Fig.3 Process map of the method to fit the von Mises distribution of wind direction

對(duì)于每個(gè)風(fēng)向風(fēng)速區(qū)間Cij(vi

圖4 四季風(fēng)向邊緣分布擬合圖Fig.4 Fitting plots of the marginal distribution of wind direction for foun seasons

利用擬合得到的風(fēng)速威布爾分布HV，風(fēng)向von Mises混合分布HΘ，風(fēng)速風(fēng)向相關(guān)系數(shù)ψ的von Mises混合分布s(ψ)，結(jié)合已知的靜風(fēng)頻率F0，最終得到各季度風(fēng)速風(fēng)向聯(lián)合分布fV,Θ。

決定系數(shù)R2常用來反映模型與原始數(shù)據(jù)的相關(guān)程度，R2∈[-1,1]，越接近1表示兩者相關(guān)越強(qiáng)，當(dāng)?shù)扔?時(shí)表示兩者完全相關(guān)。計(jì)算公式如下：

(11)

表1 芝罘島各季節(jié)聯(lián)合分布模型決定系數(shù)R2Table 1 Coefficients of determination R2 for Zhibudao in four seasons

(①Cumulative frequency of theoretical joint distribution；②Cumulative frequency of original joint distribution；③Throretical frequency of wind rose diagram；④Original frequency of wind rose diagram.)

圖5 各季度季風(fēng)向風(fēng)速聯(lián)合分布擬合圖
Fig.5 Seasonal Fitting results of the joint probability model

2 芝罘島季節(jié)極端風(fēng)況蒙特卡洛模擬

2.1 隨機(jī)數(shù)生成算法

采用條件概率方法根據(jù)聯(lián)合分布模型生成風(fēng)向風(fēng)速二維偽隨機(jī)數(shù)[v,θ]，算法流程圖如圖6。其中Fcondition(x)為當(dāng)可觀測(cè)風(fēng)H風(fēng)速累積頻率HV=t時(shí)風(fēng)向累積頻率HΘ小于x的概率。

首先根據(jù)式(8)～(9)，推導(dǎo)得

Hv,Θ=[G2(A)-G2(B)+G2(C)]/2π。

(12)

式中：A=2(HV；B=2((HV-HΘ)；C=-2(HΘ，G2(·)為s(ψ)貝塞爾展開式的二重積分形式，按下式計(jì)算：

G2(x)=

(13)

圖6 風(fēng)速風(fēng)向二維偽隨機(jī)數(shù)生成算法流程圖

式中，[wi，μi，Кi]為s(ψ)的von Mise混合分布參數(shù)。

令HΘ為定值等于t，HV=x，對(duì)式(12)求偏導(dǎo)數(shù)得到Fcondition(x)，即

(14)

式中，G1(·)按下式計(jì)算：

G1(x)=

(15)

圖7為擬合得到的各季節(jié)風(fēng)向風(fēng)速聯(lián)合概率模型，按上述算法生成二維偽隨機(jī)數(shù)散點(diǎn)圖，其中灰度隨數(shù)據(jù)點(diǎn)聯(lián)合概率密度變大而加深。

2.2 年極值風(fēng)速蒙特卡洛模擬

采用蒙特卡洛法模擬100年各季度的風(fēng)速風(fēng)向數(shù)據(jù)，每組隨機(jī)風(fēng)速風(fēng)向數(shù)據(jù)代表3 h平均值，篩選每年各季度中16個(gè)方位的極值風(fēng)速，組成各季度、方位的極值風(fēng)速序列，采用箱型圖(包含：極小、25%、50%、75%、極大及異常值)對(duì)模擬極值風(fēng)速數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析，并與實(shí)測(cè)的芝罘島地區(qū)1981—1992年各季節(jié)、方位極值風(fēng)速數(shù)據(jù)進(jìn)行比較，結(jié)果如圖8。

3月份NW方位1981—1992年極值風(fēng)速為40 m/s。實(shí)測(cè)數(shù)據(jù)點(diǎn)基本均位于各方位100年模擬數(shù)據(jù)上下分位之間，同時(shí)相當(dāng)一部分位于上下四分之一分位范圍，模擬得到的各方位年極值風(fēng)速波動(dòng)范圍與實(shí)測(cè)資料吻合較好。將各方位100年極值風(fēng)速模擬數(shù)據(jù)排序，根據(jù)其經(jīng)驗(yàn)累積頻率推測(cè)各季節(jié)、方位不同重現(xiàn)期下芝罘島各季節(jié)方位極值風(fēng)速值，如圖9所示，結(jié)果如表2。

結(jié)合圖7和表2分析，芝罘島地區(qū)主要受三股季風(fēng)影響，各季節(jié)風(fēng)況變化明顯。其中春季主要受NNW-N和SSE季風(fēng)控制，極值風(fēng)速方向?yàn)镹-NNW，10、20和50年風(fēng)速分別為20.4、22.4和24.4 m/s；夏季NNW-N季風(fēng)減弱，主要受SSE季風(fēng)影響，但50年一遇極值風(fēng)速最大值仍出現(xiàn)在NNW，為23.4 m/s；進(jìn)入秋季，W季風(fēng)逐漸占據(jù)主導(dǎo)，但風(fēng)速主要在8 m/s以下，極值風(fēng)速仍為N-NNW，10、20和50年風(fēng)速為21.8、24.0和25.6 m/s；進(jìn)入冬季，雖然主風(fēng)向?yàn)閃向，但極值風(fēng)速為NNW，10、20和50年風(fēng)速為22.8、23.5和25.0 m/s。

圖7 芝罘島風(fēng)速風(fēng)向二維偽隨機(jī)數(shù)散點(diǎn)圖Fig.7 Scatterplots of the bivariate pseudo-random numbers of wind speed and direction Graphs at Zhifudao in four seasons

3 結(jié)語

根據(jù)風(fēng)玫瑰圖數(shù)據(jù)對(duì)煙臺(tái)芝罘島地區(qū)四季風(fēng)向風(fēng)速聯(lián)合分布建立模型，采用決定系數(shù)將風(fēng)玫瑰圖原始與理論數(shù)據(jù)進(jìn)行比較，結(jié)果顯示：聯(lián)合累積頻率和風(fēng)頻的決定系數(shù)大于0.9，表明根據(jù)風(fēng)玫瑰圖建立的聯(lián)合分布模型與原始風(fēng)玫瑰圖相關(guān)性極強(qiáng)。同時(shí)，利用聯(lián)合分布統(tǒng)計(jì)模型，結(jié)合蒙特卡洛法，模擬了該地區(qū)100年各季節(jié)風(fēng)速數(shù)據(jù)，從中提取各季方位年極值風(fēng)速，與實(shí)測(cè)數(shù)據(jù)比較表明：模擬的季節(jié)方位年極值風(fēng)速數(shù)據(jù)可較好反映該地區(qū)極值風(fēng)速波動(dòng)范圍，利用其對(duì)不同重現(xiàn)期風(fēng)速方位極值進(jìn)行預(yù)測(cè)估計(jì)。對(duì)于缺乏風(fēng)速資料地區(qū)，可采用該方法分析通常及極端狀況下的風(fēng)況是一條有效的途徑。

圖8 芝罘島四季方位極值風(fēng)速模擬序列箱型圖Fig.8 Boxplots of simulative seasonal directional extreme wind speed

圖9 重現(xiàn)值-序號(hào)對(duì)應(yīng)圖Fig.9 Corresponding diagram of return value and sorted number

春Sring夏Summer秋Autumn冬WinterV10V20V50V10V20V50V10V20V50V10V20V50N20.4021.2522.5416.2316.9917.5421.0921.7825.6222.0422.8724.24NNE17.4918.5920.5517.0719.8820.4020.8822.7823.8720.9221.9023.70NE17.8919.9222.0516.3217.8719.1919.8421.1322.3418.9120.4622.48ENE18.3719.7222.1415.6816.8418.4818.8620.7722.5517.1519.4421.87E16.5917.5118.3816.0517.6519.1517.1017.8818.3017.0519.1020.70ESE18.1519.3820.6219.1019.9421.3618.7220.1022.7718.7419.8820.93SE20.1221.6422.9517.6318.7219.3020.2022.1622.2720.7822.2823.15

續(xù)表2

春Sring夏Summer秋Autumn冬WinterV10V20V50V10V20V50V10V20V50V10V20V50SSE20.2920.8923.5918.6319.0620.7921.2821.7322.2619.1121.2822.55S19.0020.6021.6517.7718.9520.0418.9819.9920.8218.4420.9821.13SSW19.3320.5422.2916.0517.9120.1817.3318.4420.5419.0020.5222.04SW19.0920.6021.8116.6818.0221.1018.5821.3225.1119.5020.7523.52WSW19.5020.4321.1916.4316.8218.2321.7123.0025.1421.2622.4123.11W19.4120.5021.6316.2117.5017.7520.4122.3023.2421.2222.3524.27WNW18.1519.0320.1416.3916.8017.5420.8623.1626.9420.9023.0023.95NW19.7620.0620.5917.8618.8119.6921.5123.5824.8120.4222.5623.28NNW19.8322.4224.4217.4619.0423.3621.8423.9824.7322.7923.5324.96季節(jié)極值N20.40NNW22.42NNW24.42SSE18.63SSE19.88NNW23.36NNW21.84NNW23.98N26.94NNW22.79NNW23.53NNW24.96

注：V10、V20和V50分別為10、20、50 a一遇極值風(fēng)速預(yù)測(cè)值，單位為m/s。

Note:V10，V20andV50are 10, 20 and 50 a retura wind speeds, respectively, unit is m/s.