结合SAR与光学遥感的疑似稀土矿非法开采情况识别方法

乐颖，夏元平，刘媛媛，钱文龙

(1.东华理工大学 测绘工程学院，南昌 330013；2.福州市勘测院，福州 350108)

0 引言

作为我国独特的一种稀土矿，离子吸附型稀土矿具有储量多、分布区域广、种类丰富、用途广等特点[1]。根据省第三环境保护督察组的反馈可知，赣州市证外稀土矿山的总数大约是证内稀土矿总数的2.5倍，分别占据901座与355座，需要治理的矿山数量达到了522座，同时废弃矿山的数量占据了很大比例。因此，针对稀土矿山的违法开采现象，如何实时准确监测开采是目前的研究热点。

合成孔径雷达差分干涉测量(differential interferometric synthetic aperture radar，DInSAR)是地表形变监测、灾害监测、变化监测等诸多应用领域重要信息的获取手段[2-4]。国内外部分学者对雷达变化监测进行了相关研究[5-6]，Kamila等[7]运用DInSAR和时序短基线集差分干涉测量(smallbaseline subset-insar，SBAS-InSAR)技术反演研究区域完整的形变趋势，并用水准数据对结果进行了评价；Xia等[8]提出了一种基于概率原理的方法，利用积分法和合成孔径雷达干涉测量(interferometric synthetic aperture radar，InSAR)技术，成功反演了地下采空区的位置；马威等[9]对比分析了基于InSAR相干系数图与DInSAR对地表信息的提取效果，获得连续变化开采区39处；何曙光等[10]提出了结合InSAR相干系数的差值方法对城市进行变化监测，并证明了该方法提取实验区变化信息具有较高的准确性。同时，国内很多学者通过影像分类或建立明确的解译标志，从而研究分析矿产资源开采变化情况[11-14]。邵艳坡等[15]结合纹理与光谱信息进行分类，结果表明,该分类方法用于定南县稀土矿的总体精度为90.33%。岳建伟等[16]叠加分析遥感与矢量数据，能够自动提取矿区违法开采信息。然而，高分辨率影像一般为光学影像，既受气候环境的限制，又受观测时间的影响，在夜间无法监测。雷达数据具有全天时全天候的优势，不仅提供了反映地物特性的幅度信息和相位信息，还可以获取稀土矿开采引起的地表形变信息。高分辨率影像和雷达数据具有较高的互补性，将二者结合起来进行变化监测具有广泛的前景。

因此，本文结合雷达影像的强度信息和光学影像的光谱信息进行变化监测，并借助已有矿权信息在光学影像上对疑似非法开采区域进行圈定，最后根据地表形变信息和相干性系数图动态监测疑似非法开采区域在2015—2020年间的开采变化情况。该研究揭示了在稀土矿山开采监测中将雷达影像和光学影像相结合进行非法开采识别的能力，同时能够给矿山管理部门治理矿山周边环境提供决策依据。

1 研究方法

1.1 基于相干性系数的 DInSAR 地表形变信息提取方法

利用干涉相干性分析可以提取地表覆盖的变化信息，借助二次差分方法从干涉相位图中去除地形及其他因素的影响，从而达到提取形变信息的目的，该方法被称为DInSAR技术[17-18]。

在DInSAR处理过程中考虑相干系数可以提高结果的准确性[19]，相干系数γ能够直接反映出影像之间的相干性指标[20]。

1.2 基于单波段比值法的遥感变化监测

图像比值法是遥感图像处理过程中最常用的方法之一，就是将不同波段或波段组合所对应的像元灰度值进行比值运算，是对两景影像做除法的变化监测方法。其基本思想是：通过比值运算计算待检测影像中每个像元在不同时相亮度值的比值，然后根据结果绘制出比值图像。检测原理如式(1)所示[21]。

(1)

式中：Ratio表示所需的影像变化监测结果；pwr1和pwr2分别是两个时相影像。若待影像中地类在两个时相未发生变化，其灰度值往往大致相同；若地物类别发生了改变，其灰度值也随之发生变化。因此，通过分析比值影像上地类变化区域的灰度值与背景值之间的差异，然后对比值结果进行一定的处理，最终可获得研究区域内地物发生变化的信息。

