基于遥感的地物目标检测是评价流域生态环境可持续发展的重要方法, 其检测方法多样, 检测精度仍有提高空间。土地利用变化、城市发展与生态环境变化等关系到人民生活和人类社会可持续发展的一系列问题, 集中起来可归结为地物目标的变化检测问题^[1-2]。遥感技术因具有时效性强、覆盖范围广、可重复观测等特点, 可以定量监测地表发生变化的位置、过程和幅度^[3], 而成为变化检测最主要的技术手段。变化检测也是目前遥感应用研究中的热点^[4-5]。遥感影像变化检测是从不同时期的遥感数据中, 定量地分析和确定地表变化的特征与过程^[6]。近些年来, 相关学者相继提出了许多基于遥感影像的变化检测方法, 并从不同的角度总结和归类。Singh等^[7]根据是否对遥感影像分类将变化检测方法分为直接比较法和分类后比较法；Rensink等^[8]根据研究对象层次的不同将变化检测方法分为基于像素级、基于特征级、基于决策级3个层次；Lu等^[9]按照采用数学方法的不同将变化检测技术分为算术运算法、变化法、分类比较法、高级模型法、GIS集成法、视觉分析法和其他方法等七类。然而, 从检测变化的本质策略上可将变化检测方法分为双时相变化检测和时间序列变化检测两大类^{[1, 10]}。前者是对两个时相遥感影像的变化检测, 常用的方法有影像代数法、主成分分析法、变化矢量分析法、纹理特征分析法、马尔科夫随机场模型法、数据融合法等。如：邓小炼等^[11]利用变化矢量分析法检测了华北平原的土地利用覆盖变化, 取得了较好的效果；Tsarouchi等^[12]人利用主成分分析法对TM影像进行了变化检测处理和检测结果分析；Brozzone等^[13]将马尔可夫随机场模型法引入到变化检测中, 提取了两时相影像中的空间邻域信息, 结果显示对噪声起到了很好的抑制作用, 优化了检测效果；Camps-Valls等^[14]将数据融合法应用到多时相遥感影像的变化检测中, 并客观评价了检测结果。

长时序植被覆盖的连续动态变化检测是生态学研究的重要领域和全球变化研究的热点。遥感技术与信息论的发展为时间序列的变化检测奠定了理论与数据基础。随着遥感数据的积累, 长时间序列的遥感影像完整地记录了地表的变化过程。然而, 目前的研究多为双时相变化检测^[15]。通过时间序列变化检测可以充分地挖掘地物在时间维的变化信息, 能够更好地反映其时空变化规律^[16-17]。时间序列数据在地表覆盖变化的研究中有广泛的应用, 其中最为直接的是监测地表植被的变化^[18-19]。常用的植被覆盖变化检测方法主要有回归分析法、主成分分析法、傅里叶变换、小波变换等。回归分析法是考虑变量之间统计联系的一种重要方法, 目前在植被覆盖的长时序变化趋势研究中应用最为广泛^[20-21]。该方法的优点是计算过程简单、清晰, 结果直观。但是要求研究对象需满足正态分布, 且对遥感影像的噪声、误差规避能力较差^[22]；此外, 该方法也很难客观地反映出长时序植被覆盖的变化强度特征。主成分分析法、傅里叶变换及小波变换类似, 是通过对时序遥感数据进行影像变换或数学处理来增强变化信息, 然后从变换分量中检测变化信息^[23], 此类方法可以有效地突出变化信息, 能够较好地表达植被覆盖的宏观变化特征, 但通过变换处理得到的影像意义不够明确, 且对长时序植被覆盖在时间维的具体变化特征难以体现出来^[24]。鉴于常用分析方法在长时间序列植被覆盖变化检测中存在的不足, 本文通过引入信息论, 提出了时间信息熵的概念及计算方法, 希望为该领域的研究提供新思维、新方法。

