杜 勇，李建柱，牛凱杰，馮 平，胡慶芳，鄭彥辰

（1. 天津大學(xué) 水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300072；2. 永定河流域投資有限公司，北京 100193；3. 南京水利科學(xué)研究院 水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210029）

1 研究背景

植被作為陸地生態(tài)系統(tǒng)的重要組成部分，具有連接大氣、土壤、水文等生態(tài)要素的功能，在自然生態(tài)系統(tǒng)和人類(lèi)生產(chǎn)生活中發(fā)揮著不可替代的作用［1-3］。歸一化植被指數(shù)（Normalized Difference Vegeta?tion Index，NDVI）是公認(rèn)的反映大尺度地表植被覆蓋和生長(zhǎng)狀況的有效指標(biāo)，與光合有效輻射、綠葉密度、植被生產(chǎn)力以及累積生物量等均具有密切的關(guān)系，通常用來(lái)表征植被活動(dòng)的強(qiáng)弱［4-6］。近年來(lái)，由于氣候變化和人類(lèi)活動(dòng)的影響，區(qū)域乃至全球尺度的植被均呈現(xiàn)顯著增加趨勢(shì)［7-10］，植被變化不僅能影響氣候，而且可以通過(guò)改變植被蒸騰、冠層截流和土壤蒸發(fā)等水文過(guò)程而改變河川徑流，故能顯著影響流域水文循環(huán)過(guò)程［11-14］。因此，分析植被的時(shí)空演替特征，揭示在不同時(shí)空尺度下導(dǎo)致植被變化的驅(qū)動(dòng)因子的空間異質(zhì)規(guī)律，并解析以植被變化為重要特征的下墊面變化對(duì)徑流變化的貢獻(xiàn)及影響，可為應(yīng)對(duì)氣候變化和人類(lèi)活動(dòng)下生態(tài)系統(tǒng)的適應(yīng)能力提供理論依據(jù)，同時(shí)有助于深刻理解植被活動(dòng)對(duì)水循環(huán)過(guò)程的影響，并對(duì)完善流域水資源管理措施具有重要理論和實(shí)踐意義。

作為植物生長(zhǎng)發(fā)育所必需的外界環(huán)境因子，氣候要素（水分和溫度等）對(duì)植被活動(dòng)起重要作用，同時(shí)人類(lèi)活動(dòng)對(duì)植被活動(dòng)也具有不同程度的促進(jìn)或抑制作用［15-18］。以黃土高原為例，已有研究表明過(guò)去30年氣候變化對(duì)黃土高原的植被恢復(fù)有一定促進(jìn)，同時(shí)植被恢復(fù)工程等人類(lèi)活動(dòng)也有利于NDVI的增加，因此氣候變化和人類(lèi)活動(dòng)共同驅(qū)動(dòng)黃土高原的植被變化，但兩者的驅(qū)動(dòng)作用呈現(xiàn)出顯著的空間異質(zhì)性［2，19-20］。同時(shí)大量研究證實(shí)，植被覆蓋作為下墊面的重要性質(zhì)，其變化將導(dǎo)致水循環(huán)過(guò)程變化，影響河川徑流過(guò)程［21-23］。目前有關(guān)植被變化等研究已在黃土高原等生態(tài)脆弱區(qū)和人類(lèi)活動(dòng)密集區(qū)全面開(kāi)展，但這些研究多采用相關(guān)性分析等方法對(duì)植被變化的動(dòng)因進(jìn)行定性研究，關(guān)于不同時(shí)間尺度下各種驅(qū)動(dòng)因素影響的空間異質(zhì)性的定量解析還不夠深入［24-25］，且當(dāng)前關(guān)于NDVI 變化對(duì)天然徑流量影響的文獻(xiàn)較少，不能清晰地從源頭理解植被變化對(duì)河川徑流過(guò)程的影響。

本文以黃土高原上的典型流域永定河山區(qū)為代表，基于氣溫、降水、徑流以及長(zhǎng)時(shí)序NDVI 數(shù)據(jù)集，通過(guò)構(gòu)建地理加權(quán)回歸模型，研究流域內(nèi)NDVI 和降水、平均氣溫的回歸參數(shù)、殘差均值及其趨勢(shì)的空間異質(zhì)性，并進(jìn)一步探討以植被變化為突出特征的下墊面變化對(duì)流域典型斷面天然徑流的影響。

2 研究方法

2.1 NDVI 與影響因素的空間趨勢(shì)分析方法NDVI 及其影響因素的空間趨勢(shì)分析可采用基于柵格尺度的普通最小二乘法（OLS），通過(guò)模擬每個(gè)柵格NDVI 及其影響因素在相應(yīng)時(shí)段的變化趨勢(shì)，從而反映研究區(qū)在不同研究時(shí)段內(nèi)的NDVI 及其影響因素空間特征?？臻g趨勢(shì)分析的計(jì)算公式為：

式中：θslope為線性回歸的斜率，表示不同時(shí)間尺度下NDVI、回歸參數(shù)、殘差等研究對(duì)象在流域或像元內(nèi)的年際變化趨勢(shì)；n 為研究時(shí)段年份數(shù)；Yi為第i年研究對(duì)象數(shù)值。

