水体动态变化信息遥感监测方法研究:以秦皇岛桃林口水库为例

崔名赫，潘军，蒋立军，刘祉余，许鹤，麻洪川

1.吉林大学 地球探测科学与技术学院，长春 130026； 2.中国人民解放军31441部队，沈阳 110001； 3.长春建筑学院 交通学院，长春 130699； 4.中国人民解放军31009部队，北京 100088

0 引言

遥感动态变化监测技术是基于同一地区的多时相遥感数据，对地表地物类型的转化/改变进行识别、属性变化信息进行获取以及遥感制图与分析的一种方法。动态变化监测作为遥感图像分析中的一项重要应用，为环境监测、土地利用和资源勘探等方面提供了有效的技术手段[1--4]。

为探究有效的动态监测方法，相关学者提出了不同的研究方法。例如，李丹等[5]以遥感图像解译技术作为支撑，将遥感影像中目标地物的各种影像特征(光谱特征，纹理特征和形状特征等)为标志，运用4个年份的多时相Landsat遥感影像实现对耕地资源空间信息的提取，分析研究区耕地分布及变化。Li et al.[6]对北京市1984—2013年的城市变化状况进行了分类与变化检测，通过目视解译选取了准确、稳定不变的训练样本，随后利用该训练样本进行了分类，并将分类结果使用到变化检测中。乔芳敏等[7]采用最大似然法对影像进行分类，再进行分类后比较，进而得出变化信息。纪超南[8]采用数据分组和直方图匹配算法对长时序影像进行分类，进而提高了变化检测的精度。

目前已有的研究方法大体上可分为两类：直接比较法和分类后比较法。前者是直接利用多时相遥感影像/数据对地物类型的光谱、空间变化特征进行增强与分类；后者是在不同时相的遥感影像上单独进行分类，然后比较图像同一位置的分类结果，进而确定变化信息的位置和类型。分类后比较的方法不必考虑不同时相图像间辐射校正和匹配等问题，但是它的精度主要取决于每个时相图像的分类精度，存在严重的误差累积问题。直接比较法大致分为以下几类：①彩色合成法：将后一时相某一波段赋予R通道、前一时相同一波段同时赋予G、B通道，变化信息一般表现为红色、青色色彩，易于识别。该方法操作简单，但精度较低，适于初步分析。②光谱特征变异法：在多时相遥感数据中，变化像元的光谱特征发生了显著变化，由此可计算两个时相数据中同一像元样本之间的欧式距离、夹角余弦等参数来作为光谱特征发生变化的度量指标，但不同地物类型变化可能会有相同或类似的度量值。③差值法、比值法，即：将两个时相的遥感数据对应波段进行差值、比值波段运算，则变化信息在差值、比值数据中具有显著区别于未变化像元的数据特征。与上述“光谱特征变异法”同理，不同的变化信息可能会具有相同或相似的“差值光谱特性曲线”或“比值光谱特性曲线”，此种特殊情况会导致信息混淆。

随着多时相图像数据源可获取性的增加，人们对遥感图像变化检测及趋势分析的要求也在不断提高。首先希望高精度获取变化结果并尽可能减少预处理误差;其次希望了解变化过程,获得变化区域更丰富的属性信息;最后希望充分挖掘多时相图像的时序信息,并全面结合变化检测结果,对未知时相进行趋势预测分析[9]。一个好的变化检测方法应该回答3个方面的问题：变化的区域和程度、变化类型的分布以及变化的类型[10]。因此，如何有效利用多时相数据的光谱信息来识别和提取某种地物类型的动态变化信息，以及如何减小误差累积、提高监测精度并减少分析的复杂度，是动态变化监测中有待解决的问题[11--13]。

笔者以秦皇岛桃林口水库为例利用多时相波段序列因子分析法进行水体动态变化信息遥感监测等研究，把前后两个时相的全部波段组合成一个新的波段序列作为原始变量，对其进行因子分析，从中提取出具有明确指示意义的波段组合作为因子，并利用方差分析方法对各因子进行显著性检验，筛选出适合提取水体变化信息的因子，根据因子构建水体变化指数，对水体变化信息进行定量识别和提取。

