RWEQ模型中風(fēng)因子的改進(jìn)及應(yīng)用

鞏國麗， 劉 鑫， 要 玲， 任麗霞， 王 敏

(1.山西能源學(xué)院，030600，山西榆次；2.生態(tài)環(huán)境部信息中心，100029，北京；3.北京維盛沃德科技有限公司，100029，北京)

土壤風(fēng)蝕是我國北方地區(qū)土地退化的主要原因。為了更好地防治區(qū)域水土流失，需要定量估算區(qū)域風(fēng)蝕狀況及其原因。RWEQ模型是近幾年應(yīng)用比較普遍的能適應(yīng)大區(qū)域定量估算風(fēng)蝕量的模型。該模型是基于大量野外實驗的一種經(jīng)驗?zāi)Ｐ蚚1]，在實際測定風(fēng)力導(dǎo)致的土壤侵蝕量以及當(dāng)?shù)氐臍庀?、地表植被、土壤濕度、地表的結(jié)皮和地表的可蝕性等因子的基礎(chǔ)上得出的一個經(jīng)驗方程。風(fēng)是產(chǎn)生風(fēng)蝕的重要氣候驅(qū)動力，也是RWEQ模型中最重要的輸入數(shù)據(jù)，風(fēng)速輸入數(shù)據(jù)的質(zhì)量對土壤風(fēng)蝕量的估算結(jié)果影響較大。實際的風(fēng)速時刻在變化，模型中代入不同時段的平均風(fēng)速會產(chǎn)生不同的模擬結(jié)果，如日均風(fēng)速和月均風(fēng)速本身存在差異，風(fēng)因子的計算結(jié)果偏差會很大。此外，代入平均風(fēng)速會使得模型計算所得的風(fēng)蝕量偏小，因為風(fēng)蝕量的產(chǎn)生前提是風(fēng)速大于臨界起沙風(fēng)速，而平均風(fēng)速可能將瞬間的大于臨界起沙的風(fēng)速過濾掉。為此模型對計算所需風(fēng)速的時間分辨率要求較高，要求有不少于500個的風(fēng)速值用以估算模型中的風(fēng)因子，但除非點位觀測，大尺度的風(fēng)因子值的模擬很難獲得高分辨率的風(fēng)速值。為解決這一問題，有學(xué)者利用分形等方式進(jìn)行風(fēng)速時間序列的插值研究[2]。模型中也內(nèi)嵌一個風(fēng)速生成器，該風(fēng)速生成器的輸入?yún)?shù)為韋伯系數(shù)，且很多學(xué)者也提到可以采用Weibull分布參數(shù)擬合估算風(fēng)速[3-4]。但也有許多研究[5]表明該方法并不能很好地擬合我國近地面風(fēng)速。郭中領(lǐng)[6]利用日均、日最大風(fēng)速擬合法來進(jìn)行模擬，根據(jù)與實測風(fēng)速計算結(jié)果比較，二者存在一定的關(guān)系，但其擬合精度仍然需要提高，為此應(yīng)選擇合適的降尺度方法提高日均風(fēng)速時間分辨率，進(jìn)而更準(zhǔn)確地估算我國風(fēng)蝕量。

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

1.1 RWEQ模型

該模型可用于定量評估土壤風(fēng)蝕量(Fryerar D W et al,1998)，一般用土壤風(fēng)蝕模數(shù)Ls(kg/m2)表示，即單位面積的沙通量。該模型是一種經(jīng)驗?zāi)Ｐ?，與氣候因子W、土壤可蝕性E、土壤結(jié)皮S、土壤糙度K′和綜合植被因子C有關(guān)。

Ls=Qx/x；

(1)

(2)

式中：x為地塊長度，m；Qx為地塊長度x處的沙通量，kg/m；Qmax為風(fēng)力的最大輸沙能力，kg/m；s為關(guān)鍵地塊長度，m。

1.2 風(fēng)速觀測數(shù)據(jù)收集與處理

1)風(fēng)速觀測數(shù)據(jù)的收集。筆者利用國家氣象信息中心的中國氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.cn/)下載的2001—2010年內(nèi)蒙古自治區(qū)近10年數(shù)據(jù)量充足且DEM落差較小的45個國家氣象站點的一日四時次的實測風(fēng)速數(shù)據(jù)來進(jìn)行韋伯分布參數(shù)擬合，并將計算結(jié)果與實測風(fēng)速進(jìn)行比較來驗證我國北方地區(qū)的韋伯分布函數(shù)對風(fēng)速的擬合結(jié)果。由于風(fēng)力侵蝕多發(fā)生于春季，本文利用下載的我國北方384個國家氣象站點的3月份的風(fēng)速值代入RWEQ模型計算出風(fēng)因子值，用來分析我國北方的風(fēng)蝕驅(qū)動力的分布基本特征。