2.2 NDVI 與影響因素的空間異質(zhì)性表征方法地理加權(quán)回歸（Geographically Weighted Regression，GWR）是由Brunsdon 等［26］提出的一種空間分析方法，有助于揭示一定的研究區(qū)域內(nèi)自變量與因變量之間的非平穩(wěn)空間相關(guān)性。GWR 模型是對(duì)普通線性回歸的重要拓展，該方法的參數(shù)是空間位置的函數(shù)，該模型的表達(dá)形式為：

式中： yi、 xik、 εi分別為空間上i 點(diǎn)的因變量、自變量和模型殘差；(ui，vi)為i 點(diǎn)的空間位置；m為樣本量i 點(diǎn)處自變量的個(gè)數(shù)，本文取2，分別表示降水量和氣溫； βk為i 點(diǎn)上的回歸系數(shù)； β0為截距。在本文中，年均NDVI 為因變量，降水、年平均氣溫為相應(yīng)自變量，用以分析NDVI 對(duì)氣候要素的空間響應(yīng)關(guān)系。

GWR 模型采用局部加權(quán)最小二乘法求解，其中的權(quán)重是待估點(diǎn)到各觀測(cè)點(diǎn)的空間距離的函數(shù)，具有距離衰減效應(yīng)。參數(shù)公式為：

式中：β( ui，vi)為回歸系數(shù)的無(wú)偏估計(jì)；W (ui，vi)為權(quán)重矩陣，一般采用高斯模型作為權(quán)重函數(shù)，具體參見(jiàn)文獻(xiàn)［27］；X 和Y 分別為自變量和因變量的矩陣。

本文通過(guò)建立GWR 模型得到年降水量和年平均氣溫對(duì)應(yīng)的回歸參數(shù)，并計(jì)算得到NDVI 的模擬值（NDVIcc），用來(lái)表征氣候因素對(duì)NDVI 的影響；然后計(jì)算NDVI 實(shí)際觀察值（NDVIobs）與模擬值（NDVIcc）之間的差值，表示人類(lèi)活動(dòng)對(duì)NDVI 的影響：

2.3 植被變化對(duì)徑流影響的研究方法流域水量平衡方程為：

式中：R 為徑流深，mm；P 為降水量，mm；ET 為實(shí)際蒸散發(fā)量，mm； ΔS 為流域儲(chǔ)水量的變化，mm。當(dāng)研究時(shí)間尺度較大時(shí)， ΔS ≈0，本文為了減少永定河山區(qū)流域ΔS的影響，取各變量時(shí)間尺度為10年的滑動(dòng)平均值。

在流域尺度上，降雨量可直接通過(guò)雨量站點(diǎn)觀測(cè)數(shù)據(jù)進(jìn)行空間插值獲得，實(shí)際蒸散發(fā)可基于Bu?dyko 假設(shè)的水熱耦合平衡方程計(jì)算：

式中：ET0為年均潛在蒸散發(fā)量，mm；n 為反映流域下墊面特征的參數(shù)（簡(jiǎn)稱(chēng)下墊面參數(shù)），與地形、土壤和植被等一系列因素有關(guān)。對(duì)于流域下墊面參數(shù)n 的估計(jì)，諸多學(xué)者提出了不同的估算方法，這些方法均認(rèn)為植被是影響n 的重要因素。Li 等［27］針對(duì)全球典型大流域研究分析發(fā)現(xiàn)，參數(shù)n 與NDVI存在較好的擬合關(guān)系：

式中：a、b 為回歸參數(shù)，對(duì)特定流域?yàn)槌?shù)，a 的大小表征了n 對(duì)NDVI 變化的敏感性程度，b 則反映了NDVI 以外的其它因素對(duì)n 的影響。

3 研究區(qū)域和數(shù)據(jù)

3.1 研究區(qū)域永定河流域位于海河流域，地理位置介于東經(jīng)111.95°—116°22′，北緯38°90′—41°16′，流經(jīng)內(nèi)蒙古、山西、河北、北京、天津五個(gè)省級(jí)行政區(qū)，集水面積約4.7 萬(wàn)km2，占海河水系的14.7%。2016年流域總?cè)丝诩s813 萬(wàn)，城鎮(zhèn)化率為58.3%。永定河流域上游一級(jí)支流為桑干河和洋河，兩支合并后始稱(chēng)為永定河。永定河流域?qū)侔霛駶?rùn)半干旱季風(fēng)氣候，多年平均降水量介于360 ~650 mm，6—8月降水量約占年均降水量的80%。干流官?gòu)d水庫(kù)以上流域?qū)儆谏絽^(qū)，面積為4.51 萬(wàn)km2，占全流域總面積95.8%，是主要產(chǎn)流區(qū)。冊(cè)田水庫(kù)、響水堡水庫(kù)分別為永定河一級(jí)支流桑干河、洋河的主要水庫(kù)（圖1），集水面積分別為1.71 萬(wàn)km2和1.45 萬(wàn)km2。

圖1 永定河山區(qū)氣象水文站點(diǎn)和河道斷面分布

3.2 數(shù)據(jù)來(lái)源與處理

