深度卷积网络在航空高光谱岩性识别中的应用<br/>——以塔木素铀矿床北部地区为例

深度卷积网络在航空高光谱岩性识别中的应用
——以塔木素铀矿床北部地区为例

2024-03-14张川易敏童勤龙叶发旺徐清俊李泊凇

世界核地质科学 2024年1期

张川，易敏，2，童勤龙，叶发旺，徐清俊，李泊凇

1 核工业北京地质研究院遥感信息与图像分析技术国家级重点实验室，北京 100029

2 北京智信遥感地理信息技术有限公司，北京 100032

高光谱遥感是利用窄且连续的光谱通道对地物进行连续遥感成像的技术，其本质上改变了传统的遥感监测方式，能够直接检测在宽波段遥感影像中无法检测的物质组分。矿物识别是迄今为止高光谱遥感最为成功的应用方向之一，特别是对于指示找矿的蚀变矿物，如含铁矿物、白云母族矿物、绿泥石族矿物和碳酸盐矿物等［1-6］。在区域地质调查中，岩性识别亦是一项基础而又十分重要的工作，对于研究一个地区的成矿潜力和矿产远景规划具有重要意义。然而，当前高光谱遥感的岩性识别远不及矿物识别应用广泛，主要有以下方面的原因：1）岩石是多种矿物的自然混合物，然而，岩矿光谱混合机理目前仍不明确；2）岩石的光谱主要体现的是主导矿物的光谱特征或光谱组合特征，故传统岩石类型分类主要依据不同端元矿物及其含量差异的三角图解进行，然而长石、石英等许多造岩矿物在高光谱遥感常用的0.4～2.5µm 波长范围内缺少明显的诊断性光谱特征；3）尽管热红外高光谱能够识别造岩矿物，但数据噪音大且数据源较少。上述原因导致高光谱遥感在岩性识别中的应用受限。

近年来，机器学习越来越多地引入到遥感地质勘查中［7-9］。这类方法是数据驱动的，能够自动学习遥感数据与所需特征之间的关系。在机器学习领域，深度学习是目前最受关注的技术，它是一种通过深度神经网络（Deep Neural Networks，DNNs）进行特征学习的方法，一些学者将深度学习引入到高光谱遥感图像分类中，获得了超过90 %甚至95 %以上的高分类精度，展示了深度学习在高光谱图像特征挖掘方面的巨大潜力［10-11］。卷积神经网络（Convolutional Neural Networks，CNNs）是深度学习的代表算法之一，它是一类包含卷积计算和深层结构的前馈神经网络，由交替的卷积层和池化层组成，在网络的末端通常使用几个全连接层作为分类器，十分适合处理遥感图像［12］。尤其是针对传统方法对高光谱遥感影像“图谱合一”特征挖掘的不足，近年来发展了图（图像）-谱（光谱）联合的CNN 方法，已在多个公开的高光谱影像数据集上获得了较好的地物分类识别效果［13-15］。然而这些公开的高光谱影像数据集中的地物类型大多为均匀稳定的农田、林地、水体以及人工地物，它们与地质环境存在差异，是否适用于岩性识别还需要进一步评估。

因此，本文将引入多种深度CNN 架构，将其应用于高光谱遥感影像的岩性识别，探讨其应用效果，为高光谱遥感的岩性识别应用提供参考。试验区为巴丹吉林盆地东部，塔木素铀矿床的北部地区，近年在这一地区获取了大面积的航空高光谱遥感影像数据，是良好的数据源。

1 试验区与数据源

1.1 试验区概况

试验区位于巴丹吉林盆地东部盆山结合带，属于巴丹吉林盆地与巴音戈壁盆地相接地带，巴音戈壁盆地西部的宗乃山—沙拉扎山构造岩浆带中段。大地构造上属于塔里木板块、哈萨克斯坦板块、西伯利亚板块和华北板块四大板块在古生代发生陆-陆碰撞的结合部位，具有复杂多变的地质构造背景［16］。试验区主体位于塔木素特大型砂岩铀矿的北部，属于塔木素铀矿床的蚀源区宗乃山隆起（图1a）。受区域性EW 向走滑断裂活动影响，形成一系列NE、NW 向次级断裂。区内岩浆多期次侵入，大面积出露晚古生代侵入岩，包括华力西期和印支期的花岗岩、华力西期的花岗闪长岩和石英闪长岩等（图1b）。试验区北部和东南部除大面积第四系盖层外，还发育下白垩统和中下侏罗统，前者主要为河流相杂色碎屑岩与沼泽相暗色泥质岩，后者下部为砾岩、砂砾岩、粉砂岩夹煤线，上部为火山角砾岩、流纹岩、英安质晶屑岩晶粗凝灰岩和粉砂质凝灰岩。区内少量石炭系地层发育，主要分布在试验区的北部。受岩浆和热液活动的影响，区内发育不同程度的蚀变，主要有绢云母化、高岭土化和硅化等。

