跨流域调水工程区域内的生态环境变化，本质上是生物、土壤、气候、地形等生态因子的演变^[1]。因此，基于对水生态环境质量及生态因子演变的研究，有利于工程对生态影响的分析，更有助于环境问题的对症治理。遥感技术作为调水工程水生态环境研究中的重要手段，具有多空间性^[2]、多时间性^[3]、高效性^[1]。不仅能对水生态环境进行有效的实时监测，还便于展开长期动态工作，快速、准确地反映大面积水环境的生态质量现状及发展趋势。目前，国内外不少学者对遥感技术在调水工程中的应用作了深入的研究。包洪福^[4]以南水北调中线工程丹江口水源地为研究区，基于遥感技术对区域内水生态环境中的动植物变化规律进行了研究，探究其生物多样性。Yang等^[5]使用无人机技术结合遥感数据对南水北调工程地下水位上升引起的生物量改变和缓冲带宽度进行了安全评估。Gao等^[6]结合了geographic information system (GIS)和remote sensing (RS)技术，以生态系统服务价值法为手段，对江苏省东部干线水源地的生态系统服务价值进行了计算并对该地的生态补偿情况展开了进一步研究。朱长明等^[7]综合地表水、地下水等数据和多源遥感长时间序列数据，对塔里木河下游生态输水前后区域水环境变化和生态响应进行监测和分析。

引黄济青调水工程通水以来，有效缓解了青岛及沿线城市水资源供需不均的矛盾^[8]。但在其建设和运行过程中，也对沿线的水生态环境产生了一定影响^[9]。本文以工程沿线为研究区，基于遥感技术对土地利用/覆盖、植被指数及覆盖度、地表温度、土壤干湿度等生态因子的时空变化规律展开讨论，进一步探究远距离调水工程沿岸生态环境的演变规律，为沿线生态的保护和治理提供了数据支撑和科学依据。

1 研究区与数据 1.1 研究区概况

引黄济青工程是中国“七五”期间重要的水利工程，位于118°07′—120°12′E，36°22′—37°16′N的山东半岛^[10]，由山东省滨州市打渔张闸引黄河水到青岛市崂山区白沙河净水厂。该地区地势平坦，受海洋影响，属暖温带季风气候区，具有大陆性气候的明显特征，年平均气温12.0~12.6 ℃，年平均降雨量在600~770 mm，常年平均蒸发量在1 000 mm以上，且起伏变化较大。地带性土壤为典型棕色森林土，分布在缓坡地和排水良好的平地。天然植被类型为暖温带落叶阔叶林，主要树种有栎类，以麻栎居多。

调水工程线路全长291.14 km，其中棘洪滩水库以上为水源工程，长253.14 km，以明渠为主，水库以下为供水工程，长38 km，以管道为主。本文以水库以上水源工程输水段为研究对象，综合考虑研究结果的代表性、准确性，建立了干渠向外左右各2 km的缓冲区，以此展开讨论。

1.2 数据来源

研究区范围广阔、地形复杂，适合卫星遥感进行大面积、多角度、高时效的生态环境监测。考虑到研究周期和影像精度的综合要求，最终选用Landsat系列影像作为初始数据。

工程通水至今已有30 a余，为研究其沿程生态因子的历史演变规律，选择以4个时期作为时间研究对象，分别是1990，2000，2010，2020年。此外，考虑到春夏两季植被生长旺盛，识别度高，而夏季云层覆盖度广，可见度低等多种因素，最终选择了5, 6月日期较接近的八景影像进行数据的处理和分析，具体影像资料见表 1。

为解决光学遥感数据获取时的误差问题，需先对下载好的影像进行预处理，包括了辐射定标和大气校正过程。其中辐射定标是为了消除传感器本身误差，大气校正是为了消除大气散射、吸收、反射引起的误差^[11]。校正后的影像再根据研究区的范围进行裁剪和镶嵌处理。

2 研究方法 2.1 技术路线

将研究区按地级市范围划分为4个区段，分别为滨州市、东营市、潍坊市及青岛市区段，基于生态遥感监测技术，选择了土地利用/覆盖、植被指数及覆盖度、地表温度、土壤干湿度等生态因子，探究其时空演变规律，讨论工程运行后沿线生态环境的变化。

2.2 土地利用/覆盖的量化

