运城市盐湖区31年植被与水体指标变化及驱动力分析

刘畅，马超，2

（1.河南理工大学 测绘与国土信息工程学院，河南 焦作 454000；2.河南理工大学 自然资源部矿山时空信息与生态修复重点实验室，河南 焦作 454000）

0 引言

运城盐湖，又名“银湖”“解池”和“河东盐池”，诞生于新生代第三纪喜马拉雅构造运动时期，是世界三大硫酸钠型内陆盐湖之一［1］。运城盐湖位于晋陕豫三省交界处，是中国最早开发食用石盐的古老盐湖，它的产盐技术——硝板上产盐，是全世界最古老最独特的制盐工艺技术。盐湖在中国化工工业有着十分特殊的地位，许多学者从盐湖的化学工业、资源开发、美容护肤等方面进行了研究，王玉萍等［2］对盐湖苦卤制取精制硫酸镁的工艺进行了研究；孙培霞［3］对运城盐湖的资源开发利用情况进行了研究；杜首英等［4-5］研究了盐湖水中的黑泥、矿物质、盐藻等对皮肤的作用。但重视资源开发的同时，对区域山水林田湖草的关注较少，运城盐湖的植被生态环境及盐湖水域变化的气候响应研究较少。

随着社会快速发展和城镇化进程加快，生态问题成为制约社会发展的重要问题。植被和水体是生物圈的重要组成部分，是地面生态系统状况以及区域和全球环境变化的直接指标［6］。遥感数据具有宏观、快捷、信息量丰富的特性，与传统测量方法相比，可以更方便，快捷提取植被生物量、水体等信息［7-10］。自20世纪七八十年代起，遥感技术已经开始应用于地表水域的动态监测。马新萍等［11］研究了秦岭山地归一化植被指数（normalized difference vegetation index，NDVI）与温度相应关系；刘少华等［12］研究了中国植被NDVI与≥10℃积温和降水量的年际变化相关关系；杨桂山等［13］结合历史调查和研究资料，查清了我国面积大于1 km2的湖泊数量，系统分析了近几十年我国湖泊面临的生态环境问题，提出了我国湖泊保护的策略；张飞等［14］分别利用归一化水体指数模型（normalized difference water index，NDW I）和修正归一化水体指数模型（modified normalized difference water index，MNDW I）对新疆艾比湖水域面积进行提取分析。Landsat系列卫星由美国地质调查局（United States Geological Survey，USGS）于1972年发布，具备较短时段、重复覆盖同一地区的能力，已连续运行近50年，积累了丰富的遥感存档数据，主要包括Landsat 1—3 MSS数据，Landsat 4—5 TM 数据，Landsat 7 ETM＋数据和Landsat 8 OLI数据，已经成为湿地植被及水体的时空变化分析的主要数据源［15-18］。

运城市农业区高度发达，盐湖区具有丰富的盐化工资源，近年来盐池化工资源的过度开发，对植被和水环境造成了很大的破坏。本文拟利用Landsat TM／ETM＋／OLI数据，针对密切联系农业种植业的归一化差值植被指数NDVI和盐业化工高度关联的归一化差值水体指数NDWI进行长时序分析，研究时间序列陆地植被生物量变化和运城盐湖面积及变化趋势，对运城市盐湖区水域、陆表植被生态环境进行评估，以期为协调盐湖区的开发与保护提供信息支撑。

1 研究区与研究数据

1.1 研究区概况

运城盐湖位于涑水盆地东南部，运城盆地的最低处和中条山北麓山前洪积平原与涑水河冲积平原的交汇处。涑水盆地为运城盐湖的汇水区，其北、东、南三面分别为孤山、稷王山、紫金山以及中条山等山脉，标高比湖北侧城区低20余米，比西侧的黄河低14 m，是典型的内陆闭流湖泊，见图1。

图1 运城市盐湖区地理位置Fig.1 Geographical location of Yanhu district in the study area