1 多时相数据组合法及因子分析原理

1.1 多时相数据组合法

多时相数据组合法是将两个时相数据的所有波段组合为新的数据，所生成的组合数据中，各像元的反射率曲线为两个时相波谱特性曲线的组合，称为“地物波谱特性时相变化曲线”(图1)。

图1 地物波谱特性时相变化曲线示意图Fig.1 Schematic diagram of spectral characteristic time-phase change curves of ground object

其前半部分(波段1～波段7)为同名像元前一时相(变化前)的地物波谱特性曲线，后半部分(波段8～波段14)为后一时相(变化后)地物波谱特性曲线，整个曲线代表了前后两时相地物类型的变化信息。

单一时相遥感数据中的地物波谱特性曲线代表了某种具体的地物类型，而“地物波谱特性时相变化曲线”代表某种地物类型变化信息。通过地物波谱特性时相变化曲线，可以了解每种地物类型变化的光谱特征，对比不同的地物类型变化的光谱特征，可以发现每种地物类型的变化都有其各自的独特性，根据这种独特性建立解译标志，从而对此地物类型变化进行识别和提取。

1.2 因子分析

在分析地物波谱特性时相变化曲线的基础上进行因子分析，确定能够代表某种地物类型变化的因子，根据因子识别和提取此类地物类型变化信息，最后得出水体的动态变化信息。

因子分析是帮助人们对大量观测数据进行综合分析解释的一种多元统计方法[14]。其基本思想是从协方差矩阵或相关矩阵出发，根据相关性大小把原始变量分组，使同组内的变量间相关性较高，而不同组的变量间相关性较低。每组变量用一个不可观测的假想变量表示，这个假想变量称为公共因子，几个公共因子能反映原始变量的主要信息[15]。因子分析分为R型与Q型两种。R型因子分析用于研究变量间的关系，Q型因子分析用于研究样本间关系[16]。本文采用R型因子分析。其矩阵形式为：

X=FAT+EC

(1)

式中：X代表标准化数据矩阵，X的(i,j)元为多时相遥感影像中像元i波段j的反射率，X(m×n)为n行m列矩阵，n为像元数(样本数)，m为波段数(变量数，全部时相的所有波段数目)。

A为公因子负载矩阵：

A=(a1,a2,…,ap)=(aij)(m×p)

(2)

式中：aij表示第i个变量在第j个公因子上的负载，即公因子j所承载的变量i的信息量，反映公因子j与变量i的关联密切程度。

公因子得分矩阵F表示为：

F=(f1,f2,…,fp)=(fij)(m×p)

(3)

式中：fij为第j个公因子在第i个样本上的取值，即公因子j所承载的样本i的信息量，反映公因子j与样本i的关联密切程度。

另外，C和E分别表示单因子负载矩阵和单因子得分矩阵：

C=diag(c1,c2,…,cm)

(4)

E=(e1,e2,…,em)=(eij)(m×p)

(5)

式中：对角线元素ci为变量i在单因子ei上的负载(i=1,2,…,m)，即单因子ei所承载的变量i的信息量，反映单因子ei与变量i的关联密切程度。而元素eij表示单因子ej在样本i中的取值，即单因子ej所承载的样本i的信息量，反映单因子ej与样本i的关联密切程度。

根据因子负载矩阵A的求解方式，确定最佳取值。为了合理解释公因子数，在实际应用中一般将多个原始变量归结为少数公因子。其中，公因子数的选取原则如下：①信息损失无几原则。方差贡献代表了因子的贡献程度，一般选定方差贡献排序靠前的少数因子，以使得在信息损失很少的前提下解决实际问题，通常将方差贡献累计百分比阈值设为80%，85%，90%等。②因子具有实际意义。公因子负载矩阵反映公因子与变量的关联密切程度，可用以因子的专业意义进行理解，通常选择能解决那些具有明确专业意义的少数公因子来对实际问题进行合理解释。

