侯 帥，張吳平，賈若男

(山西農(nóng)業(yè)大學(xué)資源環(huán)境學(xué)院，山西 太谷 030801)

0 引 言

長期以來，煤炭資源都是中國工業(yè)發(fā)展的重要支撐力量，特別是20世紀(jì)80年代以來，煤炭資源為我國的經(jīng)濟建設(shè)和發(fā)展做出了巨大的貢獻，但長期的煤炭開采造成水循環(huán)途徑變化、土地和植被損毀、生態(tài)環(huán)境退化[1]。煤炭開采改變了地表形態(tài)，破壞了地表原有的水文過程，研究煤炭開采對地表水文過程的影響，可為煤礦區(qū)的環(huán)境保護、生態(tài)重建、科學(xué)開采等提供重要的理論依據(jù)。

自20世紀(jì)50年代以來，水文模型的研究取得了巨大的進展，其中分布式水文模型可更好地分析自然和人為因素對流域水環(huán)境的影響[2]，通常水文模型對水文過程的模擬主要包括植被截留、土壤入滲、徑流匯流等過程，需要大量的水文、氣象觀測資料，而無觀測站且受采煤擾動的長河流域難以滿足這些要求。

植被截留降雨是水文循環(huán)中極其重要的一部分，植被截留受到降雨強度、降雨歷時、植被類型、環(huán)境因素等諸多因素的影響，是一個復(fù)雜的混合過程[3]。國外的Horton、Rutter和Gash模型的應(yīng)用較多，國內(nèi)研究和應(yīng)用較多的多為半經(jīng)驗性理論模型，如李崇巍、李玉霞等以儀垂祥模型為基礎(chǔ)，建立了岷江上游林冠截持降雨遙感模型，結(jié)果較為理想[4-8]。但對農(nóng)作物冠層截留的研究相對于林冠截留的研究較少且其結(jié)果差異較大[9]。SWAP模型中的作物截留計算采用了由Von Hoyningen-Hune和Braden提出的農(nóng)田作物冠層截留量的一般方程，并且也被應(yīng)用到其他研究水分運動的模型中去，如張圣微采用該方法建立了農(nóng)田水分運動模型的降雨截留子模塊[10]。大多數(shù)水文模型中的植被截留過程模擬以經(jīng)驗?zāi)Ｐ蜑橹?，且沒有區(qū)分天然植被和種植植被，對此過程的研究較為簡單。

土壤入滲指降雨或灌溉的水分從地表滲入土壤的運動過程，土壤入滲在水文循環(huán)過程中起著重要的紐帶作用。土壤入滲的研究模型主要有經(jīng)驗?zāi)Ｐ?，通過土壤入滲曲線得到模型參數(shù)如Kostiakov模型；半經(jīng)驗?zāi)Ｐ?，采用簡化的連續(xù)方程如Hortan、Holtan等模型；物理模型，基于Darcy定律和Richard方程，具有明確的物理意義，如Green-Ampt、Philip等模型,其中Green-Ampt模型廣泛應(yīng)用于水土侵蝕、降雨入滲等領(lǐng)域的研究中，如Mein 等將其改進應(yīng)用于恒定降雨條件下的土壤入滲研究，也被應(yīng)用到SWAT模型中計算下滲量[11,12]。

降雨經(jīng)植被截留、土壤入滲后，在地形作用下形成地表徑流。而數(shù)字高程模型(DEM)可反映地形特征，借助一定的算法可獲取水系特征，為徑流模擬提供了基礎(chǔ)數(shù)據(jù)，因此已被應(yīng)用到SWAT、SHE、TOPMODEL等諸多水文模型中[13]。

筆者以30 m×30 m的柵格為計算單元，基于3S技術(shù)構(gòu)建了長河流域地表水文過程模型，分析了煤炭開采下模型參數(shù)的空間分布特征，并模擬了不同降雨強度下的地表水文過程。

1 材料與方法

1.1 研究區(qū)概況

山西省長河流域位于晉城市澤州縣(圖1)，地理坐標(biāo)北緯35°30′～35°38′，東經(jīng)112°37′～112°46′，包括47個村莊，面積約為113 km2，海拔725～1 167 m，是典型的黃土丘陵區(qū)。氣候?qū)倥瘻貛Т箨懶约撅L(fēng)氣候，年平均氣溫11.2 ℃，年平均降雨量576 mm，長河于中部自北向南流過，長河以西分布著年產(chǎn)均在45 萬t以上的12座煤礦，長期的煤炭開采造成土地破壞、水資源污染、空氣污染等一系列問題，已在很大程度上影響了當(dāng)?shù)厝嗣竦纳a(chǎn)生活和社會經(jīng)濟的健康發(fā)展。

圖1 研究區(qū)地理位置及高程Fig.1 Geographical location and elevation of the study area