1.3 结合光学和雷达影像的稀土矿非法开采识别方法

根据以上理论，本文提出一种结合SAR与光学遥感的疑似稀土矿非法开采情况识别方法，具体流程如图1所示。根据雷达影像的强度信息和光学影像的光谱信息，将赣州市定南县作为研究区域，对离子型吸附型稀土矿的开采情况进行变化监测。本文变化监测结果可通过虚警率(false alarm，FA)、总体错误率(overall error，TE)、总体准确率(overall accuracy，OA)和Kappa系数进行精度评价[22]。具体操作步骤如下。

步骤1：分别对SAR和光学遥感影像进行数据预处理。雷达影像经过影像配准、滤波、地理编码和辐射定标可得到强度图，再利用基于相干性系数的 DInSAR 地表形变信息提取方法得到相干性系数图和形变干涉图。光学影像需经过大气校正、重采样、波段合成、影像镶嵌和裁剪等数据预处理，并对其进行真彩色图像增强。

步骤2：对同一年份的两景影像，根据雷达影像的强度信息和光学影像的光谱信息，分别进行基于单波段比值法的遥感变化监测，并对差异图进行真彩色合成。

步骤3：依据变化标准判断出可能变化区域。若在两个变化监测结果中都未变化，该区域没变化；当在雷达变化监测结果图上变化而光学变化监测结果影像上不变，则该区域没变化；只有当雷达变化监测结果图变化，光学变化监测结果影像也变化才能说明该区域发生了变化。

步骤4：识别出开采区域。通过上一步的变化监测结果可以得到研究区域内的变化区域，结合矿区典型地物(沉淀池、稀土矿采场、矿山建筑和引水管道等)和形变干涉图判断是否为稀土矿开采变化区域，借助已有矿权信息圈定疑似非法开采区域。

步骤5：结合形变图和相干系数图对比分析圈定区域的稀土矿非法开采情况，探究相干性系数和地物之间的关系。

2 研究区域与数据处理

2.1 研究区域概况

赣州市地处江西省南部，占据整个江西省总面积的23.6%，是中国重要的“稀土王国”[23]。定南县作为赣州市内离子型稀土主产县之一，位于114°47′49″E～115°22′48″E，24°33′37″N～25°03′21″N，地处南岭与武夷成矿带的中间部位，拥有得天独厚的成矿条件。定南县全县总面积达到1 321 km2，矿产资源非常丰富，其中包含了木子山、甲子背、长坑尾等11个稀土矿区。近几年以来，定南县稀土矿存在违法开采现象，累积了大量的矿山环境问题。因此，为了帮助当地实现可持续发展和恢复生态功能，精确、可靠地识别非法开采区域显得尤为重要。

2.2 数据处理

光学遥感影像数据本文使用了分辨率为10 m的哨兵2号数据，可在欧空局官网(https：//scihub.copernicus.eu/)免费下载。哨兵2号是高分辨率多光谱成像卫星，从可见光和近红外到短波红外，具有不同的空间分辨率，用于陆地监测，可提供植被、土壤和水覆盖、内陆水路及海岸区域等图像，还可用于紧急救援服务[24]。由于数据质量和资源的限制，选取了2017—2020年10月份的8景影像。首先需要对影像进行预处理，可根据不同地物在光学影像上表现的颜色、大小、形状等建立解译标志。

由于雷达数据与光学影像不同，属于主动遥感，影像包含的是振幅、相位和极化等多种信息。本文使用的雷达影像数据为C波段Sentinel-1A卫星的SLC数据，分辨率为30 m。为了探寻地表形变与稀土矿区开采之间的关系，选择了时间跨度为2015—2020年的影像数据，具体的雷达影像数据基本参数如表1所示。雷达数据经过预处理后可得到强度图，1)雷达影像变化监测。一般的变化监测方法只考虑处理两幅SAR图像，忽略了这些SAR图像的相位信息。作为一个相干成像系统，SAR具有更多的潜力，其相位信息应该被挖掘出来[25]。应用比值法将每一年的主影像与副影像做除法运算，可得到2015—2020年间的比值结果。图2是以疑似非法开采区域C为例的比值结果，图中的深红色区域为可能减少的信息，深绿色区域为可能增加的信息。由于地物在两景影像之间有所变化，从而造成该区域在两个时相的SAR影像的强度有所差异，同时由于地物的散射单元结构发生了变化，使得其两个时相的SAR影像相干性较弱，因此，导致了该区域的比值法结果图和相干性系数图发生变化。

