8.5.1 方法与流程
8.5.1.1 岩性信息产品生成业务化流程
采用的高光谱岩性信息分类填图方法：先将成像高光谱数据进行辐射定标、大气校正、光谱重建等过程得到光谱反射率数据，然后根据矿物特征吸收峰分布情况对高光谱反射率数据进行波谱降维，对照已有地质图选取其中的典型岩性并结合像元纯净度指数图在高光谱数据区域中建立感兴趣区，确定岩性分类的先验样区，最后用合适的方法进行岩性信息分类填图，并对结果进行优化操作。岩性信息分类填图流程图如图8.20。
图8.20 高光谱岩性信息分类填图流程
8.5.1.2 岩性分类信息提取的高光谱数据预处理
无论是高光谱成像仪还是传统的多光谱传感器，它们所记录的数据都是地面观测目标的反射或辐射能量的光谱辐射绝对值，与地物目标的光谱反射率或光谱辐亮度值是不一致的。因此，辐射定标和光谱重建是地物识别和定量分析不可缺少的环节。只有经过辐射标定、辐射校正和大气校正，剔除由于大气散射、吸收、地形起伏及传感器本身不稳定带来的各种失真，将记录的图像值转换为地面的反照率值，重建像元地面光谱，才能根据光谱特征，有效地识别地物，反演地物成分。
8.5.1.3 端元选择
对于一个地区的未知岩性分类信息提取，一般要进行岩性端元信息的提取。所谓端元，指的是成分单一的岩性像元。只有提取出端元才能进一步对研究区域进行岩性信息的分类分析。目前端元选择的方式概括起来主要有两种：①根据野外波谱测量或从已有的地物波谱信息库中选择端元。通过这种途径选择的端元称为“参考端元”。②直接从待分类岩性的图像上选择端元，然后不断对其修改、调整，确定端元，这种图像上选择的端元区域称为先验样区。
野外获取岩性端元信息一般要经过实地勘察，先选择好样区，然后选择合适的时间进行量测。一般情况下，要求在获取影像的同时进行地物波谱量测，但这种难度很大，实际应用中很少能做到。通过野外测量方式获取的参考端元理论上比较精确，但遥感图像上地物的波谱曲线受到大气、地形和传感器等的影响，这些与野外实地测量的地物波谱曲线存在很大差别，即使对影像进行各种纠正进一步消除这些因素的影响，也不会与野外测量的地物波谱曲线很好地吻合。对于大部分研究区域而言，有与研究区对应的实际地物波谱库的区域很少，也很少投入大量人力、物力进行野外波谱测量。因此在缺乏野外波谱测量数据情况下，从影像本身获取端元是目前获取端元的主要方式。目前，在端元从影像本身的像元获取方式中，除了对遥感影像像元的色调、波谱特征进行目视解译直接分析外，还可借助以下分析方法进行。
基于图像的端元选择是假定图像中存在有基本上仅反映一种岩性或矿物光谱的“纯像元”或“非混合像元”，用数学方法自动或交互地从图像中提取这些“纯像元”作为端元，用同一类“纯像元”的典型光谱或平均像元光谱作为端元光谱。PPI方法由于推出较早，计算简便直观，并且处理高光谱数据中使用最广泛的遥感图像处理系统ENVI中有该功能，因而得到较广泛的应用。但是在进行PPI之前，需先对高光谱数据进行降维处理。我们常用的降维处理手段主要是最大噪声分离（MNF）变换。
（1）最大噪声分离（MNF）变换
成像光谱的光谱分辨率很高，波段多，数据海量，且波段之间的相关性很强，数据冗余度高，需要在端元选取和矿物识别之前对数据作减维，并弱化噪声。最常用的方法是“最大噪声组分变换”。
MNF变换是利用图像的噪声组分矩阵（∑N∑-1）的特征向量对图像进行变换，使按特征值由大到小排序的变换分量所包含的噪声成分逐渐减小，图像质量顺次提高。∑为图像的总协方差矩阵；∑N为图像噪声的协方差矩阵。MNF有两个重要的性质，一是对图像的任何波段作比例扩展，变换结果不变；二是使图像矢量、信息分量和加性噪声分量互相垂直，达到信息分离的目的。乘性噪声可通过对数变换转换为加性噪声，变换后可针对性地对各分量图像进行去噪，或舍弃噪声占优势的分量，以达到成像光谱数据减维和去噪的目的。
（2）像元纯度指数（PPI）分析
基于上述MNF变换，排序低的MNF波段被暂时忽略，仅选择高序次波段进一步处理。PPI设计指定光谱极值像元，对应为混合光谱端元。通过反复投影n维散点图到随机单元矢量来计算PPI。记录每次投影的极值像元，注记每个像元被标定为极值的总次数。PPI图像产生，其中，每个像元的DN值对应像元被记录为极值的次数。这些图像的直方图显示被PPI“击中”（hit）的分布。从直方图中选择阈值，用于选择最纯的像元以保证被分析的像元数最小。这些像元被输入到分离特定光谱端元的交互式可视化算法中。
（3）N维可视化（n-Dimensional Visualization）
由于PPI算法本身确定的并不是最终的端元，而是从图像众多像元中选出包含所有端元像元的较小像元子集。所以，一般要将PPI的处理结果输入到N维可视化（n-Dimensional Visualization）工具中，选择出最终的端元像元。
在N维可视化中，光谱可视为n维散点图中的一个点，n是波段数。对给定的像元，n维空间中组成n值的点的坐标是每个波段的光谱反射率。这些点在n维空间的分布可用于估计光谱端元数和它们的纯光谱特征。在两维空间，如果只有两个端元混合，混合光谱将落入直方图的线中。纯端元将落入混合线的两端；如果三个端元混合，混合像元将落入四面体中；余类推。混合的端元落在纯端元之间，处在纯端元勾画的多面体中。这种混合光谱的凸面几何特征，可用于确定端元光谱数并估算它们的光谱特征，经过PPI选择出的潜在端元光谱输入n维散点图中进行反复旋转以识别出纯端元。根据前面的分析，较好的端元通常会出现在n维散点图的顶点和拐角处，当一系列的端元点被确定后，就可以将其输入到图像中的感兴趣区（ROI），从图像中提取每个感兴趣区平均反射率光谱曲线作为成像光谱矿物填图的候选端元。
鉴于PPI算法得出的结果是从图像众多像元中选出包含所有端元的较小像元子集（图8.21），在已有的地质图中结合PPI结果图可选取尽量纯净的端元区域作为后期监督分类的感兴趣区（ROI，Region of Interest）。
图8.21 选择出的纯净像元的波谱曲线
8.5.2 结果与分析
目前，从是否需要先验样区可将岩性分类技术分为两大类：非监督分类和监督分类。非监督分类是指人们事先对分类过程不施加任何的先验知识，而仅凭数据遥感影像地物的光谱特征的分布规律，即自然聚类的特性，进行“盲目”的分类；其分类的结果只是对不同类别达到了区分，但并不能确定类别的属性；其类别的属性是通过分类结束后目视判读或实地调查确定的。非监督分类也称聚类分析。一般的聚类算法是先选择若干个模式点作为聚类的中心。每一中心代表一个类别，按照某种相似性度量方法（如最小距离方法）将各模式归于各聚类中心所代表的类别，形成初始分类。然后由聚类准则判断初始分类是否合理，如果不合理就修改分类，如此反复迭代运算，直到合理为止。监督分类又称训练场地法，是以建立统计识别函数为理论基础，依据典型样本训练方法进行分类的技术。即根据已知训练区提供的样本，通过选择特征参数，求出特征参数作为决策规则，建立判别函数以对各待分类影像进行的图像分类，是模式识别的一种方法。要求训练区域具有典型性和代表性。判别准则若满足分类精度要求，则此准则成立；反之，需重新建立分类的决策规则，直至满足分类精度要求为止（图8.22）。
图8.22 地质图
对原始数据的211个波段进行目视解译，去掉其中明显的噪声波段得到剩下的178个波段，对这178个波段的数据信息用上述方法进行分类得到所需的结果，然后结合各种岩性地质特征，得知大多数矿物特征吸收峰均位于1300 nm以后的波段中，故在这178个波段中继续裁剪得93个波段，并对裁剪结果用相同的方法及参数重新分类，将得到的结果与第一次分类结果进行对比。
8.5.2.1 非监督分类填图结果
（1）K-Means
预选待分类数据可分为八类，最大迭代次数为2，且设变化阈值为5%，即当每一类像元数变化小于阈值时结束迭代过程，最大允许标准差和最大允许距离误差不输入，即所有像元都参与分类。从而，当达到阈值5% 或迭代达到两次时则分类结束（图8.23 ，图8.24）。
图8.23 所用数据：dts_178，K-Means填图结果
图8.24 所用数据：resize_dts_93，K-Means填图结果
（2）ISODATA
预选待分类数据可分为5～10类，最大迭代次数为2 ，变化阈值为5%，每一类最少含有像元数为500 ，最大允许标准差为10 ，即如果一类的标准差大于10 ，则该类被拆分为两类。类均值间允许最小距离为5，能够被合并成对的最大对数为5，即当类均值间距离小于5 时，这一类就会被合并，而合并后的成对类的最大值为5（图8.25 ，图8.26）。
图8.25 所用数据：dts_178，ISODATA 填图结果
图8.26 所用数据：resize_dts_93，ISODATA填图结果
8.5.2.2 监督分类填图结果
首先，对照地质图和PPI图像选取两种岩性（Cgammabeta：石炭纪黑云母花岗岩，Cdelta：石炭纪闪长岩）的感兴趣区。
（1）平行六面体
设置最大允许标准差（Max stdev from Mean）为1.4（图8.27，图8.28）。
图8.27 所用数据：dts_178，平行六面体填图结果
图8.28 所用数据：resize_dts_93，平行六面体填图结果
（2）最小距离
设置最大标准差为10 ，最大允许距离误差为2500 ，则分类过程中由两者中较小的一个判定像元是否参与分类，若一旦大于任何一个值则该像元不参与分类，归属为无类别（图8.29 ，图8.30）。
图8.29 所用数据dts_178，最小距离填图结果
图8.30 所用数据resize_dts_93，最小距离填图结果
（3）光谱角制图
设置最大允许角度为0.05°，即当像元波谱与终端端元波谱间夹角大于0.05°时，则不参与分类（图8.31 ，图8.32）。
图8.31 所用数据：dts_178，光谱角制图填图结果
图8.32 所用数据：resize_dts_93，光谱角制图填图结果
（4）光谱信息散度
设置最大散度阈值为0.002（图8.33）。
图8.33 所用数据：dts_178，光谱信息散度填图结果
（5）二值编码
设置最小二进制阈值为0.95 ，决定了哪些像元参与分类（图8.34 ，图8.35）。
图8.34 所用数据：dts_178，二进制编码填图结果
图8.35 所用数据：resize_dts_93，二进制编码填图结果
（6）最小距离（Hymap数据）
该方法需要设置两个阈值：Max stdev from Mean；Max Distance Error。通过不同的尝试，针对两种阈值的设定如图8.36所示。
图8.36 局部填图效果
（7）光谱信息散度（Hymap数据）
该方法需要设置阈值：Maximun Divergence Threshod。通过不同的尝试，参数设置和填图效果如图8.37所示。
图8.37 Threshod＝0.005 时的谱信息散度填图结果
8.5.2.3 岩性信息填图结果分析
为了分析高光谱星载模拟数据的岩性填图效果（图8.38），可对比研究区的岩性信息分布图与地质图逐像素进行对比，如果研究区总的点数为N，地质图中与岩性信息分布图的结论相同就认为该点取值1，最后统计结果中1的个数n，这样岩性信息识别精度可定义为
高光谱遥感技术原理及矿产与能源勘查应用
图8.38 岩性提取结果图
笔者发现岩性分类结果与地质图吻合度达到了89%，对出现误差的主要原因分析如下（图8.39）：
1）实验中所用数据是高光谱星载模拟数据，其空间分辨率为30m，随着分辨率的降低，单个像元所对应的地面面积将增大，导致每个像元中包含更多的岩性类型，岩性间的影响性也会增加。一方面，由于混合像元的平均效应，目标岩性在像元中的等效丰度会下降，光谱信息减弱。当像元中目标岩性的等效丰度下降到检出限以下时，岩性将不能被识别，而造成岩性分布区的外围含量较低的地段和含量较低分布区填绘面积的减少，点状集群分布区的漏识别，线状分布区的断续分布。这种效应相当于检出限的下降。另一方面，当像元中目标岩性的等效丰度在检出限以上时，会使面状岩性分布区的范围扩大、空洞的充填、相邻小区的连接，点状集群分布区的成片，线状区域的斑点效应。这两种效应的综合作用结果，空间分辨率的减小则会使强异常区（包括高丰度区和光谱反衬度较高的矿物分布区）范围的扩大、小区的相连、点群的结合，而使异常更加醒目，但其细节特征会因此消失；弱异常区（包括低丰度区和光谱反衬度较低的那些岩性分布区）和小异常区面积会缩小或漏检；线状异常可能会形成断续分布的小斑块而使线状特征和其走向变得不清晰。因此，在填图的细致程度上，分辨率越高的数据填图效果越好，这是造成填图效果差异的最主要原因之一。
2）辐射校正的精确程度。无论是高光谱成像仪还是传统的多光谱传感器，它们所记录的数据都是地面观测目标的反射或辐射能量的光谱辐射绝对值，与地物目标的光谱反射率或光谱辐亮度值是不一致的。因此，辐射校正和光谱重建是地物识别不可缺少的环节。但校正后数据必然丢失一定的信息，故再进行后续处理也会有一定的误差。
3）感兴趣区的选取是进行监督分类的一个重要环节，但由于一般情况下并不能得到纯度足够高的感兴趣区，造成选取的样区含有多种岩性，从而对后期处理造成不可避免的误差。
4）在数据一定的情况下，各种分类方法中参数的选取决定了分类效果，但由于无法遍取各种参数进行尝试，实验过程中仅是进行有限尝试后选择了具有较好的结果的参数，但并不能保证所用参数是最适合的参数。
总之，用高光谱星载模拟数据进行矿物填图是可行的，它可以在一定程度上对岩性的种类及分布进行识别。
图8.39 对比地质图（上图是下图黑框中部分）