1.2 地表水文過程模型

研究區(qū)地處黃土高原，進行短時間徑流模擬可以不考慮蒸散發(fā)和壤中流的影響，因此本文重點模擬地表水文過程中的植被截留降雨、土壤入滲、徑流匯流過程。地表水文過程模型建立的條件為：以分水嶺劃定研究區(qū)邊界，研究區(qū)內(nèi)居民以井水為主要水源，在情景模擬中無外界河流匯入研究區(qū)。

本文構(gòu)建的地表水文過程模型包括植被截留模型、土壤入滲模型、徑流匯流模型，模型的總體結(jié)構(gòu)及各模型間的聯(lián)系如圖2所示。由給定的降雨類型，首先根據(jù)各類型植被的截留降雨模型可得到截留的降雨量，落到地表的降雨再由土壤入滲模型得到入滲量，最后地表凈雨量經(jīng)地表匯流過程模擬得到匯流累積量的空間分布。本文首次在長河流域構(gòu)建地表水文過程模型，由于研究區(qū)內(nèi)缺乏水文觀測資料，模型具體應(yīng)用于下述降雨情景模擬中。

圖2 地表水文過程模型結(jié)構(gòu)示意圖Fig.2 Schematic of the surface hydrological process model structure

1.2.1 植被截留模型

(1)天然植被截留模型。本部分采用儀垂祥提出的植被截留模型，對于給定的植被類型，植被冠層存在一個最大截留量，當(dāng)降雨量P(mm)大于最大截留量時，最大截留降雨量PJ(mm)為[7]：

PJ=αVFCLAI

(1)

式中：α為葉表面最大持水深度。mm；VFC為植被覆蓋度；LAI為葉面積指數(shù)，m2/m2。

當(dāng)降雨量小于最大截留量時，植被截留量由植被覆蓋度決定，其截留量PJ為：

PJ=VFCP

(2)

(2)種植植被截留模型。根據(jù)SWAP模型中計算作物截留的一般方程來模擬作物的截留量，隨著降雨量的增加，其截留量也逐漸達到最大值aLAI[14]，截留量可由下式計算：

(3)

式中：PJ為降雨截留量，mm；a為經(jīng)驗系數(shù)；b為植被覆蓋度。

1.2.2 土壤入滲模型

采用以下由Mein和Larson將Green-Ampt模型推廣至降雨入滲情況的模型來對入滲過程進行模擬[15]。

假定有穩(wěn)定雨強p，在降雨的開始階段，降雨全部滲入土壤，當(dāng)p大于土壤的入滲能力時，地表開始形成積水，由模型可得出積水時的累積入滲量Fp(cm)：

Fp=Sf(θs-θi)/(p/Ks-1)

(4)

由此可得開始積水時間tp=Fp/p，只有當(dāng)p>Ks時才可能積水，所以整個過程的入滲率f(cm/d)為：

(5)

式中：Ks為飽和導(dǎo)水率，cm/d；Sf為濕潤峰處的土壤水吸力值，cm；θs、θi分別為土壤的飽和含水率和初始含水率，cm3/cm3；F為開始積水以后的累積入滲量，cm。F的計算采用下面修正后的公式：

t>tp

(6)

式中：ts表示假設(shè)從t=0時開始積水，到入滲量為Fp時所需的時間，計算方法如下：

(7)

在實際應(yīng)用中，首先根據(jù)降雨強度和降雨時間判斷地表是否可以形成積水，積水前的累積入滲量由式(4)得到，積水后的累積入滲量由式(6)、(7)得出。

1.2.3 徑流匯流模型

建立基于DEM的徑流匯流模型，DEM填洼處理后，采用單流向算法中應(yīng)用較多的D8法確定水流方向，根據(jù)合適的閾值提取出河網(wǎng)，應(yīng)用Strahler法對河網(wǎng)分級，并將整個流域劃分成若干子流域，從而建立研究區(qū)的數(shù)字水系[16]。

以植被截留、土壤入滲過程得到的地表凈雨量為賦值柵格，對生成的數(shù)字水系賦值可得到地表匯流量的空間分布及各子流域的匯流累積量。

1.3 數(shù)據(jù)來源與處理

1.3.1 植被截留模型參數(shù)

(1)遙感數(shù)據(jù)。來源于美國地質(zhì)勘探局(http:∥www.usgs.gov/)的2000-2015年的Landsat影像，首先需要對原始影像進行預(yù)處理，主要是輻射校正和幾何校正。

(2)植被類型。采用研究區(qū)遙感影像結(jié)合當(dāng)?shù)赝恋乩靡?guī)劃圖及《澤州縣林地保護利用規(guī)劃》對遙感影像進行解譯，并利用研究區(qū)1∶2 000正射影像對地物判別和修正得到。

