基于測風塔數據的地表粗糙度計算方法研究

項 雯，彭秀芳，施 晨，馮 浩

(中國能源建設集團江蘇省電力設計院有限公司，江蘇 南京 211102)

0 引言

風能資源評估作為建設風電場可行性的基礎，是發展風力發電的重要前提條件。采用數值模擬手段可以獲取開發利用高度處風能資源量的區域分布狀況，而地表粗糙度對近地面風速的數值模擬結果影響較大，因此，在風能資源評估中合理確定地表粗糙度值是一個非常重要的問題。

理論上，地表粗糙度z0是平均風速隨高度減小到零時的高度[1]。趙曉松[2]等根據牛頓迭代法計算零平面位移，進而得到地表粗糙度；Bruin[3]等根據質量守恒定律計算地表粗糙度。Martano[4]根據Monin Obukhov相似理論，利用單層超聲風溫資料，將計算零平面位移和空氣動力學粗糙度問題簡化為一個可由最小二乘法求單變量的過程，單獨求出z0。這些方法基本上需要利用多層風溫資料或者一層三維超聲數據，而在實際工程項目中，一般收集到的測風塔數據通常只有一層溫度數據，基本沒有三維超聲數據，這就給運用上述類似方法估算地表粗糙度帶來極大阻礙。

基于此，本文對平坦地形地表粗糙度計算方法展開研究，提出一種僅以測風塔實測數據為基礎計算地表粗糙度值的方法，并通過實例進行驗證。本文方法所需資料簡單易獲取，為今后設計人員估算地表粗糙度值、準確評估風能資源提供幫助。

1 粗糙度計算方法

1.1 風廓線法

在近地面層(即產生地點為貼地層至其上方50～100 m高度處[5])風速隨高度是否遵從對數規律分布主要由當時的大氣層結決定。大氣層結可分為不穩定、穩定和中性三種，在靜力中性條件下，風速隨高度呈現對數變化規律。

在描述近地面層結構時最常用的就是Monin Obukhov相似理論，即近地面層的無量綱風速[6]：

式中：z為距離地面高度，m；u為z高度上平均水平風速，m/s；κ為von Karman常數，通常取值0.4；u*為摩擦風速，m/s；L為Monin-Obukhov長度，m。

在中性層結時，L→∞，φ(Z/L)→1，對式(1)進行高度積分，得到各層高度風速：

式中：z0為地表粗糙度。zi為測風塔第i個風速儀安裝高度(以下簡稱“第i層高度”)，m；為 第i層高度在中性大氣層結下的風速平均值(以下簡稱“第i層高度中性風速平均值”)，m/s。考慮到近地面層高度通常不超過100m，因此本文規定zi≤100 m。

將測風塔第一層高度中性風速平均值作為參照風速，可以由式(2)得：

相鄰高度風速之間兩兩相減，得：

由式(3)可以看出若粗糙度z01、第一層高度z1和與之對應的中性風速平均值確定，則與lnzi呈線性相關關系。采用最小二乘法擬合散點(lnzi,)得到斜率a，即：

考慮到參照風速也可選為其他層高度風速，因此分別計算第i層高度風速作為參照風速時與之對應的粗糙度z0i，再對所有z0i求平均值得到最終所求粗糙度z0，即：

1.2 中性風速選取方式

根據測風塔實測數據，可得到不同高度處風速與湍流強度相關曲線，中性溫度層結對湍流發生發展基本沒有影響[5]，即當湍流強度基本保持不變時，認為熱效應的影響消失，大氣穩定度可考慮為中性，此時把湍流強度開始變穩定的風速定義為中性大氣強風閾值[7]，把湍流強度基本穩定時對應的風速區間本文在此稱為中性風速組。

運用式(2)～(7)求取粗糙度的關鍵點在于各層高度中性風速的選取。測風塔各層高度(其中，zmax≤100 m；zmin建議取20～30 m，若測風塔周邊地表環境非常簡單，對10 m高度測風數據影響程度很小，則zmin可取10 m)的中性風速確定方式為：

1)選取一個或多個完整年時段的10 min實測數據并剔除其中無效數據得到測風塔有效數據，各高度有效數據完整率均需超過90%。

2)根據測風塔有效數據，繪制不同高度處風速與湍流強度曲線，查找相鄰風速段平均湍流強度變化率最先滿足|ε|≤5%時對應的風速值，該值即為中性大氣強風閾值Ur，同時參考我國法規中推薦的修訂帕斯奎爾分類法[8]，規定各層高度的Ur不小于6 m/s。為保證中性風速組內對應的湍流強度基本穩定，需在u≥Ur范圍內進一步剔除湍流強度變化率|ε|超過10%對應的風速段，此外，考慮到較高風速段(尤其是12～15 m/s以上)因數據量過少可能會造成湍流強度計算結果存在一定不確定性，由此引起的湍流強度變化率過大往往不具有代表性，因此，本文對于風速出現頻率小于0.5%的風速段不管湍流強度變化率是否超過10%均予以保留。綜上，得到的風速區間即為某層高度的初選中性風速，對應的測風時段即為初選中性風速時段。

