莊會(huì)波，季 妤，高振勇，譚秀翠

（1.山東省水文局，濟(jì)南250002；2.山東農(nóng)業(yè)大學(xué)水利土木工程學(xué)院，山東泰安271018）

植被作為陸地生態(tài)環(huán)境中的重要組成部分，具有明顯的時(shí)空變化特點(diǎn)，其覆蓋變化可以反映區(qū)域內(nèi)自然因素與人為因素對(duì)植被的影響效果［1，2］。植被變化改變流域下墊面條件，并通過影響降雨截留、蒸散發(fā)、降水入滲、地下水補(bǔ)給等水文過程，進(jìn)而直接影響徑流的形成與演變［3，4］，因此，分析植被覆蓋的時(shí)空變化規(guī)律，及其對(duì)徑流過程的影響，有利于闡述流域下墊面變化對(duì)水文循環(huán)過程的影響機(jī)制。畢早瑩等［5］基于Budyko 理論定量分析窟野河流域植被變化對(duì)徑流的影響，研究結(jié)果表明植被變化對(duì)徑流的貢獻(xiàn)率達(dá)到84.9%，而氣候變化和人類活動(dòng)影響比較小。李慧赟等［6］借助改進(jìn)的水文模型模擬植被變化后的徑流過程，在Crawford 流域有25%的面積植樹造林情況下，流域年徑流量減少32.4 mm，由植被變化引起的徑流變化量占63%，由氣候變化引起的徑流變化量占37%。一般認(rèn)為植物覆蓋通過增加滲透到土壤基質(zhì)水分，進(jìn)而減少地表徑流［7］，但這并不適應(yīng)所有地區(qū)，所以，植被覆蓋對(duì)徑流產(chǎn)生怎樣的水文效應(yīng)還需要進(jìn)一步研究。丁婧祎等［8］研究顯示NDVI與徑流的相關(guān)系數(shù)在0.04～0.59 之間，13 個(gè)子流域NDVI與徑流間的正負(fù)相關(guān)性并不一致。在不同的時(shí)空尺度下，植被覆蓋對(duì)徑流的影響規(guī)律也存在差異。李春暉等［9］研究黃河流域降水、植被對(duì)徑流影響，結(jié)果表明，在年內(nèi)，NDVI變化與降水、徑流變化呈明顯的正相關(guān)，在年際，NDVI波動(dòng)趨勢(shì)與降水具有一定的相關(guān)性，但與徑流、徑流系數(shù)變化的關(guān)系相對(duì)復(fù)雜，不具有明顯的相關(guān)性。董喆等［10］研究表明，徑流量受降水與NDVI的共同影響，其中7-10月的降水與NDVI對(duì)徑流的影響顯著。史曉亮等［11］研究發(fā)現(xiàn)流域尺度上植被NDVI與徑流不具有明顯的相關(guān)性，但從像元尺度來看，植被NDVI和徑流的正相關(guān)和負(fù)相關(guān)共存，流域不同空間位置的植被變化與徑流的關(guān)系并不一致。

歸一化植被指數(shù)NDVI（Normalized Difference Vegetation Index）作為植被生長(zhǎng)狀況及植被覆蓋度的最佳指示因子，被認(rèn)為是監(jiān)測(cè)植被和生態(tài)環(huán)境變化的最有效指標(biāo)［12］。文中采用歸一化植被指數(shù)NDVI為植被覆蓋表征指標(biāo)，分析山東省植被覆蓋的時(shí)空變化規(guī)律，與水文要素進(jìn)行相關(guān)分析，并采用Copula函數(shù)構(gòu)建徑流與降水量、NDVI的二維聯(lián)合分布，系統(tǒng)分析徑流與影響因素之間的概率分布關(guān)系，研究結(jié)果可為水資源合理開發(fā)利用及下墊面因素影響分析提供技術(shù)指導(dǎo)。

1 研究區(qū)域