图1 试验区位置和岩性分布图Fig. 1 Location and lithology distribution map of test area

1.2 数据源及预处理

本研究数据源为加拿大ITRES 公司生产的航空高光谱成像仪SASI（Shortwave-infrared Airborne Spectrographic Imager），该仪器为短波红外光谱区域的推扫式成像高光谱传感器。短波红外是地质遥感常用的谱段区域，光谱范围涵盖了多种热液蚀变矿物的吸收光谱。SASI 对其垂直轨道方向的600 个成像像元同步获取在950～2 450 nm 波长范围内具有连续光谱覆盖的100 个光谱波段。飞行相对高度约为2 000 m，获取影像的空间分辨率达到2.5 m，获取的SASI 数据的技术指标如表1 所示。

对于飞行范围较大的航空高光谱影像采集来说，受飞行时间、外界条件等变化因素影响，获取不同航带的辐射信号具有一定的差异，造成航带间辐射亮度不均一，二向反射分布函数（Bidirectional Reflectance Distribution Function， BRDF）校正能消除遥感影像中由于太阳高度角、观测角与真实角度偏差引起的色差不一致，但由于机载遥感图像相对星载遥感图像的分辨率高、视场角大，对光照、环境等因素更为敏感，造成目前机载环境的遥感图像BRDF 校正效果不佳，如图2 所示。虽然还可以通过一些其他预处理减弱这些影响，但亦可能影响到地物原始的辐射特征。因此，为避免预处理过程中引入其他噪声，对试验区的航空高光谱影像仅进行了系统级辐射定标处理、大气校正和反射光谱重建和BRDF 校正等。大气校正统一采用6S 辐射传输模型，通过计算和反演飞行时的大气参数，基于大气辐射传输方程进行大气校正。通过对比校正后的影像像元光谱和野外实测蚀变岩光谱验证大气校正的效果。如图3 所示，为试验区绢云母化蚀变岩所在位置的光谱曲线对比图，绢云母的波谱形态和Al-OH 诊断性光谱特征基本一致，特征吸收峰波长位置偏移仅1 nm，表明光谱重建效果较好。

图3 绢云母化蚀变岩光谱曲线图Fig. 3 Spectral curve of sericitization altered rock

2 试验方法

本研究试验流程包括航空高光谱数据预处理、样本制作、模型构建和模型应用等，如图4 所示。其中，航空高光谱数据预处理在第1.2节中已经叙述。经过预处理后，SASI 高光谱图像的DN 值被转换为反射率值。最后，对各波段影像质量进行检查，剔除1 400、1 900 和2 400 nm 附近受水汽吸收影响的异常值。剔除的波段包括1 332～1 482 nm 的11 个波段、1 797～2 052 nm 的18 个波段和2 382～2 442 nm的5 个波段，剩余66 个波段。样本制作主要基于试验区地质图结合野外调查进行，并按照深度学习模型训练的要求进行比例划分为训练集、验证集和测试集。模型部分针对高光谱影像“图谱合一”的特点，重点设计了1 种基于光谱特征和2 种基于图-谱联合特征的CNN 模型进行对比研究，分别是一维卷积神经网络（One-dimensional Convolutional Neural Network，

图4 试验区航空高光谱数据试验流程图Fig. 4 Test flow chart of aerial hyperspectral data in test area

1D CNN）、一维+二维卷积神经网络（Onedimensional and Two-dimensional Convolutional Neural Network， 1D+2D CNN）和三维卷积神经网络（Three-dimensional Convolutional Neural Network， 3D CNN）。最后，对结果进行评价分析。

2.1 样本制作

深度学习任务需要合适的样本数据来挖掘数据中的特征信息，并抽象出高级信息映射到输出。前人经验表明：深度学习技术作为一个黑箱运行，通过大量的输入和给定的结果提供稳定的处理，训练样本数量越多，黑匣子的工作性能越稳定［19］。从理论上讲，真实数据是最可靠的。然而，深度学习对样本数量的需求较大，并且地面数据与航空高光谱数据存在显著的尺度差异。因此，本文利用影像数据来构建岩性样本数据集。