3.2.1 NDVI 數(shù)據(jù) 本文收集了1982—2015年的永定河山區(qū)逐月NDVI 數(shù)據(jù)，數(shù)據(jù)來(lái)源于GIMMS ND?VI3g 數(shù)據(jù)集（https：//ecocast.arc.nasa.gov/data/pub/gimms/），空間分辨率為8 km，根據(jù)研究區(qū)域空間范圍，提取了相應(yīng)網(wǎng)格單元逐月NDVI 數(shù)據(jù)，經(jīng)進(jìn)一步處理后得到季尺度和年尺度的NDVI 數(shù)據(jù)。

3.2.2 氣象數(shù)據(jù) 本文收集了1982—2015年的MSWEP 逐日降水?dāng)?shù)據(jù)和中國(guó)區(qū)域地面氣象要素驅(qū)動(dòng)數(shù)據(jù)集中的逐日平均氣溫?cái)?shù)據(jù)。兩者通過(guò)網(wǎng)絡(luò)共享獲取（網(wǎng)址分別為：http：//www.gloh2o.org/和http：//www.tpdc.ac.cn/zh-hans/），空間分辨率均為0.1°×0.1°。同時(shí)，在流域內(nèi)可獲取國(guó)家氣象科學(xué)數(shù)據(jù)中心（http：//cdc.cma.gov.cn/home.do）提供的7 個(gè)國(guó)家基本氣象站點(diǎn)1982—2015年逐日常規(guī)地表氣象數(shù)據(jù)（含雨量和氣溫）。

對(duì)永定河山區(qū)0.1°×0.1°的逐日MSWEP 網(wǎng)格降水?dāng)?shù)據(jù)和中國(guó)區(qū)域地面氣象要素驅(qū)動(dòng)數(shù)據(jù)集中的平均氣溫?cái)?shù)據(jù)，采用IDW 方法，分別采用國(guó)家基本氣象站的降水與氣溫?cái)?shù)據(jù)進(jìn)行校正，之后對(duì)校正后的逐日降水、平均氣溫網(wǎng)格數(shù)據(jù)進(jìn)行重采樣處理，得到與NDVI數(shù)據(jù)空間分辨率一致的逐日降水、平均氣溫格網(wǎng)數(shù)據(jù)，最后通過(guò)時(shí)間尺度上的聚合處理，得到月、季、年等尺度的降水、平均氣溫?cái)?shù)據(jù)。

3.2.3 地表徑流數(shù)據(jù) 根據(jù)海河流域水文年鑒和各省市水資源公報(bào)，收集了冊(cè)田水庫(kù)、響水堡水庫(kù)和官?gòu)d水庫(kù)所在斷面1956—2016年的實(shí)測(cè)徑流量，采用分項(xiàng)還原法，對(duì)天然年徑流量計(jì)算有影響的水量進(jìn)行還原，將被人類(lèi)開(kāi)發(fā)利用的水量（如農(nóng)業(yè)用水、工業(yè)用水、生活用水、環(huán)境用水耗損量、跨流域引調(diào)水量等）逐項(xiàng)還原到天然年徑流量中，得到永定河典型斷面的天然徑流量。

4 結(jié)果分析與討論

4.1 NDVI 的時(shí)空變化特點(diǎn)圖2 給出了1980年和2015年研究區(qū)域土地利用類(lèi)型分布，從中可知永定河山區(qū)主要土地利用類(lèi)型是耕地、林地、草地、水域、建設(shè)用地及未利用土地，其中又以耕地、林地和草地為主。2015年較1980年林地面積增加了2364.4 km2，相對(duì)增幅34.6%；耕地、草地面積分別減少1610.1 和1056.8 km2，相對(duì)降幅分別為7.7%和7.5%；建設(shè)用地增加了926.9 km2，相對(duì)增幅116.2%。

圖2 永定河山區(qū)1980年和2015年土地利用類(lèi)型

表1 為永定河山區(qū)土地利用面積轉(zhuǎn)移矩陣，反映出1980—2015年土地利用變化的主要特征是耕地、林地和草地之間的相互轉(zhuǎn)移以及部分耕地、草地轉(zhuǎn)化為建設(shè)用地；其結(jié)果是，流域內(nèi)耕地占總面積的比重由1980年的47.3%下降至2015年的43.7%，林地和草地面積之和占總面積的比重由47.3%上升至50.2%，而建設(shè)用地占總面積的比重由1.8%上升至3.9%，不同類(lèi)型土地利用變化量占流域總面積的比重較小且變化范圍介于0.4% ~ 5.3%之間，因此，總體而言，1980—2015年期間永定河山區(qū)土地利用格局變化不大。

表1 1980—2015年永定河山區(qū)土地利用變化轉(zhuǎn)移矩陣 （單位：km2）

圖3 給出了永定河山區(qū)1982—2015年年均NDVI 變化，從3 可知，年均NDVI 呈波動(dòng)增加趨勢(shì)。在研究時(shí)段內(nèi)，永定河山區(qū)年均NDVI 的波動(dòng)范圍是0.2722 ~ 0.3394，最小值和最大值分別出現(xiàn)在1985年和2013年?？傮w上，1982—2015年年均NDVI 線性趨勢(shì)率為0.0013/a（p<0.01），這說(shuō)明表明永定河山區(qū)植被恢復(fù)較明顯。圖4 表明，1982—2015年永定河山區(qū)NDVI 的年際變化趨勢(shì)具有一定的空間異質(zhì)性；年均NDVI 呈增加和減少趨勢(shì)的區(qū)域面積分別占總面積的95.36%和4.64%。其中，年均NDVI 增長(zhǎng)較快（slope≥2×10-3/a）的面積約占研究區(qū)域總面積的12.32%，主要分布在桑干河上游南部林地和永定河官?gòu)d水庫(kù)附近的濕地；年均NDVI 呈降低趨勢(shì)的區(qū)域分布比較零散，主要在城市擴(kuò)張的地區(qū)，且降低速度率較小，多在-1×10-3/a 以內(nèi)。