运城盐湖处于暖温带大陆性气候南缘，属温暖半干旱大陆性季风气候，四季分明，冬寒夏热，春冬多风，年降雨量500～600 mm，年平均蒸发量2 258 mm，年均气温13.61℃，无霜期208 d，有利于盐、芒硝等矿物质的析出［1］。运城盐湖素有“中国死海”之称，虽有死海之名，却不像死海一样寸草不生，相反，盐湖中生长着大量的浮游植物，每逢夏季，千年盐湖在阳光的照射下，化身为七彩湖，冬季则成为大量鸟类南迁的中转站。

1.2 研究数据

本次研究采用1987—2017年美国陆地卫星（Landsat-5／7／8）的遥感影像（http：／／www.usgs.gov），成像时间分别 为1987-08-30，1990-08-22，1992-07-26，1992-10-20，1995-09-21，1997-08-09，1999-08-23，2002-08-31，2005-08-31，2008-08-23，2011-09-01，2014-08-24，2017-08-16。传感器包括Landsat 5 TM，Landsat 7 ETM ＋和Landsat 8 OLI，以TM和ETM＋数据为主。其中，Landsat7 ETM＋除增加分辨率为15 m的全色波段外，其余波段与Landsat 5在空间分辨率与光谱特性上保持一致。Landsat 8 OLI包含了ETM＋传感器的所有波段，同时对波段重新进行调整，如波段5排除了0.825μm处水汽吸收特征；缩减波段8的范围较窄；新增蓝色波段和短波红外波段。卫星波段特征参数和主要用途见表1。

表1 Landsat TM／ETM＋／OLI传感器波段特征和主要用途Tab.1 Band characteristics and main uses of Landsat TM／ETM＋／OLI sensors

受卫星轨道和云层影响，同时相数据难以获取，综合考虑研究区域的植物生长状况、天气、云层覆盖等情况，从可获取的数据序列中，筛选出植被长势较好的夏秋季遥感影像，进行NDVI，NDW I等信息提取。

此外，研究还采用了1987—2015年国家气象数据中心全国的气象资料（http：／／data.cma.cn／）和国家科技基础条件平台-国家地球系统科学数据共享服务平台-黄土高原科学数据中心提供的黄土高原1 km分辨率逐月降雨量数据集（2001—2014年）（http：／／loess.geodata.cn／index.htm l）。

2 研究方法

2.1 遥感指数

归一化差异植被指数（NDVI）是反映农作物长势和营养信息的重要参数之一。在遥感影像中，近红外波段的反射率与红光波段的反射率之差比两者之和。NDVI数学模型为

NDVI取［－1，1］，负值表示地面覆盖云、水、雪等，对可见光高反射；0表示有岩石或裸土等，反射率ρnir和ρred近似相等；正值表示有植被覆盖，且随覆盖度增大而增大。

归一化水体指数（NDW I）是对遥感影像的绿波段与近红外波段进行归一化差值处理，以凸显影像中的水体信息［19］。其表达式为

式中：ρgreen和ρnir分别为绿波段和近红外波段的反射率。

2.2 变异系数

在概率论和统计学中，变异系数（Coefficient of variation，Cv）常常被用于描述数据离散程度。Cv可以认为是反映数据离散程度的绝对值，因其没有量纲，可以进行客观比较。变异系数是数据的标准差和数据平均数之比，表示数据的相对变化程度［20］。其计算公式为

式中：Cv为NDVI或NDWI的变异系数；xi为第i年的NDVI或NDWI；为某时间段NDVI或NDWI的平均值；Cv值越小，表明NDVI或NDWI波动越小；Cv值越大，表明NDVI或NDWI波动越大。利用变异系数分析，得到研究区内各像元遥感指数的波动程度，获取研究区的变化趋势［21］。

3 结果与分析

3.1 植被指数变化