2 研究区概况及数据预处理

研究区为秦皇岛桃林口水库，该水库位于河北省秦皇岛市西北部青龙县境内，主要为秦皇岛、唐山两市提供农业生产和城市用水，其水位变化直接影响两市的生活和经济发展。

由于近年来该流域及周边地区受人类活动影响较为显著，局部气候条件发生改变，流域内降水量明显下降，干旱现象较为严重，造成水位下降。因此，对秦皇岛桃林口水库的水体变化信息进行遥感动态变化监测，是该流域及水库水资源开发、利用、优化配置和运行管理的有力理论依据。本文的研究数据选用该地区2013年5月30日及2015年9月25日的两时相Landsat8卫星遥感影像，对原始数据进行了辐射定标、FLAASH大气校正等预处理，研究区如图2所示。

3 水体动态变化信息提取及精度检验

3.1 地物类型动态变化分析

通过对两时相的遥感影像进行解译得出，研究区有水体、植被、裸地和工矿用地共4种地物类型，2013—2015年主要变化包括：水体变为植被、水体变为水体、水体变为裸地、植被变为植被、植被变为裸地、裸地变为植被、裸地变为植被和工矿用地变为植被等共8种地物类型变化(表1)。

表1 土地利用类型变化矩阵

将两幅遥感影像按照时相的先后顺序进行合并，选取这8种地物变化类型的训练样本，并检验这8类训练样本的J--M距离和散度，结果表明，训练样本的可分离性满足要求[17]。经统计分析得出研究区内各类训练样本反射率的均值，绘制地物波谱特性时相变化曲线(图3)，根据曲线特征分析地物波谱特性。

a. 2013；b. 2015。图2 秦皇岛桃林口水库Landsat8卫星遥感影像Fig.2 Landsat8 satellite remote sensing image of Taolinkou reservoir in Qinhuangdao City

图3 地物波谱特性时相变化曲线图Fig.3 Spectral characteristic time-phase change curves of ground object

根据该曲线可以初步了解每种地物类型变化的光谱特征，水体在前一时相(2013年)中均为水体的光谱特征，其反射率在5、6、7波段明显低于其他地物；在后一时相(2015年)中，水体地类曲线表现为变化趋势，其反射率在波段12的大小依次为：水体<裸地<植被。但也存在某种地物类型变化不显著或与其他地物类型变化光谱特征相似的情况。从图3中可以看出：水体变化与其他地物变化类型变化在波段1和2的光谱特征相似，因此认为波段1和2就是无关波段，在监测中起到干扰作用。

地物波谱特性时相变化曲线虽然能够反映每种地物类型变化的光谱特征，但在监测中有干扰作用的无关波段，且对各波段所含有的目标地物动态变化信息量缺少定量评判标准，只凭借原始的地物波谱特性时相变化曲线很难准确地识别和提取具体的水体地物类型变化信息。

3.2 因子分析及显著性检验

为准确提取变化信息，在分析地物波谱特性时相变化曲线的基础上进行因子分析，通过因子分析可以对波段进行分类和归组，定义各波段的权重系数，确定能够代表水体地物类型变化的因子，根据水体变化因子识别和提取水体地物类型变化信息。

水体变化信息量隐含在各波段中的权重有所不同，根据各波段表达水体变化信息量的多少确定系数，构建水体变化因子。在对原始数据进行标准化处理的基础上，对研究区内的8种地物类型变化的训练样本在14个波段上的反射率进行R型因子分析，利用最大方差法提取旋转因子解。经过计算公因子方差，统计特征值、方差贡献以及累积方差贡献提取了方差贡献1的前3个公因子，3个因子一共表达了91.82%的信息。为直观分析3个因子的类别区分能力，将8种地物类型变化的训练样本的反射率在3个因子上的取值两两组合作为X轴和Y轴制成散点图(图4)。