根据试验区地质图资料及野外调查，在影像上初步圈定8 类，分别是：印支期花岗岩、华力西晚期花岗岩、华力西晚期花岗闪长岩、华力西中期石英闪长岩、石炭系碎屑岩、中下侏罗统火山凝灰岩、第四系沉积物和绢云母化蚀变岩。为评价模型对不同时期的同类岩性的区分能力，花岗岩样本分为印支期和华力西晚期的2 类花岗岩。对初步圈定的样本区域，均进行了野外调查验证，以保障样本的可靠性。另外野外调查发现，下白垩统河流相碎屑岩与第四系沉积物特征基本一致，从影像上看亦如此，故将其归入第四系类。地质图中没有的蚀变岩类，通过影像分析结合野外调查发现，图2影像中偏黄绿色的为蚀变岩（以绢云母化为主），基于此对蚀变岩进行初步圈定。采用影像光谱可视化工具，对每类岩性的像元光谱进行查验，剔除同类样本中与其他光谱差异较大的点，同步计算各类样本的平均光谱，如图5 所示，以构建光谱变化相对平稳的样本数据。从各类岩性平均光谱反映的特征来看，除绢云母化蚀变岩具有与其他岩性差异较大的光谱特征以外，其他岩性的光谱特征较为相似。因此从理论上来说，为了更好地区分这些岩性，不能仅从光谱特征挖掘的角度来构建模型，而应该从图像-光谱结合的角度来构建模型。为提高模型对每个类别的公平性和性能，从每类样本随机抽样2 万个点作为数据集样本，再按照7：2：1 的比例划分为训练集、验证集和测试集，分别用于模型训练、模型调参和模型测试。

图5 每类岩性样本的平均光谱曲线图Fig. 5 The average spectral curve of each type of lithology sample

由于3 种CNN 模型的结构不同，在数据的输入方面存在差异，1D CNN 输入为像素光谱，1D+2D CNN 输入为像素光谱和图像的空间信息，3D CNN 输入为高光谱数据立方体。因此，在样本数据集制作方面，1D CNN 在高光谱数据进行归一化之后，选取光谱曲线作为输入；1D+2D CNN 在1D CNN 的样本数据集基础上，增加了每个样本点周围的空间信息，以归一化的高光谱数据进行主成分分析后的第一主成分表示，范围依据经验设为像素周围30×30［20］；3D CNN 是以样本点为中心，取30×30的立方体。

2.2 模型构建

深度神经网络模型基于目前流行的PyTorch 深度学习框架进行设计、训练与构建。本研究重点测试和对比基于光谱特征的1D CNN、基于图-谱特征的1D+2D CNN 和基于图-谱特征的3D CNN 的岩性识别效果。

1D CNN 通过多个交替的一维卷积层和池化层提取特征，最后使用全连接层进行分类。1D CNN 将像素的光谱曲线作为一维向量输入，利用相应的一维卷积结构来提取高光谱影像的光谱信息特征，结构如图6 所示，卷积层采用的激活函数为ReLU，与其他激活函数如sigmoid、tanh 函数相比，ReLU 函数具有良好的稀疏性和较小的计算量，既避免了梯度消失，又加快了神经网络训练的速度，是目前深度学习中主流的激活函数［21］。

图6 基于光谱特征提取的1D CNN 结构图Fig.6 The 1D CNN structure based on spectral feature extraction

1D+2D CNN 通过增加一系列交替的2D 卷积层和池化层来提取图像的空间特征，并与1D卷积层提取特征进行融合，最后，使用全连接层进行分类。用于提取图-谱联合特征的1D+2D CNN 的结构如图7 所示，在光谱分支上，采用了上述相同的一维卷积模型的结构对光谱特征进行提取；空间信息分支上利用二维卷积层、激活函数以及池化层的堆叠提取图像特征。卷积层采用的激活函数为ReLU，池化层为常用的最大池化。提取获得图像特征和光谱特征后，将图像特征重构成一维张量，再采用Concate 融合法将图像特征与光谱特征进行连接，并使用全连接层对联合图-谱特征进行分类。

图7 基于图-谱特征提取的1D+2D CNN 结构图Fig. 7 The 1D and 2D CNN structure based on image and spectral feature extraction