圖3 1982—2015年永定河山區(qū)年均NDVI年際變化

圖4 1982—2015年永定河山區(qū)年均NDVI 趨勢(shì)空間分布

圖5 給出了永定河山區(qū)1982—2015年春、夏、秋季的NDVI 變化情況，從中可知春、夏、秋季平均NDVI 均呈波動(dòng)增加趨勢(shì)。在研究時(shí)段內(nèi)，永定河山區(qū)春季平均NDVI 的變化范圍是0.2196 ~0.2746，最小值和最大值分別出現(xiàn)在1985年和2014年，其趨勢(shì)率為0.0013/a（p<0.01）；夏季平均NDVI的變化范圍是0.4317 ~ 0.5512，最小值和最大值分別出現(xiàn)在1985年和2013年，趨勢(shì)率為0.0018/a（p<0.01）；秋季平均NDVI 的變化范圍是0.2626 ~ 0.3669，最小值和最大值分別出現(xiàn)在1984年和2013年趨勢(shì)率為0.0020/a（p<0.01）。這表明在春、夏、秋季永定河山區(qū)植被均有明顯恢復(fù)。由圖6 可知，1982—2015年永定河山區(qū)NDVI 變化趨勢(shì)在春、夏、秋季均具有一定的空間異質(zhì)性。其中，春季平均NDVI呈增加和減少趨勢(shì)的區(qū)域面積分別占總面積的92.03%和7.97%， NDVI增長(zhǎng)較快（slope≥2×10-3/a）的面積約占18.26%，主要分布在桑干河上游南部林地和永定河官?gòu)d水庫(kù)附近的濕地地區(qū)；夏季平均NDVI呈增加和減少趨勢(shì)的區(qū)域面積分別站總面積的90.58%和9.42%，NDVI增長(zhǎng)較快（slope≥2×10-3/a）的面積所占比例高達(dá)46.40%，主要分布在桑干河中上游地區(qū)林地、洋河流域北部的林草地和官?gòu)d水庫(kù)附近的濕地地區(qū)；秋季平均NDVI 呈增加和減少趨勢(shì)的區(qū)域面積分別站總面積的98.41%和1.59%，NDVI 增長(zhǎng)較快（slope≥2×10-3/a）的面積約占51.30%，主要分布情況與夏季類(lèi)似。

圖5 1982—2015年永定河山區(qū)春、夏、秋季平均NDVI年際變化率年際變化

圖6 1982—2015年永定河山區(qū)春、夏、秋季平均NDVI 變化率的空間分布

因此，在土地利用格局變化不大的情況下，永定河山區(qū)NDVI 在年度和春、夏、秋季均有較大幅度的增加，可能由于受降水、氣溫的氣象因素或流域內(nèi)植被恢復(fù)工程、農(nóng)林草地管理水平的提高等人類(lèi)活動(dòng)的影響，使年平均和各季度的NDVI 均有較大提高。

4.2 NDVI 與相關(guān)驅(qū)動(dòng)因子的空間非平穩(wěn)性

4.2.1 驅(qū)動(dòng)因子回歸參數(shù)與殘差的均值 為驗(yàn)證植被NDVI 與降水、氣溫的空間相關(guān)關(guān)系是否具有非平穩(wěn)性，首先對(duì)回歸參數(shù)均值進(jìn)行空間自相關(guān)分析，得到各要素GWR 回歸參數(shù)的Moran 指數(shù)與Z 值如表2 所示。表2 結(jié)果表明，NDVI 與降水、平均氣溫的GWR 模型回歸參數(shù)均值的Moran 指數(shù)均大于0，說(shuō)明參數(shù)具有正的空間自相關(guān)，且Z 值均大于2.58，說(shuō)明計(jì)算參數(shù)在0.01 顯著水平下具有統(tǒng)計(jì)意義。因此采用GWR 方法分析降水、平均氣溫等要素在不同地理位置對(duì)NDVI 的定量影響是合理的。通過(guò)GWR 模型殘差項(xiàng)可用來(lái)近似表征除降水、平均氣溫以外的其他因素對(duì)NDVI 的定量影響，包括地理位置、太陽(yáng)輻射、人類(lèi)活動(dòng)等因素，結(jié)合流域的特點(diǎn)和前人研究［28］，人類(lèi)活動(dòng)影響最大，本文用殘差體現(xiàn)人類(lèi)活動(dòng)的影響作用；另外，殘差作為NDVI 實(shí)測(cè)值與模擬值之差具有一定的不確定性，這種不確定性主要來(lái)源于輸入信息（如遙感數(shù)據(jù)、實(shí)測(cè)數(shù)據(jù)）的不確定性、模型結(jié)構(gòu)的不確定性等方面，對(duì)結(jié)果影響不大。