时间信息熵是检测植被覆盖的时序变化特征的有效方法。“熵”最早起源于物理学, 用来反映热力学系统的无序程度。1948年, 香浓将这一概念引入通信领域并称之为信息熵, 它反映了信息源的平均不确定性^[25]。如今, 信息熵已被广泛应用到生态学、地理学、社会科学等领域^[26-27]。信息熵在遥感影像变化检测领域中也有较多的应用, 但多是对两个时相变化信息的检测。如：李亮等^[28]提出了利用遥感影像的像斑差熵来检测两期地表覆盖变化情况；任晖等^[29]提出了土地利用异动熵来提取两期影像的土地利用变化信息, 取得了较好的效果, 且该方法对所用遥感数据源的配准精度、分辨率等要求较低。本文在已有研究成果的基础上, 基于信息论, 提出了利用时间信息熵方法来深入挖掘时间序列遥感影像的植被覆盖变化信息, 综合、定量化地表征其在时间维的变化规律。并以延河流域为试验区, 应用时间信息熵方法来客观反映延河流域长时间序列植被覆盖的时空变化特征, 以期为该区域生态环境的保育和治理提供决策支撑。

1 试验区概况与数据来源 1.1 试验区概况

延河流域生态环境问题较为突出, 是陕西省生态环境建设的重点区域。延河是黄河的一级支流, 发源于榆林市靖边县天赐湾乡, 由西北向东南依次流经靖边县、志丹县、安塞区、宝塔区和延长县, 最后在延长县南河沟乡凉水岸附近注入黄河, 干流总长286.9km。延河流域位于陕北黄土高原腹地, 地理位置36°23′—37°17′N, 108°45′—110°28′E, 流域总面积约7687km²(图 1)。地势起伏较大, 海拔范围475—1787m, 平均海拔1213m。截止2010年底, 流域总人口约95万, 其中农业人口占61%, GDP总量为382.7亿元^[30]。宝塔区是我国著名的红色旅游胜地。安塞区矿产资源丰富, 石油产业是其支柱产业。该流域属于典型的暖温带大陆性半干旱季风气候, 冬季寒冷干燥, 夏季高温多雨。从西北向东南, 降水、温度具有明显的梯度变化特征。年降水量少且季节分配不均匀, 年均降雨量为520mm, 年均气温约8.8—10.2℃。该区域的植被覆盖状况从南向北呈现着较为明显的变化^[31]。科学、客观地反映其长时间范围植被覆盖的变化情况, 可为该区域生态环境的建设和治理提供决策信息。

1.2 数据来源