选取近似相等间隔（2～3 a），植被生长最好时段（7—9月）的Landsat多光谱数据。通过辐射定标、大气校正等处理，将图像全部转化为大气层顶的表观反射率图像，有效降低了天气和传感器敏感度的影响，再通过波段运算获得12期NDVI图像（图2）。

图2 研究区时间序列NDVI密度分割结果Fig.2 NDVI density segmentation results of time series in the study area

对获得的NDVI结果进行密度分割，以便从各级占比获得一些宏观认识［22］。通过遥感影像与NDVI图像交互采样分析，对其分为5级，即水域、裸地及沙地、稀疏草地、草地及湿草地、灌丛及林地，分别对应于的NDVI为［－1.0，0.0］，（0.0，0.2］，（0.2，0.4］，（0.4，0.7］，（0.7，1.0］。

进一步对31年研究区域的NDVI总值进行植被类型量化统计，可得百分比柱状图（图3）。NDVI低于0的主要为水域（其灵敏度低于NDWI，仅供参考），平均占研究区域总面积的2.73%；介于（0.0，0.2］的裸地及沙地占研究区总面积的17.45%；介于（0.2，0.4］的稀疏草地占研究区域总面积的27.96%；介于（0.4，0.7］的草地及湿草地占研究区总面积的42.88%；介于（0.7，1.0］的灌丛及林地占研究区总面积的8.99%。

由图3可以看出，运城市盐湖区的植被结构趋于优化，高NDVI的植被增加，31年总NDVI值整体呈上升趋势，对总NDVI值和年份进行归一化处理，得到上升速率slope＝20.509，表明1987—2017年期间NDVI总量逐年增加，区域植被生态环境向好。

图3 NDVI密度分割百分比柱状图Fig.3 Percentage columns of NDVI density slice

31年中，水域面积最大为1987年，占研究区总面积的3.65%；裸地及沙地面积最大为1997年，占研究区总面积的35.57%，同时1997年也是水域面积占比最低的一年；草地及湿地面积最大为1990年，占研究区总面积的55.74%；2005，2008和2014年灌丛及草地占比相对较高，2014年最大，占研究区总面积的34.32%。2011年运城地区7—9月汛期发生强降雨，遭受洪涝灾害，农作物和农林牧渔业受灾严重［23］，当期NDVI偏低。其中2005，2008，2011和2014年的灌丛及林地，草地及湿地的变化幅度过大，其根本原因还是降水。考虑到植被动态对气候的响应以及降水的滞后性，提取运城盐湖2002年7月，2005年7月，2008年7月，2011年8月，2014年7月的降水量，由图4可以看出，时间滞后一个月的降水量与NDVI密切相关（SPSSⒸr＝0.777，P＜0.01），即一年中如果降水集中出现在卫星成像之前的几个月，会增加NDVI高值的比例，如果降水平均分摊在卫星成像前后的几个月，则NDVI变化幅度相对稳定，不会出现变化幅度过大的情况。

图4 2002—2014年研究区月降水量-NDVI变化趋势Fig.4 Trends of precipitation-NDVI in the study area during 2002—2014

3.2 水体指数变化

水体反射率从可见光到近红外依次降低，在近红外波段反射率几乎为0；植被在近红外波段反射率很高，用绿光波段和近红外波段反射率反差生成的NDW I进行水体提取，该NDW I可以很好地抑制植被。利用NDW I模型，以0为阈值，提取NDWI为正值的区域，获取研究区1987—2017年共12期的NDWI图像（图5）。

图5 研究区时间序列NDWI水体提取结果Fig.5 Time series NDWI water extraction results

NDW I用于对运城市盐湖区水域面积的提取，主要目的是对盐池动态的研究，通过NDWI阈值分析，获取1987—2017年31年盐湖水域面积变化情况（表2）。

表2 盐湖水域面积统计Tab.2 Area statistics of Yanhu lake area