3)為保證各高度中性風速測風時段同期，對步驟2)得到的各層高度初選中性風速的測風時段求交集得到共同的測風時段，該測風時段內各高度的風速即為符合中性大氣層結條件下的風速組。

2 實例分析

2.1 粗糙度取值對風能資源評估的影響

一般在風電場宏觀選址階段，規劃區域內基本上僅有少量已建測風塔，有的甚至沒有或在附近區域才有已建測風塔，這就很難直接反映整個規劃風電場區域的風能資源。對于平坦地形，WAsP軟件是目前比較常用的風能資源評估軟件，利用WAsP軟件推算規劃區域風能資源的主要過程為：①對場內或附近測風塔資料進行整理分析；②提取規劃風電場及附近區域地形地貌資料；③根據地形地貌資料繪制粗糙度線，并給粗糙度線兩側的地表粗糙度賦值；④將整理分析后的測風數據文件和處理后的地形圖文件導入軟件中進行相關計算，即可得出規劃區域風能資源分布。可見，地表粗糙度取值對風能資源的評估結果至關重要。為了說明此問題，本文選用某一平原區域來分析粗糙度取值對風能資源評估的影響。

2.1.1 研究區域概況

研究區域地面高程為3～5 m，地形平坦，基本地貌特征如圖1所示，整個區域由北至南地表粗糙度逐漸增大，風速由沿海向內陸地區呈衰減趨勢。

圖1 研究區域地形地貌及測風塔分布圖

區域①為大片魚塘養殖區，區域②為以農田為主的開闊地貌，區域③為有一些房屋和樹木的田野，區域④為有較多房屋和樹木的田野，區域⑤為有很多房屋接近市郊的田野，其余區域⑥均為有密集建筑群的小城鎮。現階段在區域②～⑤這四個區域內均建有一座測風塔，各測風塔基本信息如表1所示。

表1 測風塔基本信息表

2.1.2 粗糙度取值影響分析

對于平坦地形，可依照文獻[9]粗糙度劃分方法，區域①粗糙度選擇為0.02 m，區域②粗糙度選擇在0.03～0.10 m，區域③粗糙度選擇在0.05～0.20m，區域④粗糙度選擇在0.10～0.30 m，區域⑤粗糙度選擇在0.15～0.40 m，區域⑥粗糙度均選擇為0.6 m。

在上述粗糙度值劃分方法區間內，采用兩種設置方案。方案一：選擇區域②粗糙度值為0.08、區域③粗糙度值為0.1、區域④粗糙度值為0.12，區域⑤粗糙度值為0.15；方案二：區域②粗糙度值為0.03、區域③粗糙度值為0.15、區域④粗糙度值為0.2、區域⑤粗糙度值為0.4。以1#測風塔100 m高度實測資料作為輸入條件，模擬出的整個區域風速分布規律如圖2～圖3所示。

圖2 方案一的風速等值線分布

圖3 方案二的風速等值線分布

由圖2～圖3可知，隨著粗糙度取值的不同，整個區域的風能資源計算結果差別也較大。方案一情況下(即各區域的粗糙度相對變化值較小時)，區域③模擬風速變化范圍在6.1～6.3 m/s，自北向南風速衰減了0.2 m/s；區域④模擬風速變化范圍在6.0～6.15 m/s，自北向南風速衰減了0.15 m/s。方案二情況下(即各區域的粗糙度相對變化值較大時)，區域③模擬風速變化范圍在5.7～6.3 m/s，自北向南風速衰減了0.6 m/s；區域④模擬風速變化范圍在5.6～5.9 m/s，自北向南風速衰減了0.3 m/s。可見，隨著區域間粗糙度相對變化值的增大，區域③和區域④的風速變化范圍增大，整個區域自北向南風速衰減速率加快。

2.2 本文推薦方法計算結果分析

繼續以上述研究區域為例，采用本文提出的新的地表粗糙度計算方法來推求各局部區域的粗糙度值，并通過模擬結果來驗證本文推薦方法的可行性和可靠性。

2.2.1 中性風速組的選取

本文收集到的4個測風塔基本信息如表1所示，各測風塔數據均為10 min實測數據，其中1#和2#測風塔的測風年時段同期，3#和4#測風塔的測風年時段同期。采用1.2節提出的中性風速組選取方式來確定各測風塔的中性風速組。

以2#測風塔100 m和30 m兩個高度為例，100 m和30 m高度中性風速組初選過程如圖4所示。

圖4 2#塔100 m 和30 m高度中性風速組初選過程

由圖4(a)可知，100 m高度中性大氣強風閾值為8 m/s；在8～14 m/s風速段內各級風速出現頻率均大于0.5%，且該范圍內12 m/s風速段湍流強度和13 m/s風速段湍流強度的變化率超過10%，因此剔除風速區間[12 m/s, 12 m/s]后，得到2#塔100 m高度初選中性風速組為[8 m/s,12 m/s)∪(13 m/s, 20.6 m/s]。