3D CNN 在结构的组成上，分为基本的三维卷积块、光谱残差块、空间残差块和分类输出层，结构如图8 所示。三维卷积块由三维卷积层、激活函数和归一化层组成，如图9a 所示。光谱残差块和空间残差块在结构上类似，由两个三维卷积块堆叠而成，输入X 经过两个三维卷积块后与输入相加，再通过一个激活函数，如图9b 所示。分类输出层包含一个池化层与全连接层。

图8 基于图-谱特征提取的3D CNN 结构图Fig. 8 The 3D CNN structure based on image and spectral feature extraction

图9 三维卷积块和残差块结构示意图Fig. 9 The structure of 3D convolution block and residual block

光谱残差块和空间残差块用同一个结构，主要区别在于核大小。光谱残差块中的三维卷积层通过（1，1，7）的核进行光谱维度的信息提取，保持空间大小不变。空间残差块中的三维卷积层通过（3，3，1）的核进行空间维度的信息提取，保持光谱维度不变。整体结构上由浅入深，逐步增强光谱与空间的特征表示。两部分特征学习模块通过残差连接构建深层网络，提高了模型的表示能力。

上述神经网络均使用当前流行的深度学习框架PyTorch 进行搭建和训练，损失函数均使用CNNs 中最常用的交叉熵损失函数进行计算［22］。模型训练使用SGD 优化器，batch size 设置为256。为了提高模型训练的稳定性和效率，在初始训练阶段采用了学习率warmup 策略［23］。具体而言，在训练的初始几个周期（epoch），学习率从一个较小的值逐渐增加到预定的最大值。这一策略的引入旨在减少模型训练初期由于学习率过高而导致的不稳定性。学习率的初始值设定为Lrstart，在预定的Nwarmup周期内，学习率线性增加至Lrmax。因此，第t个周期的学习率可由公式（1）得出：

在完成warmup 阶段后，为了实现更加精细的模型优化，还采用了学习率的线性递减策略［23］。具体来说，从第Nwarmup+1 个周期开始，学习率从Lrmax开始逐渐减小，直至训练结束。这种递减策略有助于在训练的后期训练阶段更细致地调整模型参数，从而避免过拟合，并提高模型的泛化能力。根据公式（2）得出从第t个周期的学习率Lrt：

式（2）中：Ntotal—总训练周期；Lrend—结束时学习率。在本研究中，Ntotal—200；Nwarmup—10；Lrstart—0.000 01；Lrmax—0.000 1；Lrend—0.000 000 1。

模型训练获得收敛之后，使用测试集评估模型的泛化性能，完成整个深度卷积神经网络模型构建过程。最后，采用多个评价指标对模型的测试精度进行评估。本文采用的指标包括总体精度、生产者精度、平均精度和Kappa系数，它们都是遥感分类模型的重要评价指标［24-25］。表2展示了不同的CNN模型在试验区岩性样本集上的测试精度评价指标的计算结果。三个CNN模型的总体精度和平均精度均在80 %以上，其中，3D CNN 模型的总体精度和平均精度最高，均大于90 %，1D+2D CNN的精度优于1D CNN，Kappa 系数亦表现为3D CNN＞1D+2D CNN＞ 1D CNN，表明联合图-谱特征的提取方法优于仅考虑光谱特征的提取方法。

表2 卷积神经网络模型测试精度表Table 2 Test accuracy of the CNN models

利用混淆矩阵方法计算各个样本类别的生产者精度，各CNN 模型的混淆矩阵如图10 所示。混淆矩阵中的每一列表示一个预测类别，每一列之和表示该类别预测的数据点总数。每一行代表真实的数据类别，每一行的总和代表该类别中数据实例的数量。基于混淆矩阵，可以计算出每个类别的正确和错误分类比例。表3 是根据混淆矩阵计算的各CNN 模型在测试集上各个类别的生产者精度。

表3 卷积神经网络模型生产者精度表Table 3 Producer’s accuracy of the CNN models

图10 模型测试混淆矩阵Fig. 10 The confusion matrix of the models for the test set

从模型在测试集上的表现来看，联合图-谱特征的1D+2D CNN 和3D CNN 在各类别的精度上均优于仅考虑光谱特征的1D CNN。3D CNN在中下侏罗统火山凝灰岩、华力西晚期花岗闪长岩、华力西中期石英闪长岩、石炭系碎屑岩和华力西晚期花岗岩的精度上优于1D+2D CNN；1D+2D CNN 在印支期花岗岩、绢云母化蚀变岩、第四系沉积物的精度上优于3D CNN。可见引入高光谱图像的空间特征对于提高岩性识别模型精度是有意义的，但具体以哪种方式进行图-谱联合还有待进一步的评价分析。

