一种多源时序影像的稀土矿区土地毁损与恢复分析方法与流程

本发明涉及土地恢复分析技术领域，特别是一种多源时序影像的稀土矿区土地毁损与恢复分析方法。

背景技术：

离子型稀土矿于1969年在我国江西首次发现，具有分布广、储量丰富、放射性低、稀土配分齐全，富含中重稀土元素等优势，是我国目前特有的稀土矿产资源，也是世界上稀缺的矿种，广泛分布在我国南方的江西、福建、湖南、广东、广西等省，其中江西省所占份额最大。其开采工艺主要有池浸、堆浸和原地浸矿三种。池浸、堆浸两种开采工艺由于需要将矿体和表土植被剥离，导致矿区植被遭受了巨大破坏，并且稀土开采过程中会产生大量尾砂，随着雨水流向四周，导致土地沙化，河道堵塞等问题；原地浸矿开采工艺相比池/堆浸开采工艺，对矿区植被破坏更小，但由于开采过程中需要向山体中注入大量的浸矿液体，且浸矿液体需长时间浸泡山体，容易导致土壤酸化、山体滑坡等问题。不同的稀土开采工艺导致的土地毁损状况不一样，其恢复过程也不尽相同，因此，了解不同稀土开采工艺下的土地毁损与恢复特征，为稀土矿区生态环境治理，合理有序开采提供指导，对稀土矿区的可持续发展来说显得尤为重要。

目前，矿区土地毁损与恢复分析方法主要有基于实测数据的土地毁损与恢复评价模型和基于遥感数据的土地毁损与恢复过程动态监测两种，前一种方法不仅需要大量的实测数据，费时费力，且针对不同地点需要建立不同的模型，应用推广较为困难；受技术方法和经费限制，早期的土地毁损与恢复过程动态监测多采用少数几幅间隔较长无云的遥感影像进行，研究方法为多时相分类后比较法，该方法能够反映地物类型变化在空间上的分布情况，但地物类型变化的过程了解较少；遥感数据源的不断丰富以及遥感时序分析法的发展，为改进对更广的时空范围生态过程的理解提供了重大机遇。时序分析法通过量化分析地物在时相上的变化趋势和规律来获取地物的变化情况，它对变化的判断标准具有较高的一致性，算法效率较高。

中国发明专利申请cn105598143a公开了一种资源化矿山修复方法，将尾矿库中的尾矿资源和/或剥岩区/排土场的矿石资源进行有价元素提取；无再提取价值的尾矿资源和/或矿石资源、有价元素提取工序产生的尾矿如适合深加工，则进行深加工获得目标产品；如果不适合深加工，则充填至采空区/塌陷区；对剥岩区/排土场的占用土地和/或工业场地进行无害化处理；对无害化处理后的剥岩区/排土场的占用土地和/或工业场地、充填后的采空区/塌陷区进行生态修复。

技术实现要素：

本发明需要解决的技术问题是一种高效的多源时序影像的稀土矿区土地毁损与恢复分析方法。

为解决上述的技术问题，本发明的一种多源时序影像的稀土矿区土地毁损与恢复分析方法，包括以下步骤，

步骤s1：研究区域数据获取，通过数据获取单元获取研究区域所需的原始遥感影像；

步骤s2：数据预处理，通过对获得的原始遥感影像进行辐射校正、大气校正、几何校正、图像裁剪得到研究区域的遥感影像；

步骤s3：植被指数计算，通过植被指数计算单元计算预处理后影像的ndvi；

步骤s4：转换方程构建，并对转换方程的精度进行检验；

步骤s5：nvdi转换，确定时序分析法相关参数，土地毁损与恢复定量分析。

进一步的，步骤s1中所述原始遥感影像包括空间分辨率均为30m的landsat5tm、landsat8oli和hj-1bccd遥感影像。

更进一步的，所述步骤s4具体包括以下步骤，

步骤s41：转换方程构建，基于同时期的landsat5和hj-1bccd数据、landsat8和hj-1bccd数据，构建landsat5和hj-1bccd数据、landsat8和hj-1bccd数据的ndvi之间的转换方程；