2)風(fēng)速觀測數(shù)據(jù)的處理。中國氣象數(shù)據(jù)網(wǎng)下載的風(fēng)速值為標(biāo)準(zhǔn)高度10 m下的監(jiān)測值，而模型中所需的風(fēng)速值為2 m高度下的近地面風(fēng)速，這就需要在計算前先進(jìn)行風(fēng)速值的高度轉(zhuǎn)換。有研究者利用空氣動力學(xué)粗糙度[7]來估算風(fēng)速隨高度的變化，模型手冊中提到利用七分之一定理[8]來轉(zhuǎn)換不同高度的風(fēng)速值，筆者直接利用Eiliot的轉(zhuǎn)換方法：

u2=u1(g2/g1)1/7。

(3)

式中：u1和u2分別為在g1和g2高度處的風(fēng)速，m/s。

1.3 風(fēng)因子估算方法對比與驗證

1)RWEQ模型中的風(fēng)因子值估算。筆者利用國家氣象信息中心的中國氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.cn/)下載的2001—2010年內(nèi)蒙古自治區(qū)45個國家氣象站點的一日四時次的實測風(fēng)速數(shù)據(jù)。由于風(fēng)力侵蝕多發(fā)生于春季，本文選用下載的3月份的風(fēng)速值分別代入模型和利用韋伯分布函數(shù)計算出風(fēng)因子值，將計算結(jié)果進(jìn)行比較來驗證我國北方地區(qū)的韋伯分布函數(shù)對風(fēng)速的擬合結(jié)果。下載的點狀風(fēng)速數(shù)據(jù)利用Arcgis軟件中的克里金插值法插值為面狀數(shù)據(jù)，根據(jù)其他氣象站點的風(fēng)速數(shù)據(jù)分析，插值結(jié)果準(zhǔn)確度較高。

研究表明風(fēng)蝕量與風(fēng)速的三次方存在相關(guān)關(guān)系[9-11]。RWEQ中將風(fēng)因子值表達(dá)為：

(4)

式中：F為風(fēng)蝕力，m3/s3；F2為2 m處的風(fēng)速，F(xiàn)t為臨界起沙風(fēng)速，部分地區(qū)由于地表情況不同，臨界起沙風(fēng)速不同，如塔克拉瑪干沙漠等區(qū)域的臨界起沙風(fēng)速會偏小，但大部分區(qū)域一般為5 m/s。

2)二參數(shù)Weibull分布模型。利用最小二乘法來求取韋伯分布函數(shù)的二參數(shù)。

韋伯分布的概率分布函數(shù)[3-4]為：

(5)

式中：P(u≤ua)為風(fēng)速u小于給定風(fēng)速ua的概率；C為尺度參數(shù)；X為形狀參數(shù)。

筆者利用最小二乘法來求取韋伯分布函數(shù)的二參數(shù)。

將概率分布函數(shù)兩端取兩次對數(shù)后，公式變?yōu)椋?/p>

ln(-ln(1-P))=Xlnua-XlnC。

(6)

為此將ln(-ln(1-P))看做為B，lnua看做A，XlnC看做b時，上式便可寫為如下形式：

B=XA-b。

(7)

利用最小二乘法擬合實測風(fēng)速數(shù)據(jù)便可求得C、X值。

在風(fēng)速經(jīng)過高度轉(zhuǎn)換后，利用最小二乘法求得韋伯分布參數(shù)，并利用參數(shù)得出擬合風(fēng)速：

(8)

在此根據(jù)模型中的風(fēng)速值數(shù)量要求，設(shè)500個間隔概率數(shù)(0.002，0.004，…，0.999)。其中將實測數(shù)據(jù)的靜風(fēng)的概率與利用韋伯分布函數(shù)擬合的累積概率比較，當(dāng)累積概率小于實測數(shù)據(jù)的靜風(fēng)概率時，將擬合數(shù)據(jù)設(shè)為0。

3)日均、日最大風(fēng)速擬合法。郭中領(lǐng)[6]利用日均風(fēng)速和日最大風(fēng)速分別提出了風(fēng)速的降尺度方法，并在此方法基礎(chǔ)上，分別計算風(fēng)蝕力，求得擬合風(fēng)速計算結(jié)果風(fēng)蝕力與實測風(fēng)速計算結(jié)果風(fēng)蝕力的關(guān)系式。

(9)

(10)

式中：Wn為擬合的風(fēng)速；Wavg為日均風(fēng)速；Wmax為日最大風(fēng)速，為了滿足RWEQ手冊中最低風(fēng)速數(shù)據(jù)為500的要求，當(dāng)風(fēng)蝕力計算時間尺度為月時，n取1～12，由此獲取每小時風(fēng)速值。Wavg和Wmax來源于國家氣象信息中心的中國氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.cn/)下載的地面氣候資料日值數(shù)據(jù)。