从散点图中可以初步发现，因子1(F1)和因子2(F2)不能区分任何地物类型变化，而因子3(F3)可以明显区分出水体变化信息，且可以将3种具体的水体变化类型区分开。为定量检验各因子对水体变化信息的动态识别能力，采用方差分析方法进行显著性检验，引入可分性度量评价指标F，判断组间离差平方和是否显著大于组内离差平方和，以此衡量因子变化动态识别水平。F值计算如下：

(6)

式中：QA表示组间离差平方和；Qe表示组内离差平方和；fA为组间自由度；fe为组内自由度。

对水体整体变化与其他地物类型变化的整体两组数据(记为Ⅰ)、水体变化为水体与水体变化为植被的整体两组数据(记为Ⅱ)，对水体变化为水体与水体变化为裸地整体两组数据(记为Ⅲ)以及水体变化为植被与水体变化为裸地的整体两组数据(记为Ⅳ)分别进行方差分析，结果如表2所示。

表2 方差分析结果

由表2可知，在信度α=0.05下，因子F2在区分水体变为植被与水体变为裸地时的FFα，所以在F1和F3上水体变化信息与其他地物变化信息差异显著，但在3组检验中，F3的F值都远远大于F1的F值，表明F3对水体变化信息的动态识别效果最好。

3.3 水体变化因子选定及水体变化信息提取

根据因子分析及显著性检验的结果，确定F3为水体变化因子。并通过计算旋转成分矩阵得出F3的表达结果：

F3≈0.87×B5+0.71×B6+0.74×B12

(7)

式中：B5、B6、B12分别表示组合的多时相遥感数据的5、6、12波段。

基于公式(7)计算了研究区内8种地物类型(A～H)变化样本的因子取值(图5)。

a. F1-F2；b. F1-F3；c. F2-F3。图4 因子组合变化信息分析图Fig.4 Analysis diagrams of factor combination change information

图5 因子取值正态分布曲线Fig.5 Normal distribution curves of factor values

a. 详图；b. 全景图。图6 水体动态变化信息图Fig.6 Water dynamic change information

由图5可直观看出，水体变化因子(F3)对于3种水体变化类型区分效果较好。水体地物类型变化的值均较小，与其他地物类型变化的值有较大差异，且3种具体的水体地物类型变化也能相互明显区分，因子值大小为：水体变为水体<水体变为裸地<水体变为植被。

基于F3水体变化因子，在计算出水体变化信息识别阈值(均值加减2倍的标准差)的基础上，对研究区整幅图像进行水体变化信息提取，包括水体变为水体(蓝色)、水体变为植被(红色)、水体变为裸地(黄色)(图6)。

3.4 精度分析

利用选取的217个水体变为水体样本，162个水体变为植被样本和97个水体变为裸地样本进行提取结果的精度评价，基于提取结果并参考实际变化样本，统计提取精度(表3)。

表3 水体变化信息提取混淆矩阵

由混淆矩阵计算得到总体精度为96.8%，Kappa系数为0.949。Kappa系数在0.81～1范围内表明所分结果与实际水体变化类型具有高度的一致性。由以上结果说明本文方法的水体提取精度较高，可用于水体变化信息的识别与提取。

4 结论

(1)通过多时相数据组合法将变化前后两地物的光谱特性同时体现在“地物波谱特性时相变化曲线”上，可以初步了解每种地物类型变化的光谱特征。

(2)利用因子分析的理论手段，通过分析多时相遥感数据中波段序列对水体动态变化的响应程度，得出监测水体变化信息的3种波段组合方式(因子)，3个因子一共表达了91.82%的信息。利用方差分析方法对各因子提取水体变化信息的动态识别能力进行显著性检验，得出因子F3的F值为3 373.44，远远大于其他两个因子的F值，表明F3对水体变化信息的动态识别效果最好，可以作为水体变化识别指数。

(3)经精度验证分析，水体变化识别指数F3对水体变化信息的识别精度达到96.8 %，Kappa系数为0.949。Kappa系数在0.81～1范围内,表明所分结果与实际水体变化类型具有高度的一致性。由以上结果说明本文方法的水体提取精度较高，可用于水体变化信息的识别与提取。