步骤s42：转换方程的精度检验，基于获得的转换方程将hj-1bccd数据的ndvi转化为与之时相相同的landsat5或landsat8的ndvi,并与真实的landsat5或landsat8比较，采用均方根误差的方法检验转换方程的精度。

更进一步的，所述步骤s5具体包括以下步骤，

步骤51：ndvi转换，基于获得的转换方程将时序影像中landsat8和hj-1bccd数据的ndvi转化为landsat5标准下，并将ndvi时序影像组合成一个多光谱影像；

步骤52：时序分析法相关参数获取，通过数据挖掘的方法获得时序分析法相关参数；

步骤53：土地毁损与恢复定量分析，根据时序分析法相关参数对土地毁损与恢复类型进行分类，获得土地毁损与恢复类型空间分布，用于定量分析稀土开采扰动下的土地毁损与恢复特征。

采用上述方法后，本发明能实时、快速、大面积的获取稀土矿区的地表信息，相比传统根据实测数据建模分析方法，效率更高，应用更为简单，推广更容易。另外，本发明采用时序分析法分析稀土开采环境下的土地毁损与恢复过程的相关特征，使得相关特征的描述更为细化。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1为本发明一种多源时序影像的稀土矿区土地毁损与恢复分析方法的流程图。

图2为本发明一种多源时序影像的稀土矿区土地毁损与恢复分析方法的结构示意图。

图3a为本发明中landsat5与hj-1bccd散点图。

图3b为本发明中landsat8与hj-1bccd散点图。

图3c为本发明中landsat5与模拟影像散点图。

图3d为本发明中landsat8与模拟影像散点图。

图4a为本发明中暂未开采和无植被覆盖像元类型的ndvi轨迹图。

图4b为本发明中采后无植被和低植被覆盖像元类型的ndvi轨迹图。

图4c为本发明中采后植被覆盖达到采前水平像元类型的ndvi轨迹图。

图5为本发明中土地毁损与恢复类型分布图。

图6a-图6f为本发明中类型1对应的高分影像图。

图7a-图7f为本发明中类型2对应的高分影像图。

图8a-图8f为本发明中类型3对应的高分影像图。

图9a-图9f为本发明中类型4对应的高分影像图。

图10a-图10f为本发明中类型5对应的高分影像图。

图11为本发明中稀土开采扰动至植被恢复时长比较图。

具体实施方式

如图1和图2所示，本发明一种多源时序影像的稀土矿区土地毁损与恢复分析方法，包括以下步骤，

步骤s1：研究区域数据获取，通过数据获取单元获取研究区域所需的原始遥感影像；，包括空间分辨率均为30m的landsat5tm、landsat8oli和hj-1bccd遥感影像。

步骤s2：数据预处理，将步骤s1中获得的遥感数据通过辐射定标模块对原始影像进行辐射定标，辐射定标后使用大气校正模块进行大气校正，大气校正后使用几何校正模块校正影像，几何校正后使用图像裁剪模块获得研究区的遥感影像。辐射定标模块中landsat5和landsat8数据使用的envi自带radiometriccalibration工具进行，hj-1bccd数据使用插件工具进行；大气校正模块使用envi自带的flaash模块进行；几何校正模块中，以landsat数据为基准，hj-1bccd数据为待校正，校正方法选择二次多项式，对hj-1bccd数据进行校正；最后通过研究区域的矢量边界对校正后的影像进行裁剪得到研究区域的遥感影像。

步骤s3：植被指数计算，通过植被指数计算单元计算预处理后影像的ndvi。获得研究区域的ndvi影像。植被指数计算单元通过envi自带的bandmath工具进行，计算公式如下：

式中ρnir和ρred分别表示遥感影像近红外、红光波段的反射率。

步骤s4：转换方程构建，并对转换方程的精度进行检验。

步骤s41：转换方程构建，基于同时期的landsat5和hj-1bccd数据、landsat8和hj-1bccd数据，采用回归分析的方法，构建landsat5和hj-1bccd数据、landsat8和hj-1bccd数据的ndvi之间的转换方程。