表2 各要素GWR 回歸參數(shù)的Moran 指數(shù)與Z 值

永定河山區(qū)1982—2015年年均NDVI 與降水、平均氣溫的GWR 空間回歸參數(shù)及殘差的均值分布如圖7 所示。從圖7 可知，在1982—2015年期間，總體上NDVI 與降水、平均氣溫表現(xiàn)為正相關(guān)關(guān)系，殘差也為正值，這說(shuō)明降水、平均溫度、人類(lèi)活動(dòng)均對(duì)植被活動(dòng)有一定的促進(jìn)作用。同時(shí)，回歸參數(shù)及殘差的均值表現(xiàn)出顯著的空間異質(zhì)性。NDVI 與年降水量的GWR 空間回歸參數(shù)在桑干河、洋河中上游地區(qū)介于5.9×10-4~ 7×10-4，表明這些地區(qū)年均降水量每升高100 mm 會(huì)導(dǎo)致NDVI 增加5.9×10-2~ 7×10-2；同時(shí)NDVI 對(duì)降水變化的敏感性從上游至下游逐漸加大，下游地區(qū)的GWR 空間回歸參數(shù)提升至1.2×10-3~ 1.5×10-3。NDVI 與年均氣溫的GWR 空間回歸參數(shù)介于0 ~ 7.2×10-3，其中約45%的區(qū)域低于1×10-3，主要分布在桑干河和洋河下游以及官?gòu)d水庫(kù)以下地區(qū)。這說(shuō)明流域內(nèi)約一半的范圍NDVI 對(duì)平均溫度的變化不甚敏感，溫度每升高1℃導(dǎo)致植被NDVI 的增加值低于1×10-3。GWR模型殘差值在各像元內(nèi)均為正，表明永定河山區(qū)人類(lèi)活動(dòng)對(duì)NDVI 為促進(jìn)作用，且在流域周邊地區(qū)殘差值較高，這可能由于永定河山區(qū)周邊地區(qū)主要為林地、草地，而中心區(qū)域主要為耕地，近年來(lái)人工植被恢復(fù)措施導(dǎo)致林草地NDVI 增長(zhǎng)較耕地更快。

圖7 1982—2015年永定河山區(qū)年均NDVI 與降水、平均溫度的GWR 回歸參數(shù)和殘差均值空間分布

1982—2015年永定河山區(qū)不同季度的NDVI 與各影響因素的GWR 回歸參數(shù)均值分布如圖8 所示。從圖8 可以看出，各季度GWR 回歸參數(shù)分布與年平均情況下具有相似空間分布特征，且回歸參數(shù)均為正值。

圖8 1982—2015年永定河山區(qū)各季度NDVI 與降水、平均溫度的GWR 回歸參數(shù)和殘差均值空間分布

NDVI 與降水GWR 的空間回歸參數(shù)值在春、秋季大致介于2.5×10-3~ 7.5×10-3之間，NDVI 變化對(duì)降水敏感性最強(qiáng)的區(qū)域分布在永定河下游地區(qū)；而在夏季，GWR 回歸參數(shù)值介于7.6×10-4~ 3.68×10-3，表明這一季節(jié)NDVI 對(duì)降水的敏感性程度相對(duì)較弱，洋河流域北部地區(qū)敏感性程度相對(duì)較高。春、夏、秋季NDVI 與年平均氣溫的GWR 回歸參數(shù)大于2×10-3的面積占總面積的比重分別為27.4%、51.2%和61.2%，可以看出春季NDVI受溫度影響的作用較夏、秋季明顯偏弱；在空間上，各季度NDVI對(duì)溫度變化較敏感的區(qū)域主要分布在桑干河、洋河流域上游。不同季節(jié)GWR 模型的殘差項(xiàng)均為正，其中春、夏、秋季殘差項(xiàng)大于0.05 的范圍分別占流域面積的10.4%、22.5%、13.0%，主要分布在流域周邊林草地和濕地，這說(shuō)明人類(lèi)活動(dòng)在春、夏、秋季對(duì)NDVI 均為促進(jìn)作用，且對(duì)夏季NDVI 受人類(lèi)活動(dòng)影響作用最大，秋季次之、春季最少。

4.2.2 驅(qū)動(dòng)因子回歸參數(shù)與殘差的趨勢(shì) 永定河山區(qū)1982—2015年年均NDVI 與降水、平均氣溫的GWR 回歸參數(shù)及殘差的時(shí)間變化趨勢(shì)如圖9 所示。從圖9 可知，永定河山區(qū)NDVI 與年降水的GWR空間回歸參數(shù)呈下降趨勢(shì)的面積占研究區(qū)域總面積的比例接近17%，主要分布在洋河流域，而其余83%的面積對(duì)應(yīng)的像元回歸參數(shù)均呈上升趨勢(shì)，這表明1982—2015年流域內(nèi)多數(shù)地區(qū)NDVI 受降水的影響程度呈增加趨勢(shì)；NDVI 與年平均氣溫的GWR 回歸參數(shù)呈增加和下降趨勢(shì)的像元比例為62.9%和37.1%，其中洋河流域和官?gòu)d水庫(kù)以下區(qū)域主要呈下降趨勢(shì)，而桑干河流域主要呈上升趨勢(shì)，這說(shuō)明洋河流域和官?gòu)d水庫(kù)以下流域NDVI 受氣溫的影響程度在減弱，而桑干河流域NDVI 受氣溫的影響程度增強(qiáng)；殘差呈增加和下降趨勢(shì)的像元比例分別為39.1%和60.9%，表明人類(lèi)活動(dòng)對(duì)大多數(shù)區(qū)域NDVI的促進(jìn)有所減弱。