1.4 改進(jìn)風(fēng)因子在北方地區(qū)的模型誤差分析

利用誤差分析來驗證二參數(shù)韋伯分布函數(shù)是否適用于我國近地層的風(fēng)速擬合，設(shè)用實測數(shù)據(jù)求得的結(jié)果為fr，利用weibull分布函數(shù)擬合數(shù)據(jù)求得的結(jié)果為fw，則兩者的相對誤差E由下列公式求出：

(11)

2 結(jié)果與分析

2.1 結(jié)果驗證

結(jié)果表明利用二參數(shù)韋伯分布函數(shù)擬合風(fēng)速所求的風(fēng)蝕力與實測的風(fēng)速所求的風(fēng)蝕力存在誤差，相對誤差值從0.19到0.95，平均相對誤差值達(dá)到0.52。回歸分析結(jié)果表明二者相關(guān)系數(shù)達(dá)到了0.97，顯著性水平P<0.01。由于風(fēng)蝕力的計算結(jié)果是在風(fēng)速大于臨界起沙風(fēng)速(一般設(shè)為5 m/s)的情況下計算的，為此推測二參數(shù)韋伯分布函數(shù)沒能很好地擬合5 m/s以上的風(fēng)速值，使得計算結(jié)果偏差較大。表1顯示5 m/s以上的風(fēng)速值的擬合結(jié)果誤差均值達(dá)到0.09，二者的相關(guān)系數(shù)達(dá)到0.82，顯著性水平P<0.01。表示雖然擬合結(jié)果存在誤差，但是兩者具有很好的相關(guān)性。

表1 韋伯風(fēng)速擬合結(jié)果誤差分析

圖1為由實測的風(fēng)速值和由二參數(shù)韋伯分布函數(shù)擬合風(fēng)速計算所得的風(fēng)蝕力以及平均風(fēng)速(>5 m/s)的結(jié)果比較。結(jié)果顯示擬合所得風(fēng)蝕力被低估，斜率達(dá)到0.65。同樣由于風(fēng)蝕力是由大于臨界起沙的風(fēng)速值計算而得，為此推算擬合所得的5 m/s以上的風(fēng)速值結(jié)果偏小，對比結(jié)果也顯示高值風(fēng)速區(qū)(>5 m/s)的擬合結(jié)果偏小，斜率為0.74。且Donk等[5]的研究結(jié)果也表明利用韋伯分布函數(shù)擬合的風(fēng)速數(shù)據(jù)可能會被低估，尤其高值風(fēng)速。

圖1 實測與擬合風(fēng)速、月均風(fēng)蝕力計算結(jié)果對比Fig.1 Comparison of the measured results with the calculated results of the wind speed and monthly average wind erosion force

國家氣象信息中心的中國氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.cn/)能夠普遍下載的為日均和日最大風(fēng)速數(shù)據(jù)，很難滿足風(fēng)蝕力估算對風(fēng)速的質(zhì)量要求。而二參數(shù)韋伯分布函數(shù)擬合結(jié)果存在風(fēng)速值的低估，為此需要選用其他辦法來解決此問題，一些經(jīng)驗的方法常用于對日均風(fēng)速降尺度[12-13]，但這些方法都不能模擬可蝕性的風(fēng)速數(shù)據(jù)。本文利用郭中領(lǐng)[6]的日均風(fēng)速和日最大風(fēng)速擬合法，求得擬合風(fēng)速計算結(jié)果風(fēng)蝕力與實測風(fēng)速計算結(jié)果風(fēng)蝕力的關(guān)系式，并判定其擬合的精度。

由于中國北方氣象站點下載數(shù)據(jù)多為日均風(fēng)速，最大風(fēng)速值缺測數(shù)較多，在此情況下本文估算風(fēng)蝕量多使用了日均風(fēng)速，在日最大風(fēng)速滿足條件時，也參與了風(fēng)速擬合。為準(zhǔn)確應(yīng)用于模型，提高風(fēng)因子估算結(jié)果的準(zhǔn)確性，根據(jù)擬合風(fēng)速計算結(jié)果風(fēng)蝕力與實測風(fēng)速計算結(jié)果風(fēng)蝕力的關(guān)系式(如圖2中擬合所得的方程)反推得出相對準(zhǔn)確的擬合風(fēng)速的風(fēng)蝕力，即在實測風(fēng)速值不夠的情況下可利用模擬風(fēng)速值代入類似圖2的擬合方程進(jìn)而計算風(fēng)蝕力，進(jìn)一步應(yīng)用于RWEQ估算風(fēng)蝕量。

圖2 實測風(fēng)蝕力與模擬風(fēng)蝕力的相關(guān)關(guān)系Fig.2 Correlation between measured wind erosion force and simulated wind erosion force

2.2 中國北方風(fēng)因子時空分布與模擬

