采煤沉陷区高分遥感数据融合方法

占惠珠，尚 慧，甘智慧

(西安科技大学 地质与环境学院，陕西 西安 710054)

0 引 言

采煤沉陷区地表环境监测是针对采煤引发的地质环境问题进行治理恢复过程中的重要工作。当前，矿产资源的开发利用对中国地形地貌、土地资源分布及相关工程活动造成一定影响，而采用遥感技术对矿山地表环境进行实时、准确、有效地监测已成为必然趋势[1-5]。在遥感影像处理技术中，影像融合技术可以通过融合具有不同光谱特性、不同分辨率的影像数据来提高影像质量，以便准确获取因矿山开采被破坏和压占的土地、植被以及其他建筑物信息。因此，探寻一种适合于采煤沉陷区地质环境监测的遥感影像融合方法显得尤为重要。随着高分辨率商业遥感卫星影像获取途径的增加，传统中低分辨率遥感卫星影像在精准调查和监测等方面的局限性日益显现，与此同时，针对小范围地质环境细节信息提取精度不足的问题也逐渐获得解决：LI等综合运用HIS变换、主成分变换和小波融合等多种方法，借助Spot-5遥感卫星影像全色、多光谱及短波红外波段，针对矿区生态环境要素信息提取开展了相应研究，研究成果为矿区地质环境评价和开垦提供了基础数据[6]；HESHAM等基于 HIS变换方法，融合假彩色复合Spot比率影像及高分辨率的Spot全色波段，对岩石类别进行划分，并绘制相关矿区地质图[7]；MEZNED等以突尼斯北部地区为例，融合Spot-5全色波段与Landsat ETM+多光谱数据，突出显示矿区尾矿渣堆和尾矿库的分布情况，结果表明基于多光谱和多源数据融合的遥感技术能够有效应用于尾矿渣和尾矿库的监测研究[8]；褚进海等运用Brovey，IHS，PCA和Pan Sharpening等变换方法对QuickBird影像进行融合，并对融合效果进行评价，结果显示Pan Sharpening算法是最适合QuickBird数据在露天开采矿山工作面及尾矿库遥感监测方面的融合方法[9]；王俊芳等结合QuickBird影像和Spot-5影像全色波段，采用波段组合与图像融合技术，对露天磁铁及硐采煤矿区进行地质灾害及环境污染监测研究，结果证实采用高分辨率遥感监测技术对矿山地质环境进行调查的有效性和必然性[10]；白光宇等基于Spot-5影像对吉林省辽源市煤矿进行了室内图像处理及野外实地调研，并通过分析得出Spot-5卫星影像的最优波段组合、融合算法及地物解译标志，对其他矿区应用Spot-5卫星影像开展矿山地质环境调查具有借鉴意义[11]。

虽然利用高分影像融合技术解决了小范围地质环境细节信息提取精度问题，但是研究对象主要集中在矿区污染物、地质灾害、尾矿渣和尾矿库等，所用高分数据类型较单一，且在针对采煤沉陷区地质环境监测的多源高分遥感数据融合方法的研究成果偏少。以宁夏回族自治区石嘴山市惠农采煤沉陷区为研究区，选取2003年QuickBird卫星影像(研究区进行矿山地质环境治理恢复前)和2018年高分二号卫星影像(研究区矿山地质环境现状)，采用PC Spectral Sharpening法、NNDiffuse Pan Sharpening法、Gram-Schmidt法及Brovey Sharpening法进行融合分析处理，并对处理结果进行定性与定量评价，探讨适合于采煤沉陷区地质环境监测的最优融合方法，可为采煤沉陷区地质环境监测研究提供高质量的遥感基础数据，丰富矿山地质环境演化遥感监测的理论和方法，为资源枯竭型城市的可持续发展及环境质量的改善提供借鉴。

1 地质背景及主要矿山地质环境问题