3 应用结果与评价

3.1 模型应用结果

作为构建CNN 模型的最终目的，为了更好地评估其在高光谱遥感岩性识别中的实际效果，将训练好的各个CNN 模型用于试验区SASI航空高光谱影像进行岩性识别，并将识别结果通过模型进行输出。最后，通过ArcGIS 软件平台对输出结果进行分类和可视化制图，实现岩性识别与分类。图11 展示了分别使用1D CNN、1D+2D CNN 和3D CNN 模型对试验区SASI 航空高光谱影像进行岩性分类识别的结果图及与参考地质图的对比，各图幅的地理范围与图2 完全一致。

图11 试验区岩性识别结果与参考地质图的对比图Fig. 11 Comparison of identified lithology results to the reference geological map in the test area

3.2 结果评价

根据各深度CNN 模型在试验区SASI 航空高光谱影像上的岩性识别结果与参考地质图的对比情况，总体上看，1D CNN 模型和1D+2D CNN 模型的岩性识别结果（图11a、b）有较为明显的、与航空高光谱飞行方向一致的EW 向条带效应，使得试验区的岩性识别成图效果受到显著的影响，而3D CNN 模型的岩性识别结果（图11c）的条带效应明显减弱了许多。由于航空高光谱影像采集受飞行时间、外界条件等变化因素影响，获取不同航带的辐射信号差异造成了航带间辐射亮度的不均一属于正常现象，也是航空影像获取中较为常见的，这导致航带之间色差成为多航带大区域高光谱岩性填图的技术难点［26］。因此，鉴于试验结果中的3D CNN 模型在克服条带效应方面的显著效果，在航空影像的处理和信息识别方面具有更好的应用前景。

在岩性的识别方面，相对而言，3D CNN 模型的岩性识别结果与参考地质图的匹配程度最佳，具有最好的视觉效果；而1D CNN 模型的识别结果与参考地质图的匹配程度最差。对不同岩性来说，首先，对绢云母化蚀变岩的识别，三种模型识别的范围与图2 中黄绿色岩石发育区域基本一致，这是因为绢云母化蚀变岩具有较为明显的诊断性光谱特征（图3），因此，无论是基于光谱的CNN 模型和基于图-谱联合的CNN 模型，均能够较好地学习并提取绢云母化蚀变岩的特征。然而，对于参考地质图中的其他岩性，它们在1.0～2.5µm 波长范围内的诊断性特征不明确，且光谱较为相似（图5），因此在不同模型结果中的差异体现得较为明显。试验区西南部是岩性最复杂的区域，涉及印支期花岗岩、华力西晚期花岗闪长岩、华力西中期石英闪长岩等多种岩性，如图11 所示，1D CNN 模型和1D+2D CNN 模型将印支期花岗岩和华力西晚期花岗闪长岩大量混淆，华力西晚期花岗闪长岩和华力西中期石英闪长岩的混淆程度亦较大，相对而言，3D CNN 模型尽管有所混淆，但混淆程度较小，与参考地质图的岩性分布大体上一致。在试验区北部，1D CNN 模型和1D+2D CNN 模型将部分第四系沉积物识别为石炭系碎屑岩，尤其是1D CNN 大量错分，而3D CNN 对于该地段的第四系沉积物和石炭系碎屑岩的区分效果较好。在试验区东南部，对于大面积第四系沉积物中的中下侏罗统火山凝灰岩的识别，三种模型均有所反映，受条带效应的影响，1D CNN 模型和1D+2D CNN 模型将部分第四系沉积物识别为火山凝灰岩，1D CNN模型错分相对更为严重，对东南角的一片火山凝灰岩，各模型均未识别。在试验区东部，涉及华力西晚期花岗岩、华力西晚期花岗闪长岩的区分，从识别结果来看，只有3D CNN 模型获得了较好的区分效果。

因此，3D CNN 模型无论是对试验区航空高光谱影像条带效应的克服方面，还是对各类岩性的区分识别方面，均获得了最佳的应用效果，这与各模型的测试精度计算结果具有较好的一致性。相比之下，1D CNN 模型对各类岩性的混淆度还要甚于1D+2D CNN 模型，可见，在高光谱遥感常用的反射光谱段，仅依靠光谱特征进行岩性的区分误差较大，图-谱结合是高光谱遥感岩性识别的正确思路。