土地利用/覆盖中绿地、水域、建筑用地、裸地面积的提取，是在遥感影像经过预处理的基础上，利用ENVI软件中监督分类的方式进行划分的^[12]，结合影像和地图对分类样本进行多次调整直至样本间具有较好的可分离性。分离器的选择也是同时对比了多种分离试验结果，最终选择了支持向量机法进行分类。分类结果再以高分辨率的影像作为参照，对混淆部分进行人工目视解译，解决同物异谱或异物同谱带来的误差问题^[13]，最后利用ENVI软件中的混淆矩阵精度评价法进行验证，并对比“地理国情监测云平台”(http://www.dsac.cn/)相关数据产品，分类结果的精度达到80%以上。

2.3 植被的量化

植被的量化采用的是NDVI指数法^[14]，通常反映植被覆盖、生长等信息，其计算公式为：

$ \mathrm{NDVI}=\frac{\mathrm{NIR}-{R}}{\mathrm{NIR}+R} $

(1)

式中：NIR为近红外波段的反射率；R为红波段的反射率。

在现有NDVI值的基础上，进一步分析植被的覆盖情况，这里采用了像元二分法^[15]计算得到植被覆盖度P_V^[16]，其具体计算模型如下：

$ P_{V}=\frac{\mathrm{NDVI}-\mathrm{NDVI}_{\text {Soil }}}{\mathrm{NDVI}_{\mathrm{Veg}}-\mathrm{NDVI}_{\text {Soil }}} $

(2)

式中：P_V为植被覆盖度；NDVI为植被指数；NDVI_Soil为完全是裸土或无植被覆盖区域的NDVI值；NDVI_Veg为完全被植被所覆盖的像元的NDVI值。

2.4 地表温度的量化

本文采用大气校正法^[17](辐射传输方程)进行地表温度的反演，其本质是通过消除大气对地表热辐射的影响误差，将得到的地表热辐射强度转化为相应的地表温度。卫星传感器接收到的热红外辐射亮度值L_λ由3部分组成：①大气向上辐射亮度；②地面的真实辐射亮度经过大气层后到达卫星传感器的能量；③大气向下辐射到达地面反射的能量，其计算公式如下：

$ L_{\lambda}=\left[\varepsilon B\left(T_{s}\right)+(1-\varepsilon) L \downarrow\right] \tau+L \uparrow $

(3)

式中：ε为地表比辐射率；Ts为地表真实温度(K)；B(T_s)为黑体热辐射亮度；τ为大气在热红外波段的透过率；L↓为大气下行辐射亮度；L↑为大气上行辐射亮度。则温度为T的黑体的辐射亮度B(Ts)可表示为：

$ B(T s)=\frac{L_{\lambda}-L \uparrow-\tau(1-\varepsilon) L \downarrow}{\tau \varepsilon} $

(4)

式中：τ，L↑，L↓参数可在NASA提供的网站上(http://atmcorr.gsfc.nasa.gov/)通过提供时间、经纬度等数据信息进行获取。而地表比辐射率ε则采用NDVI阈值法计算，公式如下：

$ \varepsilon=0.004 P_{V}+0.986 $

(5)

获得所需数据参数后，利用普朗克公式计算得到地表真实温度Ts，具体公式如下：

$ T s=\frac{K_{2}}{\ln 〔\frac{K_{1}}{B(T s)}+1〕} $

(6)

式中：K₁，K₂为定标系数，可通过影像的元数据获取。

结合沿线滨州、东营、潍坊、青岛市气象站点收集到的实测均温数据，对地表温度的反演结果进行验证，结果显示两者波动变化一致，说明方法具有可行性。

2.5 土壤干湿度的量化

土壤干湿度的量化方法选用的是NDVI即温度植被干旱指数法^[18]，其原理是在光学与热红外遥感通道数据的基础上反演植被覆盖土层表面的水分情况，本质是结合归一化植被指数(NDVI)和地表温度(T_s)两个参数，构建T_s-NDVI特征空间，其公式表示为：

$ \mathrm{TVDI}=\frac{T_{s}-T_{s \min }}{T_{s \max }-T_{s \min }} $

(7)

式中：T_s为地表温度值；T_smin为NDVI相同值所对应的最低地表温度值，为T_s-NDVI特征空间的湿边；T_smax为NDVI相同值所对应的最高地表温度，为特征空间的干边。根据像元构造的特征空间，同时对干边和湿边进行线性拟合，其方程为：

$ T_{s \max }=a_{1}+b_{1} \times \mathrm{NDVI} $

(8)

$ T_{s \min }=a_{2}+b_{2} \times \mathrm{NDVI} $

(9)

式中：a₁，b₁为干边拟合方程的系数；a₂，b₂为湿边拟合方程的系数。如此TVDI指数的计算方法可以进一步转化为：

$ \mathrm{TVDI}=\frac{T_{s}-\left(a_{2}+b_{2} \times \mathrm{NDVI}\right)}{\left(a_{1}+b_{1} \times \mathrm{NDVI}\right)-\left(a_{2}+b_{2} \times \mathrm{NDVI}\right)} $

(10)

TVDI的值在0-1之间，干边对应的TVDI值为1，湿边对应的TVDI值为0，TVDI的值越趋向于0，土壤的湿润度越高，反之，土壤的湿润度越低。

3 结果与分析 3.1 土地利用/覆盖

基于监督分类和目视解译结合的方法将研究区内土地分为了：绿地、水域、建筑用地、裸地4大类。其中，绿地包括了林地、草地、耕地等；水域涵盖了江、河、湖、渠等；建筑用地指城乡住宅、公共用地等；裸地一般为无植被生长的未利用土地，总体空间分布见图 1。为探究各地类的变化规律，进一步统计了1990，2000，2010，2020年各类土地利用/覆盖信息的面积及比例情况如表 2。

由图 1可知，在时间维度上，1990年研究区内土地结构以绿地和裸地为主，建筑用地较少，水域地区更为稀缺，说明工程沿线土地资源丰富，但城镇化水平较低，城市基础设施建设较落后。1990—2010年期间，绿地面积变化幅度不大，相比之下建筑用地面积增势明显，同时裸地区域也在以肉眼可见的速度减少。2020年区域内地类主要以绿地和建筑用地为主，裸地区域被大面积开垦利用，可见当地城市化建设发展迅速。

在空间维度上，按研究区内4个地级市区段划分并对各区段地类进行统计，进一步分析得到滨州和青岛区段的水域面积较其他区段相对较大。可能受其地理位置的影响，滨州地处引黄济青工程黄河入水口，而青岛东、南濒临黄海，带来的河湖水源较充足。同时，水源和地理的优势也使得青岛区段的耕地、林地、草地发展的较好，其绿地面积比例在4个区段中一直都居于前列。建筑用地区域在4个区段中区别并不大，说明研究区内各地区城镇化发展相对较为均衡。

从表 2数据分析，1990—2020年，绿地面积比例共增长6.93%，长期占据研究区内一半以上土地；建筑用地面积比例增加18.37%，尤其以1990—2000年增势明显，增幅近15%。由此可见，这10 a的城镇化发展效果显著。相比之下水域面积增加并不明显(30 a间面积增长17.18 km²)，增幅1.61%，大部分河流湖泊位于滨州和青岛地段，主要受黄河和黄海影响。当然调水工程作为沿线重要的生态纽带也会直接促进水资源面积的增长。此外，城市的发展离不开城乡建设，工程运行以来，大面积的未利用土地被开垦翻新以开展房屋、基础设施建设，相对的裸土地区域大幅缩减，占总面积比例降低近27%，现阶段研究区内仅有2.61%的区域为裸土区。

3.2 植被指数及覆盖度 3.2.1 植被指数

植被指数的值一般在[-1, 1]区间，通常植被长势越好、覆盖程度越高对应的指数越高。整个研究区的NDVI空间分布见图 2。

利用均值法统计了不同时间研究区各地市级区段的平均NDVI指数(如图 3所示)。引黄济青工程沿线各个时期的NDVI均值分别为0.296 2, 0.335 7, 0.430 1, 0.487 4。从整体变化趋势来看，指数呈现不断上升趋势。在1990年植被指数处于较低的状态，说明当时的研究区内植被生长并不茂盛，植被的覆盖率也比较低。经过30a多的发展，植被指数均值增加了近0.2，增势明显，可见当地的生态发展进步显著。

通过对图 3的进一步分析，发现各地市级区段的NDVI均值变化趋势与沿线总趋势基本保持一致。其中，青岛区段的均值变化较为显著，大体上都高于其他地区，前文曾提及青岛区段中的水域面积也居于前列，而水资源对于植被作物的生长不可或缺。加之青岛市地理位置优越，气候适宜，对草木的生长和农作物的培育更是十分有益。此外，可以看到在1990年，滨州、潍坊、青岛3区段的NDVI均值几乎接近，而东营区段的植被指数却远低于3者，说明当时东营区段内虽然耕地、林地、草地的占地面积不少，但植被的茂密程度相对于其他区域较弱。到了2020年，4地市级区段的植被指数均值都相接近，可见当下研究区各地的生态发展均衡且逐渐趋于稳定。

3.2.2 植被覆盖度

植被覆盖度相较于植被指数，是更能衡量地表植被覆盖状况的重要指标，实际研究中常将两者结合来反映植被的生长情况。植被覆盖度的计算值在0~1之间，根据阈值范围和研究区植被覆盖特点将其分为5个等级，具体分级和各覆盖度面积比例统计情况见表 3。

据表 3可知，1990年引黄济青工程沿线面积的45.57%为中低植被覆盖程度区域，低植被覆盖程度区域次之，比例约36.91%，高植被覆盖程度区域比例最小，仅有2.08%，总体植被覆盖程度以低、中低覆盖为主。2000年，研究区内低、中低覆盖区域比例较1990年下降了约20%，中高植被覆盖区域相对增多，但总体覆盖程度仍以低、中低植被覆盖度为主。2000—2010年间，区域内的植被发育较好，覆盖度等级也发生了变化，由原先的以中低植被覆盖为主转变为以中高植被覆盖为主，其中高植被区域比例约31.01%，处于遥遥领先地位，相比之下其他覆盖等级相对均衡。现阶段，研究区已完全以高植被覆盖区域为主，且所占面积近总区域的一半。与前文中土地利用分类和植被指数的解译结果一致，说明研究区内的植被不仅区域在增加，长势也越来越茂盛。

3.3 地表温度

利用大气校正法得到的1990—2020年研究区内总体温度分布情况(如图 4所示)。为了更清晰的研究其空间分布和变化趋势，最终根据其值域范围将地表温度值划分为了6个区间，分别为低于20，20~25，25~30，30~35，36~40 ℃，高于40 ℃。

对各区域的地表温度均值进行统计(见图 5)。从图 5地温均值变化趋势可知，4个时期研究区地表温度均值分别为30.263，33.825，31.909，34.865 ℃，整体上呈现波动式变化。从滨州、东营、潍坊、青岛4市区段的地温变化来看，与总体的波动变化趋势大体一致。其中，平均温度最低的地区通常为青岛区段。原因在于，青岛临靠黄海，受季风性和海洋性气候共同作用。夏季温湿多雨，但无酷暑，较其他地区温度相对低；冬季风大温和，且持续时间较长。相比之下，潍坊市、东营市、滨州市区段地理较为内陆，受海洋性气候影响小，因此地表温度表现的比青岛区段高。

结合土地分类结果，发现地表温度值在空间分布上也有一定的差异性，水域地区温度最低，大多处于10~25 ℃区间内。植被的蒸腾对温度也有一定的降低作用，使得研究区内农田、林地、草地的地温都在25~35 ℃之间。而其余的建筑用地和裸地区域受太阳直射，较为炎热干燥，因此温度大都在35 ℃以上。

3.4 土壤干湿度

利用温度植被指数公式可以得到1990—2020年研究区总体土壤干湿度(即TVDI指数)的空间分布(图 6)。如图 6所示，1990，2000，2010，2020年研究区内的TVDI均值分别表现为0.758，0.694，0.568，0.433，由于TVDI指数的值与土壤湿度呈负相关，因此区域内的土壤呈现越来越湿润的状态。

图 7统计了1990—2020年研究区和各地级市区段土壤干湿度均值的变化趋势，分析可知，滨州、东营、潍坊、青岛市区段的TVDI均值变化呈现和总区域相同的下降趋势，即土壤湿润度越来越高。其中，青岛区段的土壤湿润度最高，滨州次之，潍坊、东营区段相对于沿海的青海和与黄河接壤的滨州，土壤都较为干旱。1990年，东营区段的TVDI指数高达0.866，土壤湿润程度远低于其他3区段，可能受到了水域、温度、降雨等多方面因素的影响，由于土壤的干湿性对于植被的生长有很大的作用，因此这也在另一方面解释了同年该区段内NDVI指数值不高的原因。

现阶段，研究区内的土壤干湿程度已相对均衡，各地级市区段的平均TVDI指数值保持在0.39~0.45的区间内。

4 结论

本文以引黄济青工程输水段沿线为研究区，采用遥感技术对生态因子的时空分布特征进行探究，得出结论如下。

(1) 研究区生态环境质量改善明显，城镇化水平发展迅速。30 a间，水域面积比例增长1.61%，绿地面积比例增长6.93%，裸地面积比例由29.52%缩减至2.61%，多用以建设基础设施和发展生态、农业等。此外，各地级市区段土地结构变化与研究区总趋势大致相同。

(2) 研究区内土壤湿润度及植被覆盖度显著提升，且现阶段发展趋于均衡。工程运行期间，土壤含水量增加，直接促进着区域内植被的繁茂生长及覆盖程度的大幅提高。其中，青岛区段植被生长状况最好，土壤也最湿润。

(3) 各生态因子间彼此不独立。温度的波动性、水域面积的增加、土壤的湿润程度必然都影响着植被的生长，而植被的生长又促使着绿地面积的增加和裸土区域的利用，说明因子间存在一定程度上的相关性。

(4) 利用遥感技术研究较大区域的生态变化规律是一种可行、实用、便捷的方法，可以推广应用于类似研究和实践。