惠农采煤沉陷区地处宁夏回族自治区石嘴山市，是由石嘴山矿区(包括一矿和二矿)井下煤炭开采形成，总面积约40 km2，位于东经106°45′39″～106°47′21″，北纬 39°13′50″～39°15′44″之间(图1)。沉陷区出露地层涵盖长城系、蓟县系、二叠系、新近系、第四系，含煤地层为上石炭统太原组和下二叠统山西组，可采煤层总厚度为28.75 m。经过60 a的煤炭开采，存在的主要矿山地质环境问题包括：压占和破坏土地资源、地貌景观；矿山次生地质灾害；矿山固体废弃物对环境的污染。分布广、影响最大的矿山地质环境问题主要是占用、破坏土地资源和地貌景观，其次是矿山地质灾害，其中，地质灾害发育最严重的为地面塌陷和地裂缝，沉陷区共发育7个较大的塌陷坑，地表最大下沉值超过26 m，一般深8～10 m，塌陷坑总面积达9.1 km2[12]。除此之外，区内还有由采矿引起的地面塌陷22处，面积约1.3 km2；地裂缝140余条，总长15 462 m，裂缝影响面积0.8 km2。沉陷区内地表岩体破碎，植被稀疏，砂石散布，污水横流，甚至一些区域还成为了城市垃圾堆积场，环境污染异常严峻，给沉陷区及其周边人居生活环境造成极大影响。伴随城市的进一步扩张，采煤沉陷区日愈趋近城市中心，污染严重、管理缺位的采煤沉陷区与惠农城区形成了强烈对比，为当地居民的生产生活带来不便。从2004年开始，惠农采煤沉陷区进行了近11 a的治理，治理面积约40 km2(图1)，目前矿区的环境得到很大改善，基本消除了矿区地质灾害隐患，大幅度提高了植被覆盖率，减少了水土流失，使荒弃的土地成为城市绿色屏障，建成了集湖面、园林和塌陷区地貌景观为一体的特色园区，改善了矿区居民生存环境，提升了矿区居民生活质量。

图1 采煤沉陷区位置及范围

2 数据源及其预处理

2.1 数据源选择

高分辨率遥感影像相比于传统中低分辨率遥感影像具有空间分辨率高、数据规模大等特点，可从影像中获取更多、更准确的地物信息[13]。因此，选用了符合实验标准且含云量较少的惠农采煤沉陷区2003年QuickBird影像(研究区进行矿山地质环境治理恢复前)及2018年高分二号影像(研究区矿山地质环境现状)进行融合方法研究。其中，QuickBird影像包括了分辨率为0.61 m的全色影像和2.44 m的多光谱影像，全色波长范围0.405～1.053 μm，多光谱影像的波长范围为0.430～0.545 μm(蓝色波段)、0.466～0.620 μm(绿色波段)、0.59～0.71 μm(红色波段)、0.715～0.918 μm(近红外波段)；高分二号全色影像的空间分辨率为0.80 m，波长范围：0.45～0.90 μm；多光谱影像的空间分辨率为3.20 m，包含了4个波段，波长范围分别为：0.45～0.52 μm(蓝色波段)、0.52～0.59 μm(绿色波段)、0.63～0.69 μm(红色波段)、0.77～0.89 μm(近红外波段)。

2.2 图像预处理

与传统中低分辨率影像相比，相同面积的高分辨率影像包含了更多、更复杂的空间信息，同时影像上同物异谱或者同谱异物的问题也更严重，因此在融合前，对原始多光谱波段与全色波段影像进行了几何纠正和多光谱影像的重采样等预处理，其中QuickBird影像的配准误差为0.215个像元，高分二号影像的配准误差为0.238个像元，满足影像融合的精度要求。

3 高分辨率遥感影像融合

3.1 融合算法

通过对比分析传统遥感影像融合方法[14-18]，根据QuickBird影像和高分二号影像的特点，结合采煤沉陷区地质环境监测应用的需要，分别选取光谱信息和空间信息保持相对较好的PC Spectral Sharpening，NNDiffuse Pan Sharpening，Gram-Schmidt和Brovey Sharpening4种融合算法对研究区2003年QuickBird影像和2018年高分二号影像进行了融合处理，通过定性评价和定量评价得到最优融合算法。