归一化植被指数(NDVI)是本研究的重要数据源, 其能够很好地反映地表的植被覆盖情况。NDVI是遥感卫星监测植被应用最为广泛的参数^[32]。研究中所采用的NDVI数据来源于美国NASA(https://ladsweb.nascom.nasa.gov/data)的MODIS植被指数产品MOD13Q1, 空间分辨率为250m, 时间分辨率为16d。该数据产品统一采用SIN投影并已进行去云、辐射校正等预处理。研究所选用的时间范围为2000年2月到2010年12月。利用MODLAND提供的MRT(MODIS Reprojection Tools)软件对该数据产品进行格式转换和投影转换, 将HDF格式转换为Tiff格式, 并将SIN投影转换为WGS84/Albers Equal Area Conic投影, 重采样方法为最邻近法, 重采样分辨率为250m。然后利用MVC(Maximum Value Composites)法对各年23期数据进行最大值合成, 来进一步消除云、大气和太阳高度角等因素对遥感影像产生的影响^[33], 得到试验区2000—2010年逐年的NDVI年最大值数据。

2 研究方法 2.1 信息熵

信息熵(也称狭义信息熵)的计算方法基于离散型随机变量, 其计算公式见式(1)。

(1)

式中, x_i表示随机变量X的状态或取值；p(x_i)为相应的概率。

连续型随机变量的信息熵不是对上式简单地取极限值, 而是以微分熵的形式定义的^[34], 见式(2)。

(2)

式中, f(x)为概率分布函数。

由上面的计算公式可知, 求信息熵最关键的是确定概率密度或者具体的概率分布函数。Vasicek等^[35]提出了直接基于观测数据, 而无需求概率密度函数或概率分布函数来计算信息熵的估计值, 并证明了该估计值依概率收敛到理论值。将式(2)表示为式(3)的形式：

(3)

式中, F^-1(p)的微分可由样本点x_i处的斜率值来近似表达, 见式(4)：

(4)

式中, y₁≤y₂≤…≤y_n, 是样本点x₁, …, x_n按照从小到大的顺序排列得到的。m为正整数, 且。进而, 信息熵可以表达为式(5)：

(5)

Ebrahimi等^[36]对上式进行了改进, 并提出了新的计算方式, 见式(6)：

(6)

其中,

式中：y₁≤y₂≤…≤y_n, 是样本点x₁, …, x_n按照从小到大的顺序排列得到的。m为正整数, 且。

同时, 证明了该方法依概率收敛到理论值H(X), 并且与式(5)相比有较小的偏差和均方根误差^[37]。

2.2 时间信息熵

在借鉴信息熵概念的基础上, 并结合遥感数据源的特点, 本文提出了新的概念和方法, 即时间信息熵、时间序列信息熵来充分挖掘植被覆盖在时间维的变化特征, 表征其在某一时期的变化强度和变化趋势信息。

时间信息熵(temporal information entropy)以公式(6)为基础, 用来反映植被覆盖在时间维的变化强度特征。应用时应考虑不同数据源的特点, 以本研究中使用的NDVI为例, 其具体计算如式(7)：

(7)

式中, y₁≤y₂≤…≤y_n, 是样本点x₁, …, x_n按照从小到大的顺序排列得到后得到的(x_i即表示像元第i年的NDVI值)。

下面对计算公式中的参数进行说明：

m是“时间频率”因子, 通过选取不同的m值(m为不超过n/2的正整数)可以反映研究对象在不同时间尺度上的变化特征。以本研究为例, 由于试验区的植被覆盖情况是用每年NDVI数据的最大值来反映的, 故研究对象变化的时间尺度为年, 计算时取m=1；

Δ是“缩放系数”, 表示对不同数据源进行标准化处理, 使得时间信息熵的计算结果具有可比性(一般Δ值可以取为基础变化单位；以本文采用的NDVI数据为例, Δ即由时序NDVI数据统计得到的基础变化单位Δ=0.02)。

通过时间信息熵可以反映某一时段内植被覆盖在时间维的变化强度情况：H值越大反映出这一段时间内其变化强度越大, 反之则表示变化强度越小。分析时间信息熵的计算方法, 由于计算时对遥感观测数据进行了排序处理, 结果只能反映出植被覆盖的变化强度信息, 而无法体现出其在时间维的变化趋势信息。基于此, 本文提出了时间序列信息熵(time-series information entropy)来反映植被覆盖在某一时间范围内的变化趋势信息。以研究中所用的NDVI数据为例, 具体的计算见式(8)：

(8)

式中, x_i表示像元第i年的NDVI值；sgn表示符号函数, ；其他参数含义同上。

H′值为正表示在某一时段内植被覆盖的变化呈增加趋势；反之则表示植被覆盖呈减少趋势。H′绝对值越大, 表明增加(或减少)的趋势越明显。

利用上述时间信息熵、时间序列信息熵可以挖掘植被覆盖在时间维的变化特征, 定量、客观地反映出长时序植被覆盖的变化强度和变化趋势信息。同时, 在计算得到上述熵值的基础上, 对其直方图进行分割, 进而得到植被覆盖的变化等级分布图, 可为区域生态环境的保育和治理提供更为方便、宏观的参考。

2.3 线性回归分析法

目前的研究中多采用线性回归分析法来反映区域长时间序列植被覆盖的变化情况^[21-22], 本研究将新建立的时间信息熵及时间序列信息熵的结果与线性回归分析方法对比, 以分析新方法的优势。

线性回归分析法以线性拟合得到的斜率值(slope)表达植被覆盖的变化强度和变化趋势。具体而言, 斜率的绝对值越大, 表示植被的变化强度越大, 反之则变化强度越小；斜率值为正, 表示植被覆盖呈增加趋势, 反之则呈减少趋势。其计算方法见式(9)：

(9)

式中, n为监测年数；NDVI_i为任意像元第i年的NDVI值。

2.4 时间信息熵方法的验证比较

为验证本研究所提新方法的科学性和有效性, 本文选取延河流域为试验区, 将此方法与目前长时序植被覆盖变化检测中最常用的线性回归分析法进行比较。长时间序列植被覆盖的变化情况与区域的生态环境等密切相关。随着人类活动、气候变化等综合因素的影响, 区域植被覆盖的变化情况越来越复杂：不仅有变化趋势比较明显的区域, 还有变化趋势虽然不显著但在时间维的变化波动性较大的区域。对于长时间序列植被覆盖的变化检测而言, 其重点内容是客观、科学地反映出研究区某段时间内的变化强度和变化趋势特征, 进而为区域的生态环境建设和管理提供辅助决策信息。基于此, 本文将以长时间序列植被覆盖的变化强度和变化趋势检测为标准, 选取代表性样本点(表 1、图 6), 来综合对比回归分析法与本文提出算法的计算结果, 验证和说明时间信息熵方法的客观性和科学性(详细的内容见3.3节)。

3 研究结果与讨论 3.1 线性回归分析法的结果

利用线性回归分析法计算得到2000—2010年延河流域植被覆盖的变化情况。结果发现, 此时间范围内研究区的植被覆盖总体呈上升趋势(图 2)。流域北部和东南部植被覆盖的上升趋势较大, 西南部的上升趋势较小。植被覆盖呈下降趋势的区域主要分布在流域的中部、西北部, 且中部植被的下降趋势较大。

利用线性回归分析法的判定系数R²来检验所得斜率值的可信程度。计算结果发现, 判定系数较大的区域主要分布在延河流域的北部、东南部, 表明这些区域NDVI回归分析斜率值有较高的可信度(图 3)。判定系数较小的区域主要分布在流域的中部、西南部、西北部, 表明这些区域的计算结果可信度较低, 并不能客观地反映植被覆盖的变化情况。

3.2 时间信息熵方法的结果

利用时间信息熵方法, 计算得到2000—2010年延河流域植被覆盖的时间信息熵、时间序列信息熵的分布情况, 分别见图 4、图 5。

结果表明, 此时间范围内流域的植被覆盖总体变化强度较大(图 4)。延河流域的土地利用类型多为耕地、草地^[30], 近十年来由于“退耕还林、还草”等政策的实施, 该区域的植被覆盖整体变化较大。而流域的南部和西南部为林地, 此时间范围内植被变化强度较弱, 故其时间信息熵值也较小。本文研究结果与延河流域的实际情况相符合。

研究表明, 2000—2010年间植被覆盖呈增加趋势的区域主要分布在流域的北部、东北部和东南部, 反映出这些区域的植被覆盖情况变好(图 5)。而流域的中部、西北部和西南部的植被覆盖主要呈减少趋势, 需要重视此类区域植被覆盖的保育工作, 避免进一步恶化。

3.3 算法对比

以长时间序列植被覆盖的变化强度和变化趋势检测为标准, 来综合对比回归分析法与本文提出算法的计算结果, 验证和说明时间信息熵方法的客观性和科学性。研究依据延河流域植被覆盖的斜率值和判定系数的分布情况, 分别从变化趋势平缓、上升和下降区域中, 选取4个代表性样本点的NDVI数据(图 6)。

对于上述样本点的实际数据(表 1), 分别利用线性回归分析法和本文提出的时间信息熵方法计算, 所得结果如表 2所示。

对于前3个样本点的实际数据, 可以容易地看出, 2000—2010年间样本点1的NDVI数据变化非常平稳, 样本点2增加的趋势明显, 样本点3减少的趋势显著。相应的, 从表 2中, 我们可以得到：由回归分析法得到的斜率值其判定系数较大。故利用线性回归分析法得到的斜率值有较高的可信度, 能够反映出这段时间内的植被覆盖变化特征。同样的, 本文提出的时间序列信息熵也反映出这些变化趋势：时间序列信息熵的绝对值越大, 表示变化趋势越明显, 且时间序列信息熵为正, 表示植被覆盖呈增加趋势；反之, 表示呈减少趋势。

对于样本点1、4, 由回归分析法得到的斜率值相同, 但其变化的波动性却有着很大的差别。对于植被覆盖的此种变化情况, 线性回归分析法难以客观地反映出时间维的变化特征。以样本点4数据为例, 利用该方法得到的斜率值的判定系数非常小(0.106), 也说明其可信度极小, 不能客观地表征植被的变化特征。随着人类活动的加剧、气候变化等综合因素的影响, 导致区域植被覆盖的变化波动较大、变化日趋复杂。以延河流域为例, 统计图 3中所示的判定系数分布情况, 可知其中46.9%的区域由回归分析法得到的斜率值其判定系数较小(小于0.65), 故利用回归分析法得到的斜率值难以客观地反映出此类区域的植被覆盖的变化规律。利用本文提出的时间信息熵方法能够很好地表征出此类变化特点：样本点4的时间信息熵值大, 反映出这段时间内NDVI变化的波动性较大；而样本点1的时间信息熵值只有0.3557, 表明其变化的波动性非常小。综上, 我们能够得到：本文提出的时间信息熵能够客观表征植被覆盖变化的波动性特征。

综合以上分析, 虽然线性回归分析法可以简单、直观地反映2000—2010年延河流域植被覆盖的变化情况。然而, 由于该方法的前提假设是植被在一段时间范围内呈线性变化特点, 对于植被变化波动较大的区域, 线性回归分析法却很难客观、科学地表征这段时间内植被覆盖的变化特征。而本文提出的时间信息熵方法表征长时间序列植被覆盖变化信息有着独特的优势, 能够客观地反映出研究区植被覆盖的变化强度和变化趋势特征。

更进一步的, 通过对熵值分布的直方图进行分割, 得到植被覆盖变化等级分布图, 可为区域生态环境的保护和治理提供更为方便、宏观的参考。以延河流域为例, 对反映植被覆盖变化强度信息的时间信息熵值分布图进行图像分割得到阈值A=1.68(来确定基本不变的区域)；对反映植被覆盖变化趋势信息的时间序列信息熵进行分割得到阈值B=1.96(判定增加的程度是否明显)、C=-0.73(判定减少的程度是否严重)来进一步反映研究区植被覆盖的变化等级。综合以上的内容, 2000—2010年延河流域植被覆盖变化情况可分为以下五种等级。

时间序列的遥感检测信息反映近10年流域生态环境偏好(图 7)。统计结果表明, 可以得到2000—2010年延河流域植被覆盖的变化等级以增加为主, 占流域面积的80.7%；明显增加的区域占流域面积13.9%, 主要分布在流域的东北部和东南部；植被覆盖减少的区域占比2.4%, 主要分布在流域的西部和西北部；严重减少的区域占比1.1%, 主要分布在流域的中部和西南部, 是需要重点的生态恢复与治理区域；基本不变的区域占流域面积的1.9%。

4 结论

本文利用时间信息熵方法计算了2000—2010年延河流域的植被覆盖变化信息, 厘清了其时空变化特征。研究表明, 延河流域的植被覆盖情况总体得到改善, 反映出该区域的“退耕还林、还草”政策得到了有效实施, 与延河流域的实际情况相符合。然而, 流域的中部、西部的植被覆盖呈现恶化趋势, 应重视此类区域植被覆盖的保育工作, 加强生态环境的修复和治理力度。

本文提出的时间信息熵方法与回归分析法相比, 能够更加客观地表征长时间序列植被覆盖的变化强度和变化趋势信息, 可为区域生态环境的保护和管理提供更为科学的决策参考。研究成果可为长时间序列遥感影像的植被覆盖变化检测提供新方法、新思路。研究旨在揭示某一时间范围内区域植被覆盖的整体变化特征, 尚未涉及其中的“突变”问题, 未来研究需要利用时间信息熵方法, 以期揭示植被覆盖时空变化的影响机制。