國家氣象信息中心的中國氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.cn/)中能夠下載的風(fēng)速數(shù)據(jù)為國家基準(zhǔn)站的風(fēng)速數(shù)據(jù)，而國家基準(zhǔn)站數(shù)量偏少，氣象觀測站點的密集度不夠，且分布不均，如很多沙漠區(qū)域分布的氣象站點數(shù)量較少，塔克拉瑪干沙漠腹地就只有塔中一個氣象站點，而該氣象站點的分布區(qū)正好是風(fēng)速較小區(qū)，這就影響了塔克拉瑪干沙漠區(qū)的風(fēng)蝕量估算精度。同時中國地域遼闊，地形復(fù)雜，地表土壤質(zhì)地等條件不一，不同時期的土地覆被情況也不同，即使是同一覆被條件下，不同土壤質(zhì)地區(qū)的臨界起沙風(fēng)速也不同[14-16]，筆者在利用RWEQ模型計算中，塔克拉瑪干沙漠區(qū)使用2.5 m/s的臨界起沙風(fēng)速，對于其他區(qū)域沒能很好地考慮不同土壤質(zhì)地和不同土地覆被區(qū)域的臨界起沙風(fēng)速，統(tǒng)一使用5 m/s，為此需要進(jìn)一步研究。

由于我國風(fēng)蝕多發(fā)生于植被尚未較好覆蓋，風(fēng)速較大的春季，為此選取3月份作為典型月，求取北方風(fēng)蝕區(qū)的半月風(fēng)因子值(圖3)。從圖可知，內(nèi)蒙古高原上各大沙地的分布區(qū)的風(fēng)因子值較高，東北及青藏高原高山地帶等地區(qū)的風(fēng)因子值也較高。

圖3 多年平均(2000—2010年)風(fēng)因子時空分布(以3月份為例)Fig.3 Temporal and spatial distribution of multi-year average (2000-2010) wind factor in March

3 結(jié)論

1)Weibull分布模型擬合所得高值風(fēng)速區(qū)(>5 m/s)的風(fēng)速值偏小，為此不能利用該模型對風(fēng)速的時間分辨率進(jìn)行插值進(jìn)而代入RWEQ模型計算風(fēng)蝕力。

2)日均、日最大風(fēng)速擬合法擬合效果也不盡如意，但該方法所得結(jié)果與實測數(shù)據(jù)所得結(jié)果具有很強(qiáng)的相關(guān)性，為此可依據(jù)模擬與實測的關(guān)系式反推風(fēng)因子。

3)內(nèi)蒙古高原上各大沙地分布區(qū)的風(fēng)因子值較高，東北及青藏高原高山地帶等地區(qū)的風(fēng)因子值也較高。

4 討論

1)反推出的風(fēng)因子可用于指導(dǎo)風(fēng)蝕防治，有利于找出風(fēng)蝕防治的關(guān)鍵區(qū)域和關(guān)鍵時間。

2)鑒于目前實測風(fēng)蝕量的方法多種多樣，有風(fēng)洞實驗法，同位素分析法，插釬法，風(fēng)蝕盤法等等，各方法不同所獲取的實測風(fēng)蝕量數(shù)據(jù)本身存在誤差，且可獲取的實測數(shù)據(jù)多為點狀，區(qū)域性的風(fēng)蝕數(shù)據(jù)較少，難以對利用反推的風(fēng)因子帶入RWEQ模型估算的風(fēng)蝕量與實測風(fēng)蝕量進(jìn)行對比驗證，為此還需進(jìn)一步研究。

3)不同土壤質(zhì)地，不同土地覆被條件下的臨界起沙風(fēng)速不同，大大尺度風(fēng)因子值估算過程中如何細(xì)化分區(qū)，獲取不同分區(qū)的臨界起沙風(fēng)速是需要進(jìn)一步考慮的。

4)處于山風(fēng)口區(qū)域的氣象站點觀測的風(fēng)速值很高，這往往不代表整個區(qū)域的情況，本文是在對山風(fēng)口氣象站點觀測數(shù)據(jù)刪除的情況下模擬所得結(jié)果，但實際地形復(fù)雜區(qū)域的風(fēng)向風(fēng)速變化很大，這就使得風(fēng)速插值難度大，有待進(jìn)一步研究。

5)國家氣象信息中心的中國氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.cn/)中能夠下載的風(fēng)速數(shù)據(jù)為國家基準(zhǔn)站的風(fēng)速數(shù)據(jù)，而國家基準(zhǔn)站數(shù)量偏少，氣象觀測站點的密集度不夠，有些觀測站點所設(shè)置的區(qū)域正好是風(fēng)速較高或者較低區(qū)，為此如何獲取更準(zhǔn)確的觀測數(shù)據(jù)是值得進(jìn)一步研究的。