步骤s42：转换方程的精度检验，基于获得的转换方程将hj-1bccd数据的ndvi转化为与之时相相同的landsat5或landsat8的ndvi,并与真实的landsat5或landsat8比较，采用均方根误差的方法检验转换方程的精度。

计算公式如下:

式中xi和yi分别表示为hj-1bccd数据校正后的ndvi值、landsat5或landsat8的ndvi值，n为样本个数。

步骤5：ndvi转换，确定时序分析法相关参数，并根据相关参数对研究区域的土地毁损与恢复类型进行分类，得到不同类型在空间上的分布，从而对稀土开采扰动下的土地毁损与恢复特征进行分析。

步骤51：ndvi转换。通过获得的转换方程将时序影像中landsat8、hj-1bccd影像的ndvi转化为landsat5标准下，通过envi自带的layerstacking工具将ndvi时序影像按照时间顺序组合成一个类似多光谱影像的数据。

步骤52：确定时序分析法参数。通过目视和比较得到样本点时序ndvi的变化轨迹，如图4a-图4c。研究区域范围内包含以下5种类型，1)暂未开采，整个时期植被覆盖均保持在较高水平；2)池浸/堆浸开采，采前有植被，采后植被尚未恢复，始终维持在较低水平；3)始终无植被覆盖或植被覆盖水平较低；4)植被覆盖水平缓慢升高并恢复至周边同类区域采前水平；5)采前有植被，采后植被开始恢复，但仍低于原始水平。在时序ndvi影像上选择一定数量的五种类型的样本点，应用cart(classificationandregressiontree)决策树分类方法，将采集的样本点输入rulegenv1.02软件中，算法选择quest，得到有植被与基本无植被、已开采与未开采的分类阈值，分别为tr1＝0.5285,tr2＝0.3787。

步骤53：土地毁损与恢复定量分析。根据得到的时序分析法的相关参数tr1,tr2，按照表1的标准，对研究区域的土地毁损与恢复类型进行划分，得到土地毁损与恢复类型空间分布图及矿区毁损时间分布。

表1稀土矿区土地毁损与恢复类型划分

表格中mini表示整个观测期内时序ndvi的最低值，反映i像元所在位置是否受稀土开采扰动；maxi表示整个观测期内像元时序ndvi的最高值，反映i像元所在位置多年植被覆盖的峰值；maxi_pre表示受稀土开采扰动前时序ndvi的最高值，反映i像元所在位置采前植被覆盖的峰值；maxi_post受稀土开采扰动后时序ndvi的最高值，反映i像元所在位置采后植被覆盖的峰值；ri表示植被恢复水平，反映i像元区域采后植被覆盖是否达到或超过采前水平，ri＝maxi_post-maxi_pre,ri≥0表示采后植被覆盖水平达到或超过采前水平；ri＜0表示采后植被覆盖尚未达到采前水平。

以赣州市定南县岭北稀土矿区为例，所采用的多源时序影像如表2所示，所采用的交互比较影像对如表3所示。

表2时序影像日期、类型和标识符

表3交互比较影像对

如图4a-4c所示，类型1为暂未开采，类型2为采后无植被，类型3为无植被覆盖，类型4为采后植被覆盖达到采前水平，类型5为采后植被得到一定程度恢复，但未达到采前水平。

通过图5和图11，分析稀土开采扰动下的土地毁损与恢复过程。结合googleearth高分影像发现，类型1多为未受稀土开采影响的林地，如图6a-图6f所示；类型2为稀土开采区域且未进行矿区复垦，如图7a-图7f所示；类型3多为农田、建筑等，如图8a-图8f所示；类型4为多为退耕还林、果园种植区，如图9a-图9f所示；类型5为稀土开采区域但进行了矿区复垦，如图10a-图10f所示。

虽然以上描述了本发明的具体实施方式，但是本领域熟练技术人员应当理解，这些仅是举例说明，可以对本实施方式作出多种变更或修改，而不背离本发明的原理和实质，本发明的保护范围仅由所附权利要求书限定。