运城盐池总面积由硝池和盐池（含鸭子池）组成。由时间序列水域面积统计分析知，近30年水域面积波动明显，总面积缩小了16.1%；其中，盐池水域面积缩小了19.9%，硝池水域面积缩小了6.2%，硝池在1997年前后数年几乎干涸。由图5还可以看出，枯水年份，湖域面积缩小，分割成的水域数量增加，湖域周长增加；丰水年份，湖域面积增加，分割成的水域合并之后数量减小，湖域周长减小。湖水的盈亏直接表现在晒盐等水工建筑的完整性受到了破坏，进一步影响盐、硝的产量。

3.3 气候变化与相关性分析

3.3.1 气候变化分析

运城市盐湖区范围较小，采用单站点数据精度将优于多站点数据插值，故采取运城站的气象数据作为研究数据，气候变化分析选取运城市盐湖区的年平均气温和年降水量两个气象参数。分析1987—2015年气候数据，可得：

（1）1987—2015年运城市盐湖区年降水量总体呈下降趋势，平均下降速率为2.29 mm／a，如图6（a）所示。降水波动范围较大，有4期较为明显的相对低值，即1991年（381mm）、1997 年（285.3 mm）、2012年（333 mm）、2013年（334 mm）；有7期较为明显的相对高值，即1989年（558 mm）、1996年（603 mm）、1998年（658 mm）、2007年（637 mm）、2011年（732 mm）、2003年（849.2 mm）；在 1997 年（285.3 mm）和 2003 年（849.2 mm）达到最低值和最高值两个极值。对数据进行阶段性统计分析，其结果显示1987—1990年的降水量均值为558.1 mm，1991—2000年的降水量均值为488.8 mm，2001—2010年10年的降水量均值为505.3 mm，2010—2015年的降水量均值为468.9 mm，整体呈下降趋势。

图6 1987—2014年研究区降水和气温变化趋势Fig.6 Trends of precipitation and temperature in the study area during 1987—2014

（2）运城市盐湖区1987—2015年整体年均温度总体呈上升趋势，平均上升速率为0.013℃／a，温度波动在13.6～15.5℃间，如图（图6（b））所示。在1993年（13.6℃），和1996年（13.6℃）出现相对较低值，同时在1999年（15.4℃），2002年（15.5℃），2013年（15.4℃）出现相对较高值，无较明显极值。对数据进行阶段性统计分析，结果显示1987—1990年4年温度均值为14.0℃，1991—2000年10年温度均值为14.5℃，2001—2010年10年温度均值为14.4℃，2010—2015年5年温度均值为14.5 ℃，1991—2000年，2001—2010年和2011—2015年3个阶段温度变化波动范围较小，整体呈上升趋势。

3.3.2 相关性分析

（1）NDVI与气候变化相关性。将12期年度NDVI总值与年均气温及降水作相关分析（图7）。考虑到气候因子对植被的影响具有滞后效应，对气温和降水两个气候因子与盐湖区NDVI值的相关性进行分析时，除NDVI与同时期气候因子做相关性分析外，还对气候因子做滞后分析。

图7 1987—2015年研究区降水-NDVI和气温-NDVI变化趋势Fig.7 Trends in precipitation-NDVI and temperature-NDVI in the study area during 1987—2015

降水量与NDVI值相关性低，差异性不显著，认为运城市盐湖区降水量与NDVI值无显著相关性。对气温与NDVI值相关性进行分析，在所做对比分析中，滞后一年的气温与NDVI值相关性最高，为0.395（P＜0.05），结果表明31年运城市盐湖区气温与NDVI值弱相关。

（2）水域面积与气候变化相关性。将时间序列气温、降水两个气候因子与盐湖区面积进行相关性分析，1987—2015年气温、降水与盐湖区面积联立制图（图8）。

图8 1987—2014年研究区降水-水域面积和气温-水域面积变化趋势Fig.8 Trends of precipitation-water area and temperature-water area in the study area during 1987—2014

根据SPSS对盐湖总面积和降水量、气温进行Pearson相关性分析，分析时间序列盐湖面积变化与降水量之间的相关系数：