山東省位于中國(guó)東部沿海、黃河下游，境內(nèi)中部山地突起，西南、西北低洼平坦，東部緩丘起伏，山東省地形分布見圖1。山東水系比較發(fā)達(dá)，自然河流的平均密度每平方公里在0.7 km以上。湖泊集中分布在魯中南山丘區(qū)與魯西南平原之間的魯西湖帶。氣候?qū)倥瘻貛Ъ撅L(fēng)氣候類型，年平均氣溫11～14 ℃，降水季節(jié)分布很不均衡，全年降水量有60%～70%集中于夏季。由2000-2018年《山東省水資源公報(bào)》數(shù)據(jù)，年均降水量675.99 mm，徑流深118.53 mm，各年降水量與徑流深分布過程線見圖2。

圖1 山東省地形圖Fig.1 Topographic aap of Shandong Province

圖2 2000-2018年山東省降水量與徑流深分布Fig.2 The distribution of precipitation and runoff depth in Shandong Province from 2000 to 2018

山東省土地利用現(xiàn)狀數(shù)據(jù)由資源環(huán)境科學(xué)與數(shù)據(jù)中心（http：//www.resdc.cn/）獲取，數(shù)據(jù)格式為1 km 的柵格數(shù)據(jù)。2015年，山東省土地利用類型以耕地為主，占總面積的65%，其次為建設(shè)用地，占15%，草地、林地、水域分別占比8%、6%和4%，其空間分布見圖3。

圖3 2015年山東省土地利用數(shù)據(jù)Fig.3 Land use data of Shandong Province in 2015

歸一化植被指數(shù)NDVI能反映地表植被覆蓋變化狀況，是評(píng)價(jià)植被生長(zhǎng)狀況的重要指標(biāo)［13］。文中采用的NDVI數(shù)據(jù)來源于中國(guó)年度植被指數(shù)（NDVI）空間分布數(shù)據(jù)集［14］。1998-2018年山東省平均NDVI為0.71，植被指數(shù)NDVI空間分布見圖4。

圖4 1998-2018年山東省植被指數(shù)NDVI空間分布Fig.4 Spatial distribution of vegetation index NDVI in Shandong Province from 1998 to 2018

2 研究方法

2.1 NDVI趨勢(shì)分析

一元線性回歸趨勢(shì)分析法用于計(jì)算每個(gè)柵格的變化趨勢(shì)，反映不同時(shí)期NDVI變化趨勢(shì)的空間特征［15］，計(jì)算公式如下所示：

式中：i為年序列號(hào)；n為總年數(shù)，取值21（1998-2018年）；NDVIi為第i年NDVI值；slope為各柵格NDVI變化趨勢(shì)的斜率。

2.2 突變檢測(cè)

曼-肯德爾（Mann-Kendall）法是一種非參數(shù)統(tǒng)計(jì)檢驗(yàn)方法［16］，被廣泛應(yīng)用于氣候［17］、水文［18］、植被［19］等變量的突變分析中。

對(duì)于有n個(gè)樣本量的時(shí)間序列x，構(gòu)造一個(gè)秩序列：

在時(shí)間序列隨機(jī)獨(dú)立的假定下，定義統(tǒng)計(jì)量：

式中：UF1=0，E（sk）和var（sk）是sk的均值和方差，在x1，x2，…，xn相互獨(dú)立，且有相同連續(xù)分布時(shí)。

給定顯著性水平α，若|UFi|>Uα，則表明序列存在明顯的趨勢(shì)變化。該方法計(jì)算簡(jiǎn)便，可以明確突變開始的時(shí)間，并指出突變區(qū)域，因此是常用的突變檢測(cè)方法。

按時(shí)間序列的逆序xn，xn-1，…，x1，再重復(fù)上述過程，同時(shí)使UBk=-UFk（k=n，n-1，…，1），UB1=0。

2.3 Copula聯(lián)合分布