圖9 1982—2015年永定河山區(qū)年均NDVI 與降水、平均溫度的GWR 回歸參數(shù)和殘差趨勢(shì)空間分布

1982—2015年永定河山區(qū)不同季度下NDVI 與各影響因素的GWR 空間回歸參數(shù)趨勢(shì)分布如圖10所示。降水方面，在春、夏季流域大部分地區(qū)（除東部外）NDVI 與降水的GWR 的回歸參數(shù)斜率為正，而秋季NDVI 與降水的GWR 回歸參數(shù)斜率在全流域內(nèi)為負(fù)，表明在永定河山區(qū)大部分區(qū)域春、夏季NDVI 受降水的促進(jìn)作用不斷增加，秋季NDVI 受降水的促進(jìn)作用在全流域內(nèi)呈減小趨勢(shì)；平均氣溫方面，春季流域內(nèi)大多數(shù)區(qū)域平均氣溫與NDVI 的GWR 空間回歸參數(shù)斜率為負(fù)，其中洋河流域回歸參數(shù)斜率小于-20×10-5/a，表明春季流域總體上平均溫度對(duì)NDVI 的促進(jìn)作用呈減緩趨勢(shì)，且以洋河流域降低程度最為顯著；夏季平均氣溫與NDVI 的GWR 空間回歸參數(shù)斜率為負(fù)的區(qū)域主要分布在洋河流域中下游地區(qū)，占總面積比例為31.3%，其余占流域68.7%的地區(qū)NDVI 受平均氣溫的促進(jìn)作用呈增加趨勢(shì)；秋季NDVI 受溫度的促進(jìn)作用呈減緩趨勢(shì)的面積占44.2%，分布區(qū)域較夏季有所擴(kuò)大。殘差方面，春、夏、秋季殘差趨勢(shì)呈增加的面積占流域面積的比重分別為45.9%、40.7%和45.5%，且主要分布在永定河山區(qū)周邊的林草地地區(qū)，而中部大部分區(qū)域NDVI 受人類(lèi)活動(dòng)的促進(jìn)作用呈減弱趨勢(shì)。

圖10 1982—2015年永定河山區(qū)春、夏、秋季NDVI 與降水、平均溫度的GWR 回歸參數(shù)及殘差的趨勢(shì)空間分布

4.3 NDVI 變化對(duì)流域天然徑流的影響人類(lèi)活動(dòng)對(duì)河川徑流的影響主要包括兩部分內(nèi)容，其一主要為流域內(nèi)人類(lèi)取用水和水利工程調(diào)節(jié)的影響，其二主要為下墊面變化（土地利用和覆被變化）的影響［29］。對(duì)于天然徑流，鑒于已對(duì)社會(huì)經(jīng)濟(jì)取用水和水利工程調(diào)節(jié)進(jìn)行還原，天然徑流變化可不考慮前一種影響，僅考慮下墊面變化的的影響。圖11 為永定河山區(qū)1956—2016年典型斷面（冊(cè)田水庫(kù)斷面、響水堡水庫(kù)斷面、官?gòu)d水庫(kù)斷面）天然徑流量，可以看出流域天然徑流量的年際衰減趨勢(shì)十分明顯，這是氣候變化、土地利用和覆被變化共同作用的結(jié)果。

為了進(jìn)一步探究NDVI 變化對(duì)永定河山區(qū)天然徑流的影響，首先對(duì)天然徑流變異性進(jìn)行了診斷，采用滑動(dòng)T 檢驗(yàn)法對(duì)各斷面1956—2016年天然徑流量進(jìn)行了突變?cè)\斷（采用的顯著性水平α=0.05，統(tǒng)計(jì)量臨界值為1.67），診斷結(jié)果如圖12 所示。由該圖可知，冊(cè)田水庫(kù)、響水堡水庫(kù)、官?gòu)d水庫(kù)斷面的滑動(dòng)T 檢驗(yàn)統(tǒng)計(jì)量分別在1982年、1983年、1983年取得最大值，且均顯著高于臨界值，這說(shuō)明桑干河、洋河和永定河上游年際天然徑流量在1982年、1983年、1983年發(fā)生了突變。為便于分析，本文統(tǒng)一將1983年作為3 個(gè)斷面天然年徑流量的變異點(diǎn)，同時(shí)下文將1956—1982年記作基準(zhǔn)期，將1983—2016年記作變異期。結(jié)合圖11分析，官?gòu)d水庫(kù)斷面基準(zhǔn)期年均天然徑流量43.4 mm，變異期為21.2 mm，相對(duì)減少51.1%；冊(cè)田水庫(kù)斷面基準(zhǔn)期、變異期年均天然徑流量分別為37.8 和22.5 mm，相對(duì)減少40.5%；響水堡水庫(kù)斷面基準(zhǔn)期、變異期年均天然徑流量分別為51.0 和21.9 mm，相對(duì)減少57.1%。