说明盐湖面积与年降水量有直接关系，与年均气温无直接关系。研究表明，盐湖水资源补给是多源的：一是大气降水的直接补给保证了盐湖一半的水量；二是地下水占盐湖水量三分之一的补给量；三是盐湖周边山岭丘陵的地表水径流补给［1］。逢雨季或暴雨时节，山岭丘陵的地表积水或山洪等水必然流向地势最低的盐湖，这些水源补给保证了运城盐湖形成并存在至今。

据盐湖气象站1995—2005年11年实测记录，湖内年均降水量为446.52 mm，年均蒸发量为2 273.9 mm。盐湖气候干燥多风，少雨雪，利于湖水蒸发浓缩和芒硝等盐类的结晶析出；尤其是盐湖东南有高峻的中条山，经常发生焚风，每年5级以上大风天多达86天，形成运城盐湖独特的小气候环境，蒸发量大，有利于湖水浓缩。2002年运城盐湖面积最大，达到55.08 km2。1997年盐湖水域总面积达到最低值为22.08 km2，该年度降水量达到30年间最低值（285.3 mm），这与1997年黄河断流，山西大旱相对应［24］。

3.4 变异系数分析

根据变异系数模型（式（3）），逐像元计算了1987—2017年研究区NDVI和NDWI的变异系数并对其进行分级处理（图9）。

由图9可以看出，NDVI和NDW I波动变化剧烈的地区多集中在盐池区域与部分城镇区域。植被覆盖度区域虽有波动变化，但其变化较为平缓。城镇区及盐湖水域区NDVI和NDW I变化较为剧烈。根据Cv的大小将其分为5个等级，即低波动变化、相对较低波动变化、中等波动变化、相对较高波动变化和高波动变化（表3）。

表3 研究区变异系数统计Tab.3 Statistics of coefficient of variation in the study area

NDVI高波动变化，即Cv＞1的值占4.37%（主要集中在水域面积区域，图8（a）），相对较低波动变化，即Cv＜0.2的值占研究区面积的7.15%，NDVI集中在相对较低波动变化和中等波动变化，即0.2＜Cv≤0.6，总数占研究区面积的80.68%。

虽然运城市盐湖区NDVI和NDWI高波动变化分别仅占整个研究区的4.37%和3.39%，但却分别占水域面积的49.73%和45.19%（表3），盐湖NDVI和NDWI年际变化剧烈。

4 结论

（1）1987—2017年运城市盐湖区NDVI总值波动较为明显，总体呈上升趋势，除个别年份外，植被生态环境呈逐年改善趋势，归一化上升速率slope＝20.509。NDVI最大值出现在2008年，为68.1×104（无量纲），最小值出现在1995年，为34.2×104（无量纲），波动幅度达49.7%。

（2）采用NDW I提取运城盐湖水域面积，该值波动变化较为剧烈，在2002年出现最大值，为55.08 km2，最小值出现在1997年，为22.08 km2，波动幅度达59.9%。

（3）降水对盐湖水域面积起着决定性作用。1987—2014年盐湖水域面积与降水相关系数为0.746 4（P＜0.01），盐湖水域面积与温度不相关。31年运城盐湖区NDVI与降水无显著相关性，与气温有一定相关性（r＝0.395，P＜0.05）。

（4）NDVI变异系数波动变化剧烈的地区部分出现在城镇，与城市扩张相关；相对较低波动变化和中等波动变化，即0.2＜Cv≤0.6主要分布在农村地区，占研究区面积的80.68%，与耕作方式有关。湖面NDVI和NDW I变异系数均为高波动变化，NDVI高波动变化占盐湖总面积的49.73%，NDW I高波动变化占盐湖总面积的45.19%，盐湖NDVI和NDW I年际年变化剧烈。盐湖区NDVI关联的农业种植业和NDW I关联的盐业化工生产活跃。

因有关数据无法获取，本研究没能分析NDVI与当地粮食作物的产量以及NDW I与当地盐、硝产量的关系。