Copula 函數(shù)是定義域?yàn)椋?，1］均勻分布的多維聯(lián)合分布函數(shù)，它可以將多個(gè)隨機(jī)變量的邊際分布連接起來構(gòu)造聯(lián)合分布［20］。二維聯(lián)合經(jīng)驗(yàn)分布，計(jì)算公式如下：

式中：P為X≤xi，Y≤yi的二維聯(lián)合概率值；Nmn為X≤xi，Y≤yi的序號(hào)；N為總的數(shù)據(jù)對(duì)數(shù)目。

對(duì)應(yīng)的二維聯(lián)合重現(xiàn)期計(jì)算公式為：

在水文計(jì)算中，二維聯(lián)合理論分布函數(shù)，一般采用最多的是Archimedean Copula 函數(shù)，其中Clayton Copula、Gumbel Copu?la、Frank Copula最為常用。

Clayton Copula函數(shù)：

Gumbel Copula函數(shù)：

Frank Copula函數(shù)：

為實(shí)現(xiàn)對(duì)Copula函數(shù)的擬合優(yōu)度檢驗(yàn)，采用、AIC信息準(zhǔn)則和離差平方和最小準(zhǔn)則（OLS）來評(píng)價(jià)Copula 方法的有效性［21］，選擇AIC值越小、OLS值越小的Copula函數(shù)構(gòu)建二維聯(lián)合分布。

AIC信息準(zhǔn)則計(jì)算公式為：

式中：Pei、Pi分別為聯(lián)合經(jīng)驗(yàn)頻率和聯(lián)合理論頻率；i為樣本序號(hào)；n為樣本容量；m為模型參數(shù)的個(gè)數(shù)。

AIC值越小，表明Copula函數(shù)的擬合效果越好。

OLS的計(jì)算公式為：

3 結(jié)果分析

3.1 NDVI特征分析

采用一元線性回歸趨勢(shì)分析法，計(jì)算1998-2018年山東省NDVI空間變化趨勢(shì)，計(jì)算結(jié)果見圖5，可以看出，大部分區(qū)域NDVI呈增加趨勢(shì)，植被覆蓋良好，其占總面積的63%（slpoe>0.001），NDVI基本不變的區(qū)域占13%（0.001≥slpoe>-0.001），而NDVI呈減小趨勢(shì)的區(qū)域占24%（slpoe≤-0.001）。

圖5 1998-2018年山東省NDVI變化趨勢(shì)Fig.5 Trends of NDVI in Shandong Province from 1998 to 2018

由山東省NDVI突變檢測(cè)結(jié)果可以看出，見圖6，2004年之后，NDVI統(tǒng)計(jì)值UF>0，NDVI呈增加趨勢(shì)，但未超過顯著水平線|UF|

圖6 NDVI突變檢驗(yàn)Fig.6 The abruption tests of NDVI

由山東省不同時(shí)期土地利用數(shù)據(jù)可以看出，見圖7，2000年至2015年，土地利用結(jié)構(gòu)發(fā)生了一定變化，主要是由耕地、草地、未利用土地向建設(shè)用地轉(zhuǎn)化，其面積增大2.2%。不同土地利用區(qū)域的NDVI值也有較大差異，由表1 可以看出，NDVI值相對(duì)較低的區(qū)域集中在未利用地、水域和建設(shè)用地，其中僅建設(shè)用地的NDVI變化趨勢(shì)是輕度減少，主要受城市化擴(kuò)張影響。

表1 不同土地利用類型區(qū)域NDVI變化趨勢(shì)Tab.1 Variation trend of NDVI in different land use types

圖7 2000-2015年山東省土地利用類型數(shù)據(jù)分析Fig.7 Data Analysis of Land Use Types in Shandong Province from 2000 to 2015

3.2 水文要素特征分析