3.1.1 PCA融合算法

主成分分析法(principal component analyst，PCA)是基于像元的融合算法，又称Karhunen-Loeve(K-L)变换。首先计算多光谱影像的系数矩阵的特征值与特征向量，得到多光谱影像的各主成分分量，并对第一主成分进行提取；其次将全色波段影像拉伸到与提取的多光谱影像第一主成分的均值和方差相同；最后用经过拉伸后的全色波段影像替换第一主成分，并叠加其余主成分进行逆变换，获得融合后的影像。算法特点是可以处理3个以上的波段，省去了波段选择的过程，使融合后影像既保留了原始多光谱影像的高频信息，又具有高分辨率、高光谱的特征，同时也会对原始多光谱影像的光谱信息造成一定的丢失。

3.1.2 NNDiffuse Pan Sharpening融合算法

NNDiffuse(nearest neighbor diffusion)Pan Sharpening融合算法是基于临近扩散算法实现的[19]，主要是对高分辨率全色波段影像进行降采样，得到与多光谱影像相同的分辨率，再将每个波段相同位置的像素进行线性回归运算，得到光谱贡献量T，通过权重向量T来规范光谱，得到与原影像近似的光谱信息。算法特点是运算速度快，光谱保真度较好且对波段要求低，没有数量要求。但对融合之前的多光谱影像与全色波段影像的空间分辨率有要求，需呈整数倍，否则需对影像先进行重采样后才能实现融合。

3.1.3 Gram-Schmidt融合算法

Gram-Schmidt融合算法是基于PCA算法进行改进的一种融合算法。PCA融合影像经过主成分变换后，信息在各个主成分之间进行重新分配，第一主成分包含大部分信息，其他主成分包含的信息量依次减少，而 Gram-Schmidt融合算法以消除波段间信息的冗余性和相关性为目的，对影像的矩阵进行正交变换[20],各分量所包含的信息量大致相同，这就避免了PCA融合算法中信息过于集中在某个分量的问题。该算法在一次融合过程中可处理3个以上的波段，省去了波段选择过程，融合后的影像光谱保真度较好。由于Gram-Schmidt变换过程较为复杂，在融合过程中容易出现一些难以预料的误差。

3.1.4 Brovey Sharpening融合算法

Brovey Sharpening融合算法是基于代数的一种较为简便的融合算法，将高分辨率全色影像与归一化后的多光谱影像相乘即可完成数据融合。优点是简化了影像转换过程的系数，可以很好地提高影像上阴影、水体、城镇居民地之间的对比度，在具有高对比度的RGB影像上能够产生良好的视觉效果[21]，缺点是它只允许同时对3个波段数据进行融合运算。

3.2 定量评价指标

目前对于融合影像质量的定量评价方法没有统一的标准，不同的学者选用的评价指标不同。在前人的研究基础上，文中从融合影像的空间信息融入度[22]和光谱信息保真度[23]2个方面对不同融合算法的融合影像进行评价。在空间信息融入度方面选取了标准差、信息熵、平均梯度和相关系数4个评价指标；在光谱信息保真度方面选取了灰度均值、相对偏差和光谱扭曲度3个评价指标。

3.2.1 灰度均值

灰度均值指影像中所有像元灰度值的算术平均值，即人眼观察到的影像的平均亮度。若计算得到的融合影像的灰度均值越近似于原始多光谱影像的灰度均值，则表明该融合影像的信息保留度越好。

3.2.2 标准差

标准差代表了影像灰度级的分布状况。若计算得到的融合影像的标准差越大，则表明该影像的灰度级分布越离散，影像的反差越大，影像所包含的信息量越多。

3.2.3 信息熵

信息熵是用来反映影像信息丰富情况的一个重要指标，熵值的大小直接与影像所包含的平均信息量相关。若计算得到的融合影像的信息熵值越大，则表明影像所包含的信息量越多，融合效果越好。

3.2.4 平均梯度

平均梯度即影像的清晰度，用来衡量影像对地物的细节和纹理信息的表征能力大小，若计算得到的值越大，则表明该影像对地物的细节和纹理信息的表征能力越大，清晰度越高，融合效果越好。