由圖4(b)可知，由于30 m高度湍流強度變化率最先滿足|ε|≤5%時對應的風速值小于6 m/s，因此，30 m高度中性大氣強風閾值直接取為6 m/s；在6～12 m/s風速段內各級風速出現頻率均大于0.5%，且該范圍內沒有湍流強度變化率超過10%的風速段，因此，2#塔30 m高度初選中性風速組為u≥6 m/s。

以此類推，2#塔70 m高度初選中性風 速 組 為 [7 m/s, 11 m/s)∪(12 m/s, 19.3 m/s]；90 m高度初選中性風速組為[7 m/s, 12 m/s)∪(13 m/s, 19.8 m/s]。

在上述得到的各高度初選中性風速組的基礎上進一步篩選出同期時段的數據即為2#塔最終所求中性風速組。

同理，其余1#、3#和4#三座測風塔的中性風速組分別為：

1#塔30 m高度初選中性風速組為u≥6 m/s，50 m～90 m高度初選中性風速組均為u≥7 m/s，100 m高度初選中性風速段組u≥8 m/s，以此為基礎篩選出各高度同期時段數據即為1#塔終選中性風速組。

3#塔30～85 m高度初選中性風速組均為u≥6 m/s，100 m高度初選中性風速組為[7 m/s,11 m/s)∪(12 m/s, 19.9 m/s]，以此為基礎篩選出各高度同期時段數據即為3#塔終選中性風速組。

4#塔30 m和50 m高度初選中性風速組均為u≥6 m/s，70m高度初選中性風速組 為 [6 m/s, 9 m/s)∪ (10 m/s, 16.9 m/s]，85 m高度初選中性風速組為[6 m/s, 10 m/s)∪(12 m/s, 17.4 m/s]，100m高度初選中性風速組為[7 m/s, 11 m/s)∪(12 m/s, 18.1 m/s]，以此為基礎篩選出各高度同期時段數據即為4#塔終選中性風速組。

1#～4#四座塔各高度中性風速組的平均值如表2所示。

表2 各測風塔各高度中性風速平均值表

2.2.2 粗糙度計算結果

以表2中 4個測風塔中性風速平均值計算結果為基礎，采用最小二乘法分別擬合4個塔的散點(lnzi,)，各測風塔擬合結果圖如圖5所示。

圖5 各測風塔散點(lnzi,)擬合結果圖

根據圖5，1#塔擬合直線斜率a=1.304 3，結合式(7)和表2結果求得1#塔對應的區域②粗糙度值為0.040 m。以此類推，分別求得區域③粗糙度值為0.095 m，區域④粗糙度值為0.216 m，區域⑤粗糙度值為0.347 m。

以上各區域的粗糙度計算值均與文獻[9]粗糙度劃分方法結果基本吻合。

2.2.3 風資源模擬結果分析

考慮到1#和2#測風塔測風年時段同期，3#和4#測風塔測風年時段同期，因此，本文分別以1#和3#塔的100 m高度實測值作為輸入條件，采用WAsP軟件分別計算對比在方案一、方案二以及本文推薦方法這3種粗糙度取值方法下2#和4#塔位置處的計算結果。各方案下，1#塔推求2#塔對比結果如表3所示，3#塔推求4#塔對比結果如表4所示。

表3 1#塔推求2#塔的計算結果分析表

表4 3#塔推求4#塔的計算結果分析表

根據表3，當采用方案一時，2#塔計算風速與實測值相差0.11 m/s，風速誤差為1.82%，風功率密度誤差為5.86%；當采用方案二時，2#塔計算風速與實測值相差-0.13 m/s，風速誤差達-2.15%，風功率密度誤差達-8.11%；當采用本文推薦的方法時，2#塔計算風速與實測值僅相差-0.03 m/s，風速誤差僅為-0.49%，風功率密度誤差僅為-2.25%。

根據表4，當采用方案一時，4#塔計算風速與實測值相差0.24 m/s，風速誤差達4.31%，風功率密度誤差達14.04%；當采用方案二時，4#塔計算風速與實測值相差-0.07 m/s，風速誤差為-1.26%，風功率密度誤差為-4.49%；當采用本文推薦的方法時，4#塔計算風速與實測值僅相差-0.02 m/s，風速誤差僅為-0.36%，風功率密度誤差僅為-1.69%。

綜上，采用本文推薦方法計算的粗糙度值準確性較高，風能資源模擬結果更加接近真實值。

3 結語

地表粗糙度是影響陸上風電場風能資源評估的一項重要因素。以往在實際工程中，一般常用文獻[9]平坦地形粗糙度劃分法選取粗糙度值，但該方法是根據定性描述的地形特點來劃分粗糙度值，易受個人主觀影響，粗糙度取值存在較大不確定性，影響風能資源評估準確性。鑒于此，本文提出一種僅以測風塔實測數據為基礎計算地表粗糙度值的方法，通過實例驗證，采用新方法計算的地表粗糙度值準確性高，以此模擬的風能資源更加接近真實值，較大程度地降低了粗糙度取值的不確定性，而且計算所需基礎資料也簡單易獲取。