(3)葉表面最大持水深度。根據(jù)各鄉(xiāng)鎮(zhèn)的統(tǒng)計資料及文獻[17]，并結(jié)合實地調(diào)查確定。

(4)植被覆蓋度。利用獲取的遙感影像來提取歸一化植被指數(shù)(NDVI)來估算[18]。

(5)葉面積指數(shù)。通過遙感影像結(jié)合植被冠層多次散射模型來獲取，具體步驟見文獻[19]。

1.3.2 土壤入滲模型參數(shù)

模型的關(guān)鍵參數(shù)土壤飽和導(dǎo)水率由土壤理化性質(zhì)通過van Genuchten模型得到，濕潤峰處的土壤水吸力值根據(jù)張光輝等的研究得出[20,21]。所需的土壤的各項理化性質(zhì)通過GPS確定1 km×1 km的采樣網(wǎng)格對土壤采樣后測定，共得到117組采樣點數(shù)據(jù)。

1.3.3 徑流匯流模型參數(shù)

DEM數(shù)據(jù)來源于澤州縣國土資源局，分辨率為10 m，在ArcGIS中經(jīng)填洼、確定流向、提取河網(wǎng)、劃分子流域后可生成研究區(qū)的數(shù)字水系，其中含184個子流域，各徑流于研究區(qū)中部匯流，總體從東北流向西南，經(jīng)流域西南出口流出。

2 結(jié)果與討論

2.1 參數(shù)的空間分布特征

2.1.1 植被覆蓋度的空間分布特征

植被覆蓋度是影響植被截留的重要因素，也反映了煤炭開采對植被的影響程度。2015年研究區(qū)的植被類型的空間分布如圖3，天然植被類型主要為暖溫帶落葉闊葉林和針闊混交林，農(nóng)作物主要是玉米、小麥。草地、農(nóng)作物的分布基本相等，灌木林地、闊葉林、針闊混交林在礦區(qū)的面積遠大于非礦區(qū)。

圖3 研究區(qū)植被類型分布圖Fig.3 Distribution map of vegetation types in the study area

將2000-2015年的植被覆蓋度分為5級：①無覆蓋，VFC≤0.1；②低覆蓋，0.1

圖4 2000-2015年研究區(qū)植被覆蓋度等級的空間分布圖Fig.4 Spatial distribution of vegetation coverage grades in the study area during 2000-2015

圖5 2000-2015年礦區(qū)各植被覆蓋度等級的面積Fig. 5 Area of the vegetation coverage grades in mining area during 2000-2015

表1 2000-2015年礦區(qū)植被覆蓋度等級轉(zhuǎn)移概率矩陣Tab.1 Transfer probability matrix of vegetation coveragegrades in mining area during 2000-2015

2.1.2 土壤飽和導(dǎo)水率的空間分布特征

土壤飽和導(dǎo)水率反映了土壤的滲透能力，是土壤水分運動的重要特征之一。表2是采樣點土壤飽和導(dǎo)水率Ks的基本統(tǒng)計特征，礦區(qū)Ks的最小值、最大值、均值均小于非礦區(qū)，且Ks屬中等變異，變異性較強，但礦區(qū)變異系數(shù)大于非礦區(qū)，說明受采煤影響，礦區(qū)Ks的變化情況大于非礦區(qū)。經(jīng)Kolmogorov-Smirnov法檢驗，表明土壤飽和導(dǎo)水率符合正態(tài)分布(P=0.121>0.05)，通過地統(tǒng)計學(xué)方法對數(shù)據(jù)進行半變異函數(shù)的擬合，并通過交叉驗證檢驗其合理性和擬合精度，選取決定系數(shù)較大(0.568)、殘差較小(1.108×10-3)的最優(yōu)擬合函數(shù)模型為指數(shù)模型，根據(jù)Kriging插值方法可得到土壤飽和導(dǎo)水率的空間分布(圖6)[23]。研究區(qū)土壤飽和導(dǎo)水率的范圍為4.337～139.321 cm/d，西部礦區(qū)的平均土壤飽和導(dǎo)水率34.375 cm/d小于東部非礦區(qū)的44.426 cm/d，土壤飽和導(dǎo)水率總體上呈現(xiàn)出西部小于東部的趨勢。說明煤炭開采改變了原有的土壤性質(zhì)，土壤飽和導(dǎo)水率變小。

表2 礦區(qū)、非礦區(qū)土壤飽和導(dǎo)水率的基本統(tǒng)計特征Tab.2 Basic statistics of soil saturated hydraulic conductivityin mining area and non-mining area

圖6 研究區(qū)土壤飽和導(dǎo)水率的空間分布圖Fig.6 Spatial distribution of soil saturated hydraulic conductivity in the study area

2.2 基于多年降雨資料的降雨情景模擬