3.2.5 相关系数

相关系数反映的是原始多光谱影像与融合影像之间的相关程度。若计算得到的相关系数越接近于1，则表明2幅影像的相关程度越高，融合影像保持了较多的原始多光谱影像的光谱特征，融合的效果也越好。

3.2.6 相对偏差

相对偏差是融合影像相比原始多光谱影像在光谱信息上的偏离程度。若计算得到的相对偏差值越小则说明融合影像不仅提高了空间分辨率，且较好地保持了原始多光谱影像的光谱特征。

3.2.7 光谱扭曲度

光谱扭曲度即光谱信息丢失程度，一般用融合影像与原始多光谱影像在光谱信息上的差距大小来表示。若计算得到的光谱扭曲度越大，则说明融合影像的原始光谱信息丢失越多，与原始多光谱影像的匹配度越小。

3.3 融合结果评价与分析

3.3.1 定性评价

定性评价即通过目视判别，由人眼去观察融合后影像的整体亮度、色彩反差是否合适、色调是否均匀、细节和纹理信息是否丰富、各种地物边缘是否清晰等。尽管定性评价通常会受到显示器性能、评价主体的经验与偏见等影响，但由于人眼对色彩具有强烈的感知能力且评价过程与方法具有简单、快捷的特点，其对融合后影像的评价是其他方法所无法比拟的。

由QuickBird和高分二号影像以真彩色RGB波段组合，2%的线型拉伸显示的多光谱影像及几种算法融合后的影像对比可知(图 2、图3)，在空间分辨率方面，4种算法的融合影像与原始多光谱影像相比，采煤沉陷区内厂房、矸石山、植被等地物的纹理和轮廓结构在各幅融合影像上更清晰。尤其是PCA融合影像与NNDiffuse融合影像上的地物纹理细节特征更加清晰丰富，更易于辨别；在光谱保真度方面，前3种算法的融合影像色彩与原始多光谱影像较为接近，较好地保留了原始多光谱影像的光谱信息，其中QuickBird影像数据的GS融合影像和PCA融合影像色彩的亮度、饱和度和对比度有很大程度提高，可准确判读治理前煤矸石山的范围与形态，矿区建筑物边界与道路清晰可见。而NNDiffuse和 Brovey融合影像的融合效果相对较差；高分二号PCA与 NNDiffuse融合影像与原始多光谱影像相比，色彩的饱和度、亮度、对比度都有不同程度的提高，且地物的纹理细节与边缘特征较为突出，可看出经过治理后沉陷区内植被的颜色变化，能有效提取植物生长高度及密度等信息。GS融合影像包含的信息量较少，色彩的亮度、对比度和地物的纹理细节和边缘特征均较差。Brovey融合影像在色彩上与原始多光谱影像相差较大，根据其原理可知在对 RGB影像进行变换或计算时会丢失一部分的光谱信息，但地物的纹理与边缘特征得到了增强。

图2 几种算法融合前后的QuickBird影像对比

图3 几种算法融合前后的高分二号影像对比

3.3.2 定量评价

通过分析QuickBird和高分二号融合前后的影像参数(表1、表2):

表1 QuickBird原始多光谱与融合后影像的参数

表2 高分二号原始多光谱与融合后影像的参数

QuickBird影像数据在空间信息融入度方面，标准差代表了影像灰度级的分布状况，PCA融合影像的标准差最大(值为127.295)，灰度级分布最分散，包含的信息量最多；信息熵反映了影像信息丰富情况，PCA融合影像的信息熵最大(值为8.125)，表明该融合影像包含的平均信息量最多；平均梯度是衡量影像对地物的细节和纹理信息的表征能力大小的指标，GS融合影像的平均梯度最大(值为23.415)，其空间细节信息表现最清楚，PCA和NNDiffuse融合影像次之，Brovey融合影像的平均梯度最小；相关系数反映的是原始多光谱影像与融合影像之间的相关程度，NNDiffuse和GS融合影像的相关系数较大，分别为0.961，0.956。在光谱信息保真度方面，灰度均值反映了影像的平均亮度，GS融合影像的灰度均值(值为329.293)最接近原始影像，光谱保持最好，与目视结果相一致；PCA和GS融合影像的偏差指数均较大，NNDiffuse融合影像的光谱扭曲度最小，Brovey融合影像的光谱扭曲度最大。