根據(jù)2000-2018年《山東省水資源公報(bào)》中降水、徑流深統(tǒng)計(jì)資料，進(jìn)行突變檢測(cè)。

由徑流深突變檢測(cè)結(jié)果可以看出，見圖8，徑流深|UF|

圖8 徑流深突變檢驗(yàn)Fig.8 The abruption tests of runoff depth

由降水量突變檢測(cè)結(jié)果可以看出，見圖9，與徑流深突變分析曲線分布比較類似，降水量|UF|

圖9 降水量突變檢驗(yàn)Fig.9 The abruption tests of precipitation

受土地利用類型與植被覆蓋度的雙重影響，2000-2015年，耕地、草地面積減少，但區(qū)域內(nèi)NDVI增大，建設(shè)用地面積增加，區(qū)域內(nèi)NDVI降低。而2000-2015年，降水量由607.3 mm，減少到575.7 mm，徑流深由106.9 mm，減少到53.8。受降水與下墊面因素共同影響，2000年到2015年，山東省徑流深減少50%。

3.3 NDVI水文響應(yīng)分析

為分析NDVI與徑流的變化關(guān)系，采用相同時(shí)間尺度下的統(tǒng)計(jì)數(shù)據(jù)，分析NDVI與徑流的變化關(guān)系。由圖10（a）可以看出，在年尺度下，徑流深與NDVI呈正相關(guān)，即隨著NDVI的增大，徑流深增大，但兩者的相關(guān)系數(shù)R2值較低，僅為0.24，而徑流深與降水量的相關(guān)程度較為密切，兩者的相關(guān)系數(shù)R2值為0.90。表明，在年尺度下，山東省的徑流深主要受降水量的影響。

由徑流深與NDVI、降水量三者的分布關(guān)系可以看出，見圖10（b），隨著NDVI的增大、降水量的增大，徑流深呈增加趨勢(shì)，徑流深最小值為33.9 mm，對(duì)應(yīng)的降水與NDVI均為最小值，分別為420.2 mm 和0.65，當(dāng)徑流深達(dá)到最大值222.9 mm 時(shí)，對(duì)應(yīng)的降水量為936.3 mm，NDVI為0.71。

圖10 2000-2018年山東省徑流深與降水、NDVI相關(guān)關(guān)系Fig.10 The correlation between runoff depth，precipitation and NDVI in Shandong Province from 2000 to 2018

3.4 二維Copula聯(lián)合分布

多變量聯(lián)合分析比單變量分析能更好地描述水文事件的內(nèi)在規(guī)律和分析各個(gè)特征屬性之間的相互關(guān)系［22］。因此，本文利用Copula 函數(shù)分別構(gòu)建徑流深與NDVI、降水量的二維聯(lián)合分布，該函數(shù)不需要假設(shè)變量是獨(dú)立的或者正態(tài)分布的或者它們具有相同的邊緣分布，可以用來描述水文變量之間的相關(guān)性結(jié)構(gòu)，能夠靈活地構(gòu)造邊緣分布為任意分布的水文變量聯(lián)合分布［20］。

通過擬合優(yōu)度檢驗(yàn)，計(jì)算結(jié)果如表2所示，通過比較可以看出，Gumbel Copula 函數(shù)的AIC與OLS值均最小，擬合精度較高。由Gumbel Copula 函數(shù)的理論聯(lián)合分布與經(jīng)驗(yàn)聯(lián)合分布圖可以看出，見圖11、圖12，Gumbel Copula 聯(lián)合分布的相關(guān)系數(shù)R2分別是0.86和0.97，表明選取Gumbel Copula函數(shù)構(gòu)建二維聯(lián)合分布是合理的。

表2 Copula函數(shù)的擬合優(yōu)度檢驗(yàn)Tab.2 The goodness-of-fit test for Copula function

圖11 NDVI與徑流深理論聯(lián)合分布與概率聯(lián)合分布Fig.11 The theory joint distribution and probability joint distribution of NDVI and runoff depth