表1 雷达影像数据的基本参数

其包含了地面分辨元的雷达后向散射强度信息，然后对强度图进行比值运算，即将后一时相的强度图与前一时相的强度图进行比值。最后，对雷达数据进行DInSAR处理获得相干性系数图和形变图。相干性系数图反映了两景影像之间的相干性程度，形变干涉图可反映出地表形变信息，为后续分析提供参考。

3 实验结果分析

3.1 变化监测结果与分析

开采行为会造成地物类型发生改变，因此这部分区域在SAR影像上的灰度值会产生较为明显的变化。由于雷达影像的几何畸变、斑点噪声同样会在变化监测结果中显示出来，前者是由于雷达成像时角度变化造成的，后者是由一个分辨单元内各个散射点相干回波之间的干涉效应引起的，但是这两类变化并不属于真实的稀土矿区开采变化，对于变化结果的准确性有很大影响。本文考虑到这两个因素在光学遥感影像上并不会发生变化，故采用雷达影像与光学遥感影像相结合的方法，对比分析那些雷达成像时几何畸变和斑点噪声引起的变化，最终剔除这些因素的影响。

2)光学影像变化监测。对于光学影像的动态变化监测，大部分研究都需要进行地物分类，但是中低分辨率遥感影像的分类精度并不高。由于同一地物的光谱特征相同，不同地物的光谱特征有所差异，因此，可以根据比值法对光学遥感影像进行变化监测，疑似非法开采区域C在2017—2020年间的光学遥感影像如图3所示。将变化监测结果与两景影像进行真彩色合成能够更好地分析变化区域，其中图3中的蓝色区域为可能减少的信息，黄色区域为可能增加的信息。在图3中可以看到，疑似非法开采区域C内具有沉淀池这一稀土矿山的典型地物，且在每一时期的影像上都发生了不同程度的变化，则可以证明该区域确实进行了开采活动。图4为疑似非法开采区域C在2017—2018年、2018—2019年、2019—2020年和2017—2020年的变化情况。

3)变化监测结果分析。在不受外界因素干扰的情况下，选取合适的遥感数据类型和时相，可以发现一些非法开采行为。由于雷达影像和光学影像自身的成像特点，对变化监测都会产生一定的影响，导致出现误检的情况。因此，可以充分利用两种数据的互补性进行变化监测，从而弥补光学影像受天气影响导致的局限性和改善雷达影像在几何畸变和斑点噪声区域的检测结果，开展雷达数据和光学影像在稀土矿山监测的应用研究。当仅在雷达变化监测和光学变化监测其中一个结果中出现变化而在另一个变化监测结果中表现未变化时，则判定该区域没变化；只有当雷达变化监测结果图变化，光学变化监测结果影像也变化才能说明该区域发生了变化，变化监测精度评价如表2所示。相较于雷达影像变化监测和光学影像变化监测，本文方法的准确率分别提高了11.6%和4.46%，虚警率分别降低了38.24%和23.53%。结果表明，和单一影像相比，本文提出的结合雷达和光学影像的变化监测方法结果更准确，利用遥感技术在研究区开展矿业活动开采识别动态监测是行之有效的。

表2 变化监测精度评价