高分二号影像数据在空间信息融入度方面，GS融合影像的标准差(值为712.718)最大，NNDiffuse和PCA融合影像次之，Brovey和HSV融合影像较小；NNDiffuse融合影像的信息熵最大(值为11.258)，信息丰富程度最高，PCA和GS融合影像次之；GS融合影像的平均梯度(值为97.767)较其他融合影像大，表明其对地物的细节和纹理信息的表征能力最大，清晰度最高；而相关系数除GS融合影像较低外，其他3个融合影像都与原始多光谱影像具有高度相关性，PCA融合影像的相关系数最高，达到0.942。在光谱信息保真度方面，PCA融合影像的灰度均值(值为990.059)最接近原始多光谱影像的灰度均值(值为1 046.940)，原始影像信息保留度最好，GS融合影像和NNDiffuse融合影像次之，Brovey融合影像和HSV融合影像最差；PCA融合影像的相对偏差(值为0.124)和光谱扭曲度(值为152.030)最小，与原始多光谱影像最匹配，NNDiffuse融合影像和GS融合影像次之，Brovey和HSV融合影像的原始光谱信息丢失最多。

3.3.3 评价结果

惠农采煤沉陷区治理工程从2004年开始至2015年基本完成。因此，综合前文定性评价和定量评价的结果，对于研究区进行矿山地质环境治理恢复前的2003年QuickBird数据，采用GS融合算法最易于获取因矿山开采引起地质灾害或占用、破坏土地资源的典型地物，如地裂缝、煤矸石、煤堆、矿坑水等的范围、纹理及形态信息(图4)；而对于研究区进行矿山地质环境治理恢复后的2018年高分二号数据，采用PCA融合算法最适合提取复绿植被、观景平台、中心广场等人造景观的几何结构和纹理信息等(图5)。研究成果可推广到其他类似研究区，可为采煤沉陷区地质环境监测提供高质量遥感基础数据。

图4 GS融合处理后的QuickBird遥感影像

图5 PCA融合处理后的高分二号遥感影像

4 结 论

1)通过对4种融合算法得到的融合影像与多光谱影像进行对比分析，QuickBird影像采用GS和PCA融合算法得到的融合影像与原始多光谱影像相比，色彩的亮度、饱和度和对比度有很大程度提高，可准确判读治理前煤矸石的范围与形态、矿区建筑物边界、道路等；高分二号影像的PCA和NNDiffuse算法的融合效果优于Brovey、HSV和GS融合算法，且PCA融合影像较NNDiffuse融合影像的地物几何结构和纹理信息表现得更为明显，易于采煤沉陷区地表环境信息的提取。

2)通过对比分析融合影像的7个定量评价指标值，QuickBird影像数据的GS融合影像不仅在空间信息融入度方面较好，而且具有很强的光谱保真能力。其与原始多光谱影像相关性最高，像元灰度级分散，包含的平均信息量多，空间细节信息表现最清楚；高分二号影像数据在空间信息融入度方面，GS融合影像最优，像元灰度级较分散，影像较清晰。在光谱信息保真度方面，PCA融合影像最优，与原始影像灰度均值最接近，相对偏差和光谱扭曲度最小，较其他融合影像光谱信息丢失程度小。

3)综合融合算法定性和定量评价结果，GS融合算法最适合研究区进行矿山地质环境治理恢复前的QuickBird数据提取，如地裂缝、煤矸石、煤堆、矿坑水等的范围、纹理及形态信息；PCA融合算法最易于提取研究区矿山地质环境治理恢复后的高分二号数据中复绿植被、观景平台、中心广场等人造景观的几何结构和纹理等信息；研究成果可为采煤沉陷区地质环境监测提供高质量遥感基础数据。