宣化盆地高光谱遥感第四纪地层分类提取初步研究*

陈立泽， 徐岳仁

(中国地震局 地震预测研究所，北京100036)

宣化盆地高光谱遥感第四纪地层分类提取初步研究*

陈立泽*， 徐岳仁

(中国地震局 地震预测研究所，北京100036)

基于Hyperion高光谱遥感数据，采用像素纯度指数(PPI)和光谱角制图(SAM)等光谱分析方法对河北宣化盆地第四纪地层的分布情况进行定量提取，并对获取的结果进行实地验证和优化.同时，分析了第四纪地层的分布与宣化盆地南缘断裂带活动之间的关系.

高光谱遥感；宣化盆地；第四纪地层

20世纪80年代兴起的高光谱遥感技术能在特定的光谱范围内以细小的波段宽度获取地物完整而连续的光谱曲线，而且能同时将传感器视域内所有地物特征记录下来，在空间域上形成图像，获取的数据具有图谱合一的特征，使得遥感数据在传统的二维图像上增加了光谱维.由此，遥感应用可以针对光谱开展，通过光谱特征对地物进行定量的区分，使遥感技术的定量化水平提高到了一个全新的阶段.自该技术提出以来，美国、加拿大、澳大利亚等发达国家开展了大量的研究和应用开发，在地质调查、岩石矿物识别、土壤监测及精准农业和森林防护等方面应用广泛.其中，地质是高光谱遥感应用中最为成功的一个领域，国外很多学者在地质填图和矿物识别上开展了大量的工作[1～3]，主要基于矿物的光谱吸收特征参数获取矿物的丰度信息及空间分布，进行岩石矿物的识别和填图.受数据源的控制，我国高光谱遥感起步相对较晚，主要应用在农业和林业方面，在地质和第四纪研究中应用相对较少.但是也有不少学者在开展相关研究，努力推进高光谱遥感技术在我国的应用[4～6].第四纪地层因其形成的环境不同，在成分、结构上有所差异，它们的光谱特征互不相同.但是，这种差异往往非常细微，对它们进行探测，需要很高的光谱分辨率才能识别，只有在高光谱技术日臻成熟的今天才能实现.

本文选取位于北京西北的宣化盆地，利用高光谱遥感数据所获取的地物光谱特征，经过光谱分析处理，对该地区第四纪地层定量分析和提取，并探讨宣化盆地第四纪地层分布与构造活动性之间关系.一方面能加强对宣化盆地构造活动性的了解，另一方面是对高光谱遥感第四纪定量分析的一个有益的尝试.在我国正计划发射高分辨率的高光谱遥感卫星之际，利用高光谱遥感数据在不同领域开展应用研究，有利于发挥我国高光谱遥感卫星的运行效益.

1 宣化盆地概况

宣化盆地位于首都圈西北，鄂尔多斯块体东北缘.构造上处于山西地堑系和张家口-渤海断裂带交汇地带.晋冀蒙盆岭构造区内，是一个NW向的不规则第四纪盆地(图1)，有关盆地的形成时代和盆地内新生代松散沉积物的厚度等资料很少，对这个盆地的第四纪沉积及其形成演化过程研究程度较差.研究区内第四纪地层主要沿基岩山体山前地带以及山间沟谷分布[7].盆地南缘有一条断裂，西起水泉村，向东经龙门堡、崞村南，止于魏家沟附近，长约15 km，走向近东西，正断性质，断裂上盘可见第四纪黄土和砾石层.晚更新世以来宣化盆地南缘断裂有过活动，断层活动对第四纪地层和地貌有一定的影响[8].

2 Hyperion数据预处理

