一种采矿扰动地块边界识别方法与流程
本发明属于矿山生态环境监测领域，具体涉及一种采矿扰动地块边界识别方法。
背景技术：
伴随采矿生产过程中的挖损、压占、塌陷、污染、生态修复等活动改变了矿区土地原有的地形地貌，并形成了采矿扰动地块。其中，扰动地块的边界代表了采矿活动的影响范围，扰动地块包括损毁地块和恢复地块两种类型。对采矿区域开展生态环境影响评估、执法检查、生态修复等活动均需要对采矿扰动地块的边界进行有效的识别。
现有技术中，对于扰动地块边界的识别主要基于实地采集数据或单幅遥感影像数据进行。实地采集数据或单幅遥感影像数据的实际应用成本过高，且并不能有效地适用于采矿区域使用。同时，其识别的准确度较低，并不能为生态环境的修复及管理提供可靠的技术支持。
技术实现要素：
针对上述现有技术存在的问题，本发明提供一种采矿扰动地块边界识别方法，该方法能精准有效的识别采矿扰动地块的边界，能为进行生态环境的修复及管理提供可靠的技术支持。
为了实现上述目的，本发明提供一种采矿扰动地块边界识别方法，采用长时序遥感影像进行采矿扰动边界识别，具体包括以下步骤；
步骤一：选取采矿扰动区域作为研究区，获取该区域的landsat遥感影像；
a1：对遥感影像进行辐射定标、大气校正、几何校正和影像裁剪的预处理，得到研究区的预处理影像；
a2：根据公式(1)求出各影像的ndvi指数图，记研究区域的ndvi指数图为ndvi；
式中，nir为预处理后各影像的近红外波段的反射值，r为预处理后各影像的红光波段的反射值；
步骤二：对所有ndvi指数图进行ndvi归一化处理，得到归一化处理后的ndvi指数图，记为i；
b1：设定参考影像与待归一化影像之间的回归直线方程如公式(2)所示；
式中，x为待归一化影像的ndvi值，y为参考影像的ndvi值；
b2：采用最小二乘法进行参考影像与待归一化影像之间的直线拟合，分别根据公式(3)和公式(4)求解参数
式中，
步骤三：对归一化处理后的ndvi指数图进行运算，获取研究区的ndvi差分图和ndvi绝对差分图；
c1：分别通过公式(7)和公式(8)对各个月份的可用影像按年份大小顺序进行迭代作差，获取每一个月份的ndvi差分图和ndvi绝对差分图，并分别记为md和amd；
式中，k表示该月份可用的影像年份的数量，ii为该月份可用影像按年份排序的第i个可用的归一化处理后ndvi指数；
c2：分别通过公式(9)和公式(10)将各个月份的md和amd进行累计求和，得到研究区域的ndvi差分图imaged以及ndvi绝对差分图aimaged；
式中，
步骤四：利用otsu阈值分割法获取ndvi指数差分影像的最优分割阈值，对研究区的ndvi差分图进行阈值分割；
d1：对于背景区域更暗的影像，其图像大小为x×y，统计其灰度直方图，并分别通过公式(11)和公式(12)计算感兴趣区域像素点占整幅影像的比例ω0和背景区域像素点占整幅影像的比例ω1；
ω0＝n0/(x×y)(11)；
ω1＝n1/(x×y)(12)；
式中，n0为像素值灰度小于阈值t的像素个数，n1为像素值灰度大于阈值t的像素个数，
且n0和n1满足公式(13)；ω0和ω1满足公式(14)；
n0+n1＝x×y(13)；
ω0+ω1＝1(14)；
μ0为感兴趣区域像素点占整幅影像的平均灰度；
μ1为背景区域像素点占整幅影像的平均灰度；
g为类间方差记；
μ为图像的总平均灰度；
μ＝μ0×ω0+μ1×ω1(15)；
g＝ω0(μ0-μ)2+ω1(μ1-μ)2(16)；
d2：通过公式(15)和公式(16)结合得到公式(17)，并求解类间方差g；
g＝ω0×ω1×(μ0-μ1)2(17)；
d3：在类间方差g达到最大时，阈值t即为最优阈值；以t为界限，按照公式(18)对研究区ndvi差分图及ndvi绝对差分图进行运算，分别得到研究区ndvi差分图和及ndvi绝对差分图对应的阈值分割图，记为imaget和aimaget；
d4：对阈值分割图imaget和aimaget进行图像去噪，获取采矿扰动区域处理后图像imagen和aimagen；
步骤五：采用边缘检测法对采矿扰动地块和损毁地块边界进行自动化识别；通过公式(19)和公式(20)利用roberts边缘检测算子进行采矿扰动地块和损毁地块边界进行检测；
r(x,y)＝|f(x,y)-f(x+1,y+1)|+|f(x+1,y)-f(x,y+1)|(20)；
式中，f(x,y)为输入影像图像大小(x,y)坐标处的像素值，r(x,y)为输入影像图像大小(x,y)的输出影像，分别记为边缘检测图像imager和aimager；
步骤六：将边缘检测图像imager和aimager分别转为矢量图，分别记为采矿扰动地块边界的矢量图imagev和损毁地块边界的矢量图aimagev，通过公式(21)对二者进行几何求差，获得恢复地块边界矢量图imagef；
imagef＝aimagev-imagev(21)。
作为一种优选，在步骤二中，采用伪不变特征相对辐射归一化法进行ndvi归一化，选取m个伪不变特征点，利用回归方程完成影像的ndvi归一化。
作为一种优选，采用均值滤波法对阈值分割图imaget和aimaget进行图像去噪。
本发明采用遥感影像数据进行采矿扰动地块边界识别，实现成本低，识别精度高，且覆盖范围广。由于矿区土地扰动频繁，一些地块在损毁后可能又会快速的恢复。传统的单时相遥感影像分类方法、双时相遥感影像变化检测方法，并不能全面反映所有扰动过的区域，即，可能会漏检一些损毁后快速恢复的区域，并且，需要夏季的遥感影像来反映地表变化，当一些区域缺少夏季影像时，就无法完成扰动地块的识别。利用多时相ndvi差分和ndvi绝对差分的方法提取扰动斑块，一方面可以准确提取历史累计扰动边界，另一方面，该方法对遥感影像的获取时间要求低，可以利用所有季节的可用影像，可以适用于可用遥感影像较少的区域，从而可以有效识别出损毁后恢复的区域，能更全面地反映所有扰动过的区域。因此，本发明不仅可以识别出扰动地块，还可以进一步将采矿扰动地块划分为损毁地块与恢复地块，这可为矿区生态环境评价工作提供数据基础，对实现矿区生态环境恢复及管理具有重要意义。本发明提出的方法能精准地识别出采矿扰动地块的边界，可为矿区采矿规划提供可靠的基础依据，解决了现有技术对因采矿引起的扰动地块边界的识别不准确的问题。
附图说明
图1是本发明的流程示意图；
图2为本发明中实施例的ndvi差分图；
图3为本发明中实施例的ndvi绝对差分图；
图4为本发明中实施例的ndvi差分阈值分割图；
图5为本发明中实施例的ndvi绝对差分阈值分割图；
图6为本发明中实施例的ndvi差分滤波去噪图；
图7为本发明中实施例的ndvi绝对差分滤波去噪图；
图8为本发明中实施例的ndvi差分边缘检测图；
图9为本发明中实施例的ndvi绝对差分边缘检测图；
图10为本发明中实施例的采矿扰动地块边界矢量图；
图11为本发明中实施例的损毁地块边界矢量图；
图12为本发明中实施例的恢复地块边界矢量图。
具体实施方式
下面结合附图对本发明作进一步说明。
如图1所示，一种采矿扰动地块边界识别方法，采用长时序遥感影像进行采矿扰动边界识别，具体包括以下步骤；
步骤一：选取采矿扰动区域作为研究区，获取该区域的landsat遥感影像；
a1：对遥感影像进行辐射定标、大气校正、几何校正和影像裁剪的预处理，得到研究区的预处理影像；
a2：根据公式(1)求出各影像的ndvi指数图，记研究区域的ndvi指数图为ndvi；
式中，nir为预处理后各影像的近红外波段的反射值，r为预处理后各影像的红光波段的反射值；
步骤二：对所有ndvi指数图进行ndvi归一化处理，得到归一化处理后的ndvi指数图，记为i；作为一种优选，采用伪不变特征相对辐射归一化法进行ndvi归一化，选取m个伪不变特征点，利用回归方程完成影像的ndvi归一化；作为一个实施例，可以以某研究区landsat8遥感影像作为参考影像，并选取15个伪不变特征点，此时，m＝15；该实施例的可用影像的月份的个数为8；
b1：设定参考影像与待归一化影像之间的回归直线方程如公式(2)所示；
式中，x为待归一化影像的ndvi值，y为参考影像的ndvi值；
b2：采用最小二乘法进行参考影像与待归一化影像之间的直线拟合，分别根据公式(3)和公式(4)求解参数
式中，
在m＝15时，
步骤三：对归一化处理后的ndvi指数图进行运算，获取研究区的ndvi差分图和ndvi绝对差分图；
c1：分别通过公式(7)和公式(8)对各个月份的可用影像按年份大小顺序进行迭代作差，获取每一个月份的ndvi差分图和ndvi绝对差分图，并分别记为md和amd；
式中，k表示该月份可用的影像年份的数量，ii为该月份可用影像按年份排序的第i个可用的归一化处理后ndvi指数；
c2：分别通过公式(9)和公式(10)将各个月份的md和amd进行累计求和，得到研究区域的ndvi差分图imaged以及ndvi绝对差分图aimaged，对于本申请中的实施例，其对应图像分别如图2和图3所示；
式中，
在n＝8时，
步骤四：利用otsu阈值分割法获取ndvi指数差分影像的最优分割阈值，对研究区的ndvi差分图进行阈值分割；
d1：对于背景区域更暗的影像，其图像大小为x×y，统计其灰度直方图，并分别通过公式(11)和公式(12)计算感兴趣区域像素点占整幅影像的比例ω0和背景区域像素点占整幅影像的比例ω1；
ω0＝n0/(x×y)(11)；
ω1＝n1/(x×y)(12)；
式中，n0为像素值灰度小于阈值t的像素个数，n1为像素值灰度大于阈值t的像素个数，且n0和n1满足公式(13)；ω0和ω1满足公式(14)；
n0+n1＝x×y(13)；
ω0+ω1＝1(14)；
μ0为感兴趣区域像素点占整幅影像的平均灰度；
μ1为背景区域像素点占整幅影像的平均灰度；
g为类间方差记；
μ为图像的总平均灰度；
μ＝μ0×ω0+μ1×ω1(15)；
g＝ω0(μ0-μ)2+ω1(μ1-μ)2(16)；
d2：通过公式(15)和公式(16)结合得到公式(17)，并求解类间方差g；
g＝ω0×ω1×(μ0-μ1)2(17)；
d3：在类间方差g达到最大时，阈值t即为最优阈值；以t为界限，按照公式(18)对研究区ndvi差分图及ndvi绝对差分图进行运算，分别得到研究区ndvi差分图和及ndvi绝对差分图对应的阈值分割图，记为imaget和aimaget，对于本申请中的实施例，其对应的结构分别如图4和图5所示；
d4：对阈值分割图imaget和aimaget进行图像去噪，获取采矿扰动区域处理后图像imagen和aimagen；作为一种优选，采用均值滤波法对阈值分割图imaget和aimaget进行图像去噪；对于本申请中的实施例，其对应的结果分别如图6和图7所示；
步骤五：采用边缘检测法对采矿扰动地块和损毁地块边界进行自动化识别；通过公式(19)和公式(20)利用roberts边缘检测算子进行采矿扰动地块和损毁地块边界进行检测；
r(x,y)＝|f(x,y)-f(x+1,y+1)|+|f(x+1,y)-f(x,y+1)|(20)；
式中，f(x,y)为输入影像图像大小(x,y)坐标处的像素值，r(x,y)为输入影像图像大小(x,y)的输出影像，分别记为边缘检测图像imager和aimager；
对于本申请中的实施例，其结果分别如图8和图9所示。
步骤六：将边缘检测图像imager和aimager分别转为矢量图，分别记为采矿扰动地块边界的矢量图imagev和损毁地块边界的矢量图aimagev，对于本申请中的实施例，其结果分别如图10和图11所示。通过公式(21)对二者进行几何求差，获得恢复地块边界矢量图imagef；对于本申请中的实施例，其结果分别如图12所示。
imagef＝aimagev-imagev(21)。
- 还没有人留言评论。精彩留言会获得点赞！