2.2.1 長河流域降雨趨勢分析

分析來源于澤州縣氣象中心的流域1986-2015年間的降雨資料可知，降雨主要集中在7、8月，兩月降雨量占年總降雨量比例的均值為41.524%。按照中國氣象局降雨量等級的劃分標(biāo)準(zhǔn)，降雨類型主要為小雨(0.1～9.9 mm)和中雨(10～24.9 mm)，分別占67.982%、20.175%，其余大雨(25～49.9 mm)占9.357%，暴雨(50～99.9 mm)、大暴雨(100～249.9 mm)共占2.486%。

2.2.2 情景模擬分析

在上述降雨資料分析的基礎(chǔ)上，以降雨強度分別為2.5、8、16 mm/h，降雨時間為1 h的降雨情景對研究區(qū)2015年7月的地表水文過程進行空間模擬。對于有植被覆蓋的柵格，首先根據(jù)植被截留模型判斷3種降雨情景下的降雨量是否達到各類型植被的最大截留降雨量：截留量未達最大截留量的柵格，降雨將不會落到地表；超過最大截留量的柵格根據(jù)土壤入滲模型中地表積水形成的時間，判斷落到地表的降雨在1 h內(nèi)是否會形成地表積水，形成積水的柵格再由徑流匯流模型進行匯流累積量的計算。對于無植被覆蓋的柵格，降雨直接落到地表，同理根據(jù)土壤入滲模型判斷降雨在1 h內(nèi)是否形成地表積水，是否進行徑流匯流的計算。計算結(jié)果表明：

(1)2.5 mm/h降雨強度下，降雨大部分被植被截留，較少降雨滲入土壤，如圖7(a)、(b)植被冠層平均截留量為0.160 mm，土壤平均入滲量分別為1.955 mm；

(2)8 mm/h降雨強度下，1 h內(nèi)植被冠層截留量達到最大，其中農(nóng)作物平均截留量為0.194 mm；闊葉林平均截留量為0.656 mm；針闊混交林平均截留量為0.442 mm；灌木林平均截留量為0.331 mm；草地平均截留量為0.197 mm。由于降雨強度小于土壤飽和入滲率，其余降雨經(jīng)土壤入滲后難以形成地表徑流，如圖7(c)，(d)植被冠層平均截留量為0.198 mm，土壤平均入滲量分別為6.672 mm；

(3)16 mm/h降雨強度下，1 h內(nèi)植被冠層截留量達到最大后，其余降雨一部分滲入土壤，一部分形成積水，植被截留量和土壤入滲量的空間分布如圖7(e)、(f)，土壤平均入滲量為13.457 mm。積水繼而形成徑流，經(jīng)徑流匯流模擬可得到研究區(qū)地表徑流匯流量的空間分布(圖8)，匯流量由兩側(cè)向中部、由北向南逐漸增加，在流域西南即第184號子流域出口達到最大值。

圖7 不同降雨強度下研究區(qū)植被截留量、土壤入滲量的空間分布圖Fig.7 Spatial distribution of vegetation interception and soil infiltration in the study area under different rainfall intensities

圖8 研究區(qū)降雨匯流量的空間分布圖Fig.8 Spatial distribution of rainfall accumulation in the study area

3 結(jié) 語

(1)充分運用3S技術(shù)，通過構(gòu)建植被截留模型、土壤入滲模型、徑流匯流模型，對長河流域的地表水文過程進行了空間模擬，模型具有明確的物理意義，其中充分考慮了不同植被類型對降雨的截留效應(yīng)。

(2)2000-2015年礦區(qū)植被覆蓋度有所增加，可見隨著退耕還林、土地整治等措施的實施，植被有所恢復(fù)，植被的增加提高了植被的截留作用，減弱了地表土壤的侵蝕，提升了森林涵養(yǎng)水源的能力，引起一系列積極的生態(tài)效應(yīng)。

(3)煤炭開采對土壤的破壞短時間內(nèi)難以恢復(fù)，礦區(qū)的土壤飽和導(dǎo)水率小于非礦區(qū)，土壤滲透能力的改變降低了降雨和灌溉用水的利用效率，地表水流失加劇，生態(tài)環(huán)境和農(nóng)業(yè)生產(chǎn)受到一定影響。

(4)長河流域在2.5 mm/h降雨強度下，降雨被植被截留；在8 mm/h降雨強度下，植被截留量達到最大后，其余降雨被土壤入滲，難以形成地表徑流；在16 mm/h降雨強度下，降雨會在短時間內(nèi)形成積水并進行匯流，匯流量于流域西南出口達到最大值。

(5)長河流域礦區(qū)的生態(tài)治理工作需進一步加強，本研究可為流域的科學(xué)管理與生態(tài)恢復(fù)工作提供科學(xué)的參考依據(jù)。