美国NASA于2000年11月21日把EO-1(Earth Observing-1)新型地球观测卫星送入太空，该卫星上搭载了第一个星载民用成像高光谱仪Hyperion.Hyperion数据价格昂贵，而且国内获取渠道有限.不过现在可通过USGS网站(http://edcsns17.cr.usgs.gov/NewEarthExplorer)和国际科学数据服务平台(http://datamirror.csdb.cn)获得用于科学用途的免费数据.Hyperion传感器工作谱段从356～2 577 nm，一共242个波段.空间分辨率30 m，幅宽180 km×7.5 km.本文使用的Hyperion于2004年2月3日获取，为L1级HDF格式的产品，由传感器本身造成的可见光-近红外(VNIR)波段与短波红外(SWIR)波段之间的空间错位已进行过系统校正.但是L1级产品上不正常的像元依然很多.在进行地物定量分类提取前，需要对Hyperion高光谱遥感数据进行未标定及水气影响波段的去除、坏线修复、条纹去除、Smile效应去除、大气校正等一系列的预处理[9，10]，以提高光谱分类提取的精度.

2.1 未标定及水气影响波段的去除

Hyperion数据一共242个波段，VNIR波段70个，波段号从1到70，SWIR波段172个，波段号从71到242.NASA只对VNIR8～57和SWIR77～224一共198个波段做了系统标定.VNIR56～57与SWlR77～78波长重叠，而且前者噪音比后者小，通常保留VNlR56～57波段，得到196个独立标定波段.在标定的波段中，121～127、167～178和224波段的工作波长分别为1 356～1 417 nm、1 820～1 932 nm和2 395 nm，这几个波长范围受水汽的影响较大.其中123～125(1 376，1 386，1 396 nm)三个强水汽吸收波段在用FLAASH模型做大气校正时对云层特别是高空云的估算有重要作用，121、122、126、127、167～178和224波段没有其他用途，而且包含的地面信息也较少，因此在做地物分类提取时直接删除，所以总共保留的波段为179个(表1).由于生成L1级产品时，VNIR和SWIR分别采用了40和80作为扩大因子.将像元值转换为绝对辐射值时，VNIR波段和SWIR波段要分别除以40和80进行还原.

表1 保留的Hyperion波段

2.2 坏线的修复与条纹去除

在Hyperion的L1级产品中，仍然存在少量无数据或数据值非常小的行或者列，这里称之为坏线.坏线的出现没有规律性，需要对保留的179个波段进行逐一检查，并用其相邻行或列像元的平均值取代，进行坏线修复.

Hyperion的L1级图像产品中，除了坏线之外，有些波段上还存在条纹，VNIR波段图像上条纹数目较少，但是很明显，SWIR波段图像上的条纹多，但是相对较轻.条纹的存在严重影响了图像的质量和图像的应用.因此，需要进行条纹去除或修复处理的工作.本文采用全局去条纹的方法，效果明显，见图2，而且对整个图像的像元值改变不是很大，能有效的保持光谱曲线的形状.

2.3 图像的Desmile处理

Hyperion产品存在一种比较特殊的效应.由于Hyperion探测器分光方式和推扫成像模式，导致图像在垂直飞行方向上，像元的波长中心位置向两边偏移，这就是Smile效应.VNIR波段和SWIR波段图像分别由两个不同的探测器获取，它们的Smile效应程度不一样.通常VNIR波段的Smile效应比较明显，而SWIR波段比较小.对于光谱变化平缓的水体和其他地物，Smile效应对光谱的影响不明显，但是当光谱变化剧烈时，Smile效应就会很明显，而且大气校正过程也不能消除这种影响.对Smile效应的消除处理称之为Desmile处理.

本文在大气校正前对Hyperion数据采用移动线性适应与差值的处理方式来降低Smile效应的影响.并通过直接查看762 nm两侧波段的差值图像上是否存在亮度梯度判断纠正的效果.图3是Hyperion VNIR40(752 nm)和VNIR42(772 nm)波段的差值图像，可看到Smile效应去除前，有明显的亮度梯度，去除后的差值图像几乎没有亮度梯度，说明Smile效应去除效果良好.

2.4 大气校正

大气校正是高光谱遥感数据进行地表参数定量分析的前提，其主要目的是减少或消除大气散射和吸收作用等因素对从地面反射到传感器的光谱的影响.本文采用Fast Line-of-sight Atmospheric Analysis of Hypercubes(FLAASH)基于辐射传输方程理论模型进行校正.该方法计算精度高，但是计算过程需要参数多以及计算量大.表2是本文进行FLAASH输入的参数，大部分参数可以从高光谱数据的元数中直接得到或计算得出；少量参数，比如模型的选取，根据成像条件和经验得出.图4是大气校正前后的光谱曲线对比.

表2 FLAASH大气校正输入参数

经过大气纠正后的光谱和野外用光谱仪采集的光谱具有很好的相似性，表明FLAASH大气校正是有效的[11].大气校正后把水汽吸收严重的123～125波段剔除，保留176个波段做进一步的分析处理.从校正后的光谱曲线上可以看出有很多连续锯齿噪音，需要进一步处理.首先采用MNF转换的方法来进行光谱平滑，VNIR和SWIR波段图像由两个不同的传感器获取，噪音表现特征不一样，因此分开进行MNF转换，然后选取特征值高的MNF波段的子集进行MNF反向转换，最后将两个MNF反向转换图像合在一起，生成整体的MNF去噪图像.经过MNF去噪后的光谱曲线还有较多的锯齿，再使用Savitzky-Golay滤波器进行平滑滤波，这是一种在时域内基于局域多项式最小二乘法拟合的滤波方法，在滤除噪声的同时可以确保信号的形状、宽度不变[12]，能在消除噪音的同时比较好的保持光谱形状.滤波处理时，在图像上找一种典型地物，比如水或者雪，通过调整滤波参数，保持图像滤波后的光谱与地物的理论光谱形状一致.

3 第四纪地层分析与提取

3.1 地层非监督分类

在不知道研究区内第四纪地层的种类及分布的情况下，为了快速地了解研究区内地物的种类及分布，可以采用非监督分类对高光谱图像进行分析.经过大气校正后，研究区的Hyperion保留了176个波段，其中有些波段信息冗余而且存在噪声.在这种情况下，需要一种特征降维的方法来减少特征数，减少噪音和冗余，减少计算量.本文采用主成分分析法(Principal component analysis，PCA)对高光谱数据进行降维处理[13]，然后对其进行非监督分类，获取研究区地物信息.

主成分分析法是一种数学变换的方法,它把给定的176个波段认为是一组相关变量，通过线性变换转成另一组不相关的变量，这些新的变量按照方差依次递减的顺序排列，并保持变量的总方差不变.本文选取了前五个主成分分量，其特征值总和已达到99.74%，认为基本已包含了原来176个波段的信息.

非监督分类采用Iterative Selforganizing Data Analysis Techniques Algorithm(IsoData)算法.通过设定初始参数而引入人机对话环节，并使用归并与分裂的机制，当聚类结果某一类中样本数太少，或两个类间的距离太近时，将它们合并为一类，当聚类结果某一类中样本某个特征类内方差太大，将其分为两类.根据初始聚类中心和设定的类别数目等参数迭代，最终得到一个比较理想的分类结果.用IsoData进行分类处理，需要对预期的类数、每一类中允许的最少样本数目、类内各特征分量分布的相对标准差上限、类中心间的最小距离下限、在每次迭代中最多可以进行“合并”操作的次数、允许的最多迭代次数进行初始化.

3.2 地层光谱角制图

光谱角制图(Spectral Angel Mapper,SAM)是一种比较直接的利用光谱特征进行地物分类提取的方法.进行分类提取前，需要知道所有要提取类的光谱特征.这里，我们通过计算像素纯度指数(Pixel Purity Index，PPI)获取所要进行分类的地物组分端元.因为有些噪音的光谱特征也很独特，在PPI计算中也会体现出较高的指数.为了减少噪音的影响，在做PPI运算前，先对数据做最小噪音分离(Minimum Noise Fraction，MNF).进行地层光谱角制图，首先要对高光谱数据进行最小噪音分离，根据特征向量的变化，选取信噪比高而且图像质量好的MNF波段作为PPI提取数据集，再利用PPI指数高的像元作为待分类的地物组分端元，最后通过SAM制图方法，获取研究区各种地物的分布情况.

MNF变化可以看做为两步叠加的主成分变换[14].第一步，基于图像噪声协方差矩阵对噪音进行分离变化，使变换后的数据噪音方差和波段之间的相关性都最小；第二步，对噪声分离后的数据进行标准主成分变换，变换后得到的向量中的各元素互不相关，第一分量集中了大量的信息，而且信噪比也最高，随着维数的增加，影像质量逐渐下降，信噪比也减小.经过MNF变化的图像基于特征值可以分为两部分：一部分特征值较大，图像信息丰富；其余部分特征值较小而且近似，以噪音为主.本文对经过大气校正后的176个Hyperion波段进行最小噪音分离变化后用于PPI指数提取.

进行PPI提取的目的是为了找到高光谱数据中那些光谱很纯或者独特的像元.它们往往代表某一类型具备独特光谱特征的地物，可用作高光谱分类的组分单元.PPI算法把经过MNF变化去噪降维后的图像所有像元重复的投影到一个n维随机单位向量空间里，并使处在边界位置的像元向量出现在随机单位向量边缘的概率最大.记录每次投影中落在单位向量端部的独特像元及每个像元被标记为独特像元的次数.被标记为独特像元的频率越高，认为该像元光谱可区分性越强.像素纯度指数高的像元主要代表可区分的地物类型.

光谱角制图是对地物光谱波形相似性的一种度量.光谱角分类将Hyperion数据的每个像素的光谱都视为波谱空间的一个矢量，波谱空间的维数与数据的波段数一样.在这个波谱空间里，由光照或阴影的作用等因素的影响造成同一类型地物的不同像素上的总光照增加或减小,这个像素对应的矢量长度将增大或减小,但是它的角度方向将保持不变.不同的地物其光谱角度方向是不同的.计算图像上所有像元与前面基于PPI提取出的纯像元之间的光谱矢量之间的夹角，以确定其相似度.二者之间的夹角越小，说明光谱越相似，属于同一类型地物的可能性就越大.利用光谱角进行分类，可以减小由太阳入射角、地形和观测角等因素引起的同物异谱现象.此外，光谱角对于乘性干扰具有良好的抗干扰性，不受照度变化的影响[15].

3.3 地层分类结果比较

对地层非监督分类和光谱角制图结果进行比较，见图6.用IsoData进行非监督分类，在宣化南部地区分出了6类地物.其中，盆地内部主要有4类：分别用黄、暗黄、深绿、蓝色，前三种为第四纪沉积物，最后的蓝色代表水体.南部山区比较杂乱，有绿色、黄色、暗黄色，还有蓝色存在.绿色可能代表山区的基岩.黄色、暗黄色为第四纪沉积物，在山区有一定的分布.可能由于基岩风化，形成现代土壤，或者盆地地区干燥沉积物由于风的原因搬运到了山区，导致山区的基岩分类中夹杂着第四纪沉积物的分类.在山区发现有代表水体的蓝色，说明非监督分类存在一定的错误.对此，本文分析认为山区有阴影，阴影地区的信息量相对较少，而水体的反射率低，总体的信息量也少，导致基于主成分分析后的非监督分类把信息量少的这些像元归到了同一类.

光谱角制图的结果上分的类与Isodata非监督分类的结果大致相同，但是可以看出光谱角制图分类结果更加精细，地物的分布也更为合理.盆地内主要有黄色和暗黄色两类第四纪沉积物.蓝色对应水体，需要指出的是在蓝色水体的北部还有一种粉红的分类，这是由于该地区水浅，水底泥沙的光谱参与了分类造成的.在南部山区，紫色对应基岩，深绿色对应黄土，黄色和暗黄色对应第四纪沉积物.山区阴影位置，在利用光谱角制图中没有出现像非监督分类一样把阴影当做水体的情况.阴影区的反射率尽管比较低，但光谱角分类能消除或减弱因太阳入射角、地形、坡向等因素对分类结果的影响，基于光谱特征分类能获取更好的结果.

4 结 论

为了检验基于Hyperion分类结果，对宣化盆地南缘断裂周边地区进行了实地考察验证，考察点的具体分布见图7.考察结果显示，宣化盆地出露的主要第四纪地层为晚更新统(Qp3)，主要分布于山麓地带，沿坡麓堆积，以风成马兰黄土为主，局部夹杂薄层地表片流堆积或洪积砂砾层、砂土层，厚度变化较大，对应图中晚更新世黄土台地.全新统(Qh)沿基岩山前坡麓地带、盆地以及沟谷分布，成因类型以坡洪积、冲洪积和冲积为主.其中坡洪积主要沿沟口呈冲洪积扇型分布，见图7，这些冲洪积扇已基本相互连结而形成山前缓倾斜平原，以砾块、砂、砂土和亚砂土为主，靠近坡麓黄土台地，多叠置有次生含砾黄土；冲洪积层沿接近盆地中心地带的冲洪积平原分布，对应图中的平原堆积；冲积层主要沿沟谷及其阶地分布，对应图中的河流阶地和现代河流沉积.

高光谱遥感分类比较准确的分出了晚更新统黄土，全新统(Qh)坡洪积、冲洪积分布情况，而且断层的线性陡坎在高光谱图像上也有明显的体现.说明高光谱遥感在第四纪研究和活动构造研究中能发挥重要的作用.利用高光谱遥感分类和解译结果，可以对研究区内宣化南缘断裂的构造活动性做出初步的评价.

宣化盆地南缘断裂南侧为基岩山区，与盆地相对高差500～700 m，出露地层主要为侏罗系火山岩、沉火山岩.从基岩山麓向盆地过渡，往往顺坡堆积黄土而发育黄土斜坡或黄土台地.顺坡堆积的黄土与侏罗系火山岩之间呈断层接触，由于后期经历了断层作用、黄土加积、坡积作用以及地表侵蚀作用的共同改造，线性地貌陡坎位于黄土台地前缘，高15～20 m.断层全新世活动，假设陡坎高度代表晚更新世以来断层累积位错量，因此估算断层的平均滑动速率约为0.2 mm/a，未来有发生地表破坏性地震的可能性.

[1] GERSMAN R, BEN-DOR E, BEYTH M, et al. Mapping of hydrothermally altered rocks by the EO-1 hyperion sensor, northern danakil depression,eritrea[J]. International Journal of Remote Sensing, 2008,29(13):3 911-3 936.

[2] RESMINI R G, KAPPUS M E, AlDRICH W S, et al. Mineral mapping with hyperspectral digital imagery collection experiment (HYDICE) sensor data at Cuprite, Nevada, USA[J]. International Journal of Remote Sensing, 1997, 18(7):1 553-1 570.

[3] VANE G, GOETZ A F H. Terrestrial imaging spectroscopy: current status, future trends[J]. Remote Sensing of Environment, 1993, 44:117-126.

[4] 甘甫平,王润生,杨苏明.西藏Hyperion数据蚀变矿物识别初步研究[J].国土资源遥感, 2002,4:44-50.

[5] 王青华,王润生,郭小方.高光谱遥感技术在岩石识别中的应用[J].国土资源遥感, 2000, 4:39-43.

[6] 张宗贵,王润生,郭小方,等.基于地物光谱特征的成像光谱遥感矿物识别方法[J]. 地学前缘, 2003, 10(2):437-205.

[7] 徐锡伟,吴卫民,张先康,等.首都圈地区地壳最新构造变动与地震[M].北京：科学出版社.2002.

[8] 冉勇康,陈立春,徐锡伟.北京西北活动构造定量资料与未来强震地点的讨论[J].地震学报, 2001, 23(5):502-513.

[9] BARRY P. EO-1/Hyperion Science Data User's Guide[R]. Redondo Beach, CA: TRW Space, Defense & Inform Syst, 2001.

[10] 谭炳香,李增元,陈尔学,等. EO-1 Hyperion高光谱数据的预处理[J]. 遥感信息,2005, 6:36-41.

[11] 袁金国,牛铮,王锡平.基于FLAASH的Hyperion高光谱影像大气校正[J].光谱学与光谱分析, 2009, 29(5):1 181-1 185.

[12] SAVITZKY A, GOLAY M J E. Smoothing and differentiation of data by simplified least squares procedures[J]. Analytical Chemistry, 1964, 36(8):1 627-1 639.

[13] RICHARDS J A. Remote sensing digital image analysis: an introduction[M]. Berlin：Springer-Verlag，1999:240.

[14] GREEN A A, BERMAN M, SWITZER P.et al. A transformation for ordering multispectral data in terms of image quality with implications for noise removal[J].IEEE Transactions on Geoscience and Remote Sensing, 1988, 26(1):65-74.

[15] KRUSE F A, LEFKOFF A B, BOARDMAN J B, et al. The spectral image processing system (sips)-interactive visualization and analysis of imaging spectrometer data[J]. Remote Sensing of the Environment, 1993, 44:145-163.

责任编辑：龙顺潮

Preliminary Classification of Quaternary Strata at the Xuanhua Basin Based on Hyperspectral Remote Sensed Data

CHENLi-ze*,XUYue-ren

(Institute of Earthquake Science, China Earthquake Administration， Beijing 100036 China)

Hyperspectral remote sensing technique has been developed rapidly since 80s of the 20th century,because it has the advantages of high spectral resolution and the ability to present imagery and spectral profile together. In this paper, the Quaternary strata at the Xuanhua Basin in the Hebei Province was extracted from the Hyperion data using spectral analysis methods such as Pixel Purity Index(PPI), Spectral Angel Mapper(SAM) and so on. The extracted results were verified with field work. Based on the improved results, the relationship between Quaternary strata and activities of the Southern Margin Fault at the Xuanhua Basin was discussed.

hyperspectral remote sensing;Xuanhua bains;quaternary strata

2014-08-11

中国地震局地震预测研究所基本科研业务费专项(2010IES0203)；科技部国际科技合作项目(2012DFR20440)

陈立泽(1979— ),男，湖南 隆回人，副研究员.E-mail：chenlize@hotmail.com

TP75

A

1000-5900(2015)01-0081-06