圖12 降水量與徑流深理論聯(lián)合分布與概率聯(lián)合分布Fig.12 The theory joint distribution and probability joint distribution of precipitation and runoff depth

由徑流深與NDVI的二維聯(lián)合分布圖可以看出［見圖13（a）］，隨著徑流深與NDVI的增大，兩者的聯(lián)合概率值增大。由二維聯(lián)合分布等值線圖［見圖13（b）］，可以確定不同徑流深與NDVI組合的聯(lián)合概率值，當(dāng)徑流深一定時(shí)，隨著NDVI的增大，聯(lián)合概率值增大；當(dāng)NDVI一定時(shí)，隨著徑流深的增大，聯(lián)合概率值增大。

Fig.13 徑流深與NDVI二維聯(lián)合分布與等值線圖Fig.13 Two-dimensional joint distribution and contour map of runoff depth and NDVI

徑流深與降水的二維聯(lián)合分布反映了基本相同的變化規(guī)律，由圖14可以看出，隨著徑流深與降水量的增大，兩者的聯(lián)合概率值增大。根據(jù)等值線圖，可以確定不同徑流深與降水量組合的概率值。

Fig.14 徑流深與降水量二維聯(lián)合分布與等值線圖Fig.14 Two-dimensional joint distribution and contour map of runoff depth and precipitation

由徑流深與NDVI、降水量的聯(lián)合重現(xiàn)期及等值線圖，可以看出（見圖15、圖16），隨著各要素的增大，聯(lián)合重現(xiàn)期增大，但最大不超過20年，且隨著重現(xiàn)期的增大，徑流深與NDVI、降水量的取值范圍變小。聯(lián)合重現(xiàn)期是表示為“或”的聯(lián)合分布重現(xiàn)期，即徑流深大于等于特定值或是NDVI、降水量大于等于特定值時(shí)的重現(xiàn)期情況。

Fig.15 徑流深與NDVI聯(lián)合重現(xiàn)期與等值線圖Fig.15 Combined return period and contour map of runoff depth and NDVI

Fig.16 徑流深與降水量二維聯(lián)合分布與等值線圖Fig.16 Combined return period and contour map of runoff depth and precipitation

4 結(jié) 論

本文分析了1998-2018年山東省植被覆蓋時(shí)空變化規(guī)律及其對(duì)徑流的影響，基于Copula 函數(shù)構(gòu)建了徑流與NDVI、降水的二維聯(lián)合分布函數(shù)，通過以上研究得到如下結(jié)論。

（1）在空間上，山東省63%的面積NDVI增大，13%面積NDVI不變，24%的面積NDVI減小。

（2）在時(shí)間上，山東省NDVI呈增加趨勢(shì)，但增加趨勢(shì)不顯著。

（3）徑流深與NDVI、降水量呈正相關(guān)，在NDVI與降水量的共同影響下，2011年前后，徑流深減少30%。

（4）采用Gumbel Copula 函數(shù)構(gòu)建徑流深與NDVI、降水量的二維聯(lián)合分布，可以確定不同要素組合下的聯(lián)合概率值及重現(xiàn)期。

植被覆蓋變化受自然因素與人為因素共同影響，在空間上呈現(xiàn)差異性，其產(chǎn)生的水文效應(yīng)就更為復(fù)雜。如何在空間上，有效的分析植被覆蓋的水文響應(yīng)，遙感技術(shù)提供了良好的技術(shù)手段，其可以快速獲取大范圍、高精度的水文資料，如降水［23，24］、徑流［25］等，有效彌補(bǔ)了水文氣象站點(diǎn)布設(shè)不足的缺陷。因此，可在后續(xù)的研究中，將綜合運(yùn)用實(shí)測(cè)資料與遙感衛(wèi)星資料，提高數(shù)據(jù)分析的精度與可靠性，為水文響應(yīng)研究提供有效手段。 □