要判断变化区域是否是由稀土矿开采引起的，可以根据稀土矿山开采的典型地物建立解译标志。开采的典型地物中比较容易判读的地物为沉淀池和浸矿池，也就是用来沉淀稀土的池子，主要集中在矿区旁边，形状为规则的圆形或长方形。通过对雷达影像和光学影像变化监测结果进行分析，并结合稀土矿山解译标志和矿权范围进行目视解译，圈定了9个疑似非法开采区域，疑似违法图斑圈定结果如图5所示，并以疑似非法开采区域C为例说明本文方法的有效性。图6依次为图5圈定的9个疑似非法开采区域在2017—2020年间变化监测结果真彩色合成图，变化监测结果中红色、绿色和蓝色分别代表疑似非法开采区域在2017、2018和2019年间的变化情况。

3.2 干涉相干性分析

相干图表示的是两幅影像的相关性，是最常用的相位质量图，也是最直观的质量评价图。将光学遥感影像与雷达相干性系数图进行对比可以发现，二者具有较高的一致性。在这几个干涉对中，失相干区域多为植被覆盖区域，相干性高区域多为具有稳定地物的区域。因此，不同地物在相干性系数图上表现出不同的系数值，且相干系数高低可排序为：建筑物>裸地>植被。在研究区域的相干系数图中，由于城市区域地物目标和裸地具有较高的稳定性，因此表现出相干性非常高；与此相比，林地、农田等地物覆盖区域，相干性较低，色调以蓝绿色为主；在水体等低散射回波地物区域，则表现为黑色区域。

本文通过对比疑似非法开采区域于2015—2020年间相干系数图的变化，来探究开采区域与相干性之间的关系。本文研究区域整体相干性系数范围在0.37～0.71之间，具体如图7所示。通过该图可以看出，这几处疑似非法开采区域相干性系数变化趋势具有较高的一致性，进一步证实了这些区域进行了开采活动。

3.3 地表形变干涉结果与分析

经过信号采集与数据处理，SAR影像的每一像素既包含地面分辨元的雷达后向散射强度信息，也包含与斜距有关的相位信息，将覆盖同一地区的两幅SAR影像对应像素的相位值进行差分，便可得到一个一次差分相位图，也就是干涉相位图。通过对12景覆盖研究区的SAR影像数据进行DInSAR处理，借助相干系数阈值获取相干性高的像元，并将其作为轨道精炼的控制点引入后续的处理中，最终可得到赣州市定南县在2015—2020年的地表形变信息。本文通过对比9个疑似非法开采区域在形变图中的变化，圈定了9个疑似非法开采区域于2015—2020年间的地表形变变化信息(图8)。通过该图可以看出，这几处疑似非法开采区域地表形变趋势具有较高的一致性，稀土矿的开采能够引起缓慢的地表形变，其变化范围为-25～15 mm。

4 结束语

本文在稀土矿非法开采识别过程中，结合SAR和光学遥感影像的特性，给出了一种可更有效识别疑似非法开采区域的方法。该方法综合利用了Sentinel系列影像的相干性系数、强度信息、形变信息和光谱信息等因素，对雷达和光学影像分别进行变化监测确定变化区域，并依据矿区典型地物、变化标准和已有矿权信息圈定疑似非法开采区域。相较于仅用单一影像进行矿区非法开采识别，结合SAR和光学遥感的疑似非法开采情况识别效果更优。

本文以赣州市定南县的稀土矿区为例进行实验，得到以下结论。

1)结合SAR和光学遥感的稀土矿疑似非法开采情况识别方法充分利用两种数据的互补性，不仅弥补了光学影像受天气影响的局限性，也改善了雷达影像的几何畸变和斑点噪声区域的检测结果。

2)相较于雷达影像变化监测和光学影像变化监测，本文方法的准确率分别提高了11.6%和4.46%，虚警率分别降低了38.24%和23.53%，可更准确识别出疑似稀土矿开采区域。

3)疑似稀土矿开采区域在相干系数图和地表形变图中的变化趋势具有较高一致性，相干性系数和地表形变范围分别在0.37～0.71和-25～15 mm之间。

本文揭示了雷达影像和光学影像相结合在稀土矿山开采监测中的疑似非法开采情况识别能力，同时能够给矿政管理部门治理矿山周边环境、制定矿产资源分配和确定矿产资源开发秩序等提供决策依据和技术支持。然而，本文采用的影像数据分辨率有限，未来可在本文基础上使用高分辨率影像进行下一步的研究。