圖11 1956—2016年永定河山區(qū)主要斷面地表天然徑流量

圖12 1956—2016年永定河山區(qū)控制斷面天然地表徑流量滑動(dòng)T 檢驗(yàn)統(tǒng)計(jì)量

圖13 給出了永定河山區(qū)干支流斷面在不同階段的年降水量-天然徑流量散點(diǎn)圖及線性回歸關(guān)系。從圖13 可知，基準(zhǔn)期內(nèi)降水量與天然徑流量有較好相關(guān)性，冊(cè)田水庫(kù)、響水堡水庫(kù)和官?gòu)d水庫(kù)斷面降水-天然徑流量線性回歸方程的決定系數(shù)分別為0.57、0.63 和0.59，表明基準(zhǔn)期降雨量對(duì)天然徑流有較強(qiáng)的解釋能力；而變異期內(nèi)降水-徑流量線性回歸方程的決定系數(shù)變?yōu)?.17、0.07 和0.10，降水量對(duì)天然徑流的解釋能力明顯變差。變異期年降水-天然徑流量散點(diǎn)基本位于基準(zhǔn)期散點(diǎn)下方，同等年降水量條件下永定河山區(qū)產(chǎn)流能力大幅度下降。根據(jù)兩個(gè)時(shí)期的年降水-天然徑流回歸關(guān)系，計(jì)算得到變異期內(nèi)由降水減少導(dǎo)致各斷面天然徑流累計(jì)縮減量分別為3.64、4.26 和3.34 mm，對(duì)天然徑流減少的貢獻(xiàn)率分別為16.4%、14.6%和21.8%。相關(guān)研究［30］表明，海河流域潛在蒸散發(fā)變化對(duì)徑流的變化影響微弱，降水減少和由于人類(lèi)活動(dòng)導(dǎo)致的下墊面變化是流域徑流量減少的主導(dǎo)因素。由于降水減少對(duì)永定河徑流減少的貢獻(xiàn)率僅為20%左右，因此導(dǎo)致永定河流域天然徑流減少的主要原因是土地利用和覆被等下墊面條件的變化，由前述分析可知，永定河山區(qū)土地利用格局總體變化不大，而NDVI 呈顯著變大的趨勢(shì)，這進(jìn)一步表明永定河山區(qū)天然徑流衰減是由NDVI 的增加所引起的。

圖13 永定河山區(qū)干支流斷面年降雨量-天然徑流量關(guān)系

通過(guò)式（5）與式（6）得到冊(cè)田水庫(kù)、響水堡水庫(kù)、官?gòu)d水庫(kù)控制流域內(nèi)在不同年份的下墊面參數(shù)n，發(fā)現(xiàn)永定河山區(qū)三個(gè)水庫(kù)控制流域n 的多年平均值分別為3.18、2.72、2.74，均大于傅抱璞模型推薦的n 的缺省值［31］，因此永定河山區(qū)的植被水分利用效率要高于平均水平，這是由于永定河山區(qū)處于中國(guó)半干旱地區(qū)，實(shí)際蒸散發(fā)量占降水量的比重高于90%，屬于典型的水分限制區(qū)，因此參數(shù)n 的均值大于平均值符合水熱平衡規(guī)律，且在各控制流域中，冊(cè)田水庫(kù)控制流域的持水能力最大，而響水堡水庫(kù)控制流域的持水能力最小。同時(shí)發(fā)現(xiàn)，與天然徑流的變化趨勢(shì)相反，永定河山區(qū)下墊面參數(shù)n和年均NDVI 隨時(shí)間呈增加趨勢(shì)，其中冊(cè)田水庫(kù)、響水堡水庫(kù)、官?gòu)d水庫(kù)下墊面參數(shù)n 的年均增長(zhǎng)率分別為0.0359、0.0426、0.0433，對(duì)應(yīng)各流域年均NDVI 的增長(zhǎng)率分別為0.0011、0.0013、0.0013，其中，響水堡水庫(kù)與官?gòu)d水庫(kù)的下墊面參數(shù)n 和NDVI 增長(zhǎng)情況表現(xiàn)出高度一致性，均大于冊(cè)田水庫(kù)控制流域。繪制了不同年份滑動(dòng)平均的下墊面參數(shù)n 與NDVI 點(diǎn)據(jù)分布，對(duì)下墊面參數(shù)n 與NDVI 進(jìn)行擬合，如圖14 和表3 所示，可以看出下墊面參數(shù)n 與NDVI 有較好的正相關(guān)關(guān)系，線性回歸系數(shù)a 分別為17.90、22.14、20.97，截距b 分別為-2.19、-2.84、-3.45，表明NDVI 每增加1%時(shí)，下墊面參數(shù)n將增加20%左右，這意味著永定河流域產(chǎn)流性能對(duì)植被變化有極強(qiáng)的敏感性?；貧w方程的決定系數(shù)分別為0.44、0.45、0.46，對(duì)應(yīng)相關(guān)系數(shù)擬合效果均達(dá)到顯著性水平（P<0.01）。

表3 永定河山區(qū)Budyko 模型中下墊面參數(shù)n 與植被NDVI 擬合結(jié)果

圖14 永定河山區(qū)典型流域內(nèi)年均NDVI 與下墊面參數(shù)的關(guān)系

永定河山區(qū)NDVI 的增加極大促進(jìn)了流域河川天然徑流量的衰減作用，結(jié)合前述NDVI 與相關(guān)驅(qū)動(dòng)因素的GWR 模型計(jì)算結(jié)果，NDVI 增加與降水量、平均氣溫作用下的自然生長(zhǎng)有關(guān)，也與植被恢復(fù)相關(guān)的退耕還林、風(fēng)沙源治理等大規(guī)模人類(lèi)活動(dòng)有關(guān)。根據(jù)調(diào)查，自21 世紀(jì)以來(lái)永定河流域內(nèi)開(kāi)展實(shí)施了多項(xiàng)植被恢復(fù)工程，主要有2000年開(kāi)始實(shí)施的京津風(fēng)沙源治理一期工程、2002年實(shí)施的退耕還林工程以及2013年實(shí)施的京津唐風(fēng)沙源治理二期工程等，而2000年以前流域內(nèi)無(wú)相關(guān)大型植被恢復(fù)工程，在一定程度上說(shuō)明，2000年前植被的生長(zhǎng)主要受氣候條件的影響，而2000年后植被生長(zhǎng)同時(shí)受氣候條件和人類(lèi)活動(dòng)的共同影響作用。結(jié)合圖3 永定河山區(qū)年均NDVI 變化圖可知，1982—2000年期間年均NDVI 增加值為0.012，而2000 至2015年年均NDVI 增加0.043，后者遠(yuǎn)大于前者，表明2000年后大型植被恢復(fù)工程極大促進(jìn)了永定河山區(qū)NDVI 的增加，同時(shí)也引起流域生態(tài)系統(tǒng)蒸散作用增強(qiáng)，導(dǎo)致永定河山區(qū)在同等降水條件下產(chǎn)流能力顯著降低。有研究已指出，大型植被措施對(duì)濕潤(rùn)地區(qū)的水資源影響較小，而在干旱半干旱地區(qū)往往會(huì)顯著減少水資源［32］，顯然永定河山區(qū)植被與河川天然徑流的變化情況與后者是相符的。1982—2015年永定河山區(qū)冊(cè)田水庫(kù)、響水堡水庫(kù)、官?gòu)d水庫(kù)控制流域的平均濕潤(rùn)指數(shù)（即降水量與潛在蒸散發(fā)的比值）分別為0.39、0.40、0.41，相應(yīng)下墊面參數(shù)n 的平均值分別為3.18、2.71、2.74，表明永定河山區(qū)屬于持水能力高的半干旱地區(qū)，結(jié)合文獻(xiàn)［33］給出的中國(guó)人工植被恢復(fù)的原則，建議未來(lái)永定河流域植被恢復(fù)應(yīng)以自然修復(fù)為主，對(duì)人工植被恢復(fù)工程應(yīng)進(jìn)行適當(dāng)?shù)目刂?，同時(shí)遵守植被分布規(guī)律來(lái)界定植被分區(qū)［34］。

5 結(jié)論

通過(guò)構(gòu)建地理加權(quán)回歸模型探究了1982—2015年永定河山區(qū)NDVI 和降水、平均氣溫的空間非平穩(wěn)關(guān)系，并進(jìn)一步剖析了NDVI 變化對(duì)流域典型斷面天然徑流的影響，結(jié)果如下：（1）1982—2015年永定河山區(qū)土地利用格局變化不大，但流域尺度上植被恢復(fù)明顯，年均和春、夏、秋季NDVI平均趨勢(shì)率為0.0013/a、0.0013/a、0.0018/a、0.0020/a，且流域內(nèi)絕大多數(shù)空間像元NDVI 的年際增加趨勢(shì)均顯著。（2）1982—2015年，永定河山區(qū)年、季尺度下的NDVI 和降水、平均氣溫的回歸參數(shù)、殘差均值以及兩者的年際變化率均具有明顯空間異質(zhì)性。其中回歸參數(shù)及殘差均值在空間像元內(nèi)均為正，表明流域內(nèi)降水、平均氣溫等氣候因素和人類(lèi)活動(dòng)因素對(duì)NDVI均表現(xiàn)為增加作用，但在不同空間位置的差異亦較為顯著。（3）1956—2016年永定河山區(qū)天然徑流量具有顯著的衰減趨勢(shì)，結(jié)合水熱耦合平衡方程發(fā)現(xiàn)流域天然徑流量對(duì)植被變化具有高度敏感性，以植被變化為突出特征的下墊面變化對(duì)永定河天然徑流衰減具有主導(dǎo)作用。（4）永定河山區(qū)屬于持水能力高的半濕潤(rùn)半干旱地區(qū)，流域產(chǎn)流能力對(duì)植被覆蓋狀況的變化較為敏感。大規(guī)模的恢復(fù)工程的建設(shè)對(duì)河川天然徑流衰減的可能影響應(yīng)給予高度重視，建議未來(lái)永定河流域植被恢復(fù)應(yīng)以自然修復(fù)為主，對(duì)人工植被恢復(fù)工程應(yīng)進(jìn)行適當(dāng)控制。