CN108805920A - 未利用地上土壤污染风险区域的识别方法和装置 - Google Patents

Abstract

本发明提供一种未利用地上土壤污染风险区域的识别方法和装置，其中方法包括：根据监测区域在两个时相的遥感图像中每个像素点的灰度，获得反映各遥感图像的纹理特征的纹理特征影像；根据两时相遥感图像的灰度差异以及纹理特征影像的纹理差异，获得监测区域内灰度和纹理特征变化均大于预设阈值的区域，作为目标区域；根据监测区域内的未利用地的范围，利用高分辨率影像对所述目标区域进行识别，获得未利用地上的土壤污染风险区域。本发明实施例克服对监测区域进行逐像元查找、缺乏针对性的弊端，极大地减少了工作量，提高了监测效率。

Description

未利用地上土壤污染风险区域的识别方法和装置

技术领域

本发明涉及环境监测技术领域，更具体地，涉及未利用地上土壤污染风险区域的识别方法。

背景技术

未利用地是指直接用于农业生产的耕地、园地、林地、牧草地、其他农用地以及商业、交通运输业、居民点和工矿、公共事业等建设用地以外的土地。滩涂、荒草地、裸地等未利用地大多属于生态脆弱地区，一旦遭受污染，恢复和治理难度极大，加强对未利用地的监管对于保持土地利用结构、提高土地的利用率、改善生态环境具有重要的意义。近年来，未利用地污染问题日益突出，我国大部分地区都存在未利用地土壤污染的现象，由于人类活动的影响使得局部范围内土壤污染物含量严重超标，防止土壤污染、改善未利用地土壤环境是我国现阶段发展亟需采取的措施。因此加强未利用地上污染企业、污染物排放渣场的监测具有重要意义。

传统的土壤污染调查通过去实地采集样本，对土壤和土壤周围水体、沉积物进行化学分析，然而未利用地污染往往涉及范围广阔，需要耗费大量的人力、物力和时间，很难进行实时、大区域土壤污染监测。之前的方法多通过实地采样调查、研究植被生长状况来分析是否为土壤污染风险点，或者采用目视解译方式判断土壤污染风险区域。

发明内容

本发明提供一种克服上述问题或者至少部分地解决上述问题的未利用地上土壤污染风险区域的识别方法和装置

根据本发明的一个方面，提供一种未利用地上土壤污染风险区域的识别方法，包括：

根据监测区域在两个时相的遥感图像中每个像素点的灰度，获得反映各遥感图像的纹理特征的纹理特征影像；

根据两时相遥感图像的灰度差异以及纹理特征影像的纹理差异，获得监测区域内灰度和纹理特征变化均大于预设阈值的区域，作为目标区域；

根据监测区域内的未利用地的范围，利用高分辨率影像对所述目标区域进行识别，获得未利用地上的土壤污染风险区域。

优选地，所述根据两时相遥感图像的灰度差异以及两时相纹理特征影像的纹理差异，获得监测区域内灰度和纹理特征变化均大于预设阈值的区域，作为目标区域的步骤，具体包括：

根据两时相的遥感图像中对应位置像素点的灰度，获得监测区域内灰度差值绝对值大于预设阈值的区域，作为第一区域；

根据两时相的纹理特征影像中对应位置像素点的纹理特征值，获得监测区域内纹理特征值差值的绝对值大于预设阈值的区域，作为第二区域；

将第一区域与第二区域的重合区域作为所述目标区域。

优选地，所述获得反映各遥感图像的纹理特征的纹理特征影像的步骤之前，还包括获得所述监测区域在两个时相的遥感图像的步骤：

获取覆盖监测区域的云量小于预设门限的两时相的Landsat影像；

对所述两时相的Landsat影像分别进行辐射定标、大气校正以及几何校正处理，并叠加矢量图像，裁剪出监测区域在两个时相的原始遥感影像；

对两个时相的原始遥感影像分别进行主成分分析，提取前一定个数的主成分，构成对应时段的所述遥感影像。

优选地，所述根据监测区域内的未利用地的范围，利用高分辨率影像对所述目标区域进行识别之前还包括：获取监测区域内的未利用地的步骤：

根据遥感影像的光谱特征，结合野外实测资料和高分辨率影像确定真实地物类别；

根据真实地物在遥感影像上的色调、形状、纹理及空间位置特征，建立未利用地解译标志；

根据所述未利用地解译标志进行目视解译，确定监测区域内的未利用地的矢量边界。

优选地，所述获得未利用图上的土壤污染风险区域的步骤之后，还包括：

利用POI数据对所述土壤污染风险区域进行空间验证，判别污染风险点的污染源。

优选地，所述根据监测区域在两个时相的遥感图像中每个像素点的灰度，获得反映各遥感图像的纹理特征的纹理特征影像的步骤，具体包括：

确定所述两时相的遥感图像的灰度等级，选定若干个特征量作为纹理统计量；

利用滑动窗口对所述遥感图像截取子图像，对每次截取的子图像利用灰度共生矩阵法计算所述纹理统计量，并将所述纹理统计量作为对应子图像中心像素的纹理特征值；

获取所有像素点的纹理特征值，对任意一个像素点，将该像素点的纹理特征值中的各特征量进行加权求和，作为该像素点的新灰度值，以获得纹理特征影像。

优选地，所述纹理统计量包括以下特征量：能量、对比度、逆差矩以及熵。

根据本发明的另一个方面，还提供一种未利用地上土壤污染风险区域的识别装置，包括：

纹理特征影像获取模块，用于根据监测区域在两个时相的遥感图像中每个像素点的灰度，获得反映各遥感图像的纹理特征的纹理特征影像；

目标区域获取模块，用于根据两时相遥感图像的灰度差异以及纹理特征影像的纹理差异，获得监测区域内灰度和纹理特征变化均大于预设阈值的区域，作为目标区域；

识别模块，用于根据监测区域内的未利用地的范围，利用高分辨率影像对所述目标区域进行识别，获得未利用地上的土壤污染风险区域。

优选地，所述目标区域获取模块具体用于：

根据两时相的遥感图像中对应位置像素点的灰度，获得监测区域内灰度差值绝对值大于预设阈值的区域，作为第一区域；

根据两时相的纹理特征影像中对应位置像素点的纹理特征值，获得监测区域内纹理特征值差值的绝对值大于预设阈值的区域，作为第二区域；

将第一区域与第二区域的重合区域作为所述目标区域。

优选地，所述识别装置还包括：

污染源识别模块，用于利用POI数据对所述土壤污染风险区域进行空间验证，判别污染风险点的污染源。

本发明提出的未利用地上土壤污染风险区域的识别方法和装置，根据监测区域在两个时相的遥感图像中每个像素点的灰度，获得反映各遥感图像的纹理特征的纹理特征影像，遥感图像中每个像素点的灰度值能够反应区域对应位置的地表情况，通过获取两个时相的遥感图像的纹理特征影像，也就相当于获得了监测区域在两个时相上的地物分布情况。根据两时相遥感图像的灰度差异以及两时相纹理特征影像的纹理差异，获得监测区域内灰度和纹理特征变化均大于预设阈值的区域，作为目标区域，能够快速地识别出这些变化大的地方，同时基于纹理特征影像的性质，即使两个时相的遥感图像分辨率不同，也不会影响后续获取目标区域的结果，进一步增强了本发明实施例的识别方法的适应性。根据监测区域内的未利用地的范围，利用高分辨率影像对目标区域进行识别，获得未利用地上的土壤污染风险区域，克服对监测区域进行逐像元查找、缺乏针对性的弊端，极大地减少了工作量，提高了监测效率。

附图说明

图1为根据本发明实施例所提供的未利用地上土壤污染风险区域的识别***的场景示意图；

图2为根据本发明实施例的未利用地上土壤污染风险区域的识别方法的流程示意图；

图3为根据本发明实施例的通过遥感影像获得纹理特征影像的流程示意图；

图4为根据本发明实施例的对遥感图像进行预处理的流程示意图；

图5为根据本发明实施例的通过目视解释获得监测区域内未利用地的流程示意图；

图6为根据本发明实施例的未利用地上土壤污染风险区域的识别装置的功能框图。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

本发明实施例提供一种未利用地上土壤污染风险区域的识别方法和装置。为了便于理解，在此先对本实施例及后续实施例中可能涉及到的相关概念进行解释说明：

主成分分析法：主成分分析也称主分量分析，旨在利用降维的思想，把多指标转化为少数几个综合指标。在数学变换中保持变量的总方差不变，使第一变量具有最大的方差，称为第一主成分，第二变量的方差次大，并且和第一变量不相关，称为第二主成分。依次类推，I个变量就有I个主成分。

灰度共生矩阵(GLCM)：由于纹理是由灰度分布在空间位置上反复出现而形成的，因而在图像空间中相隔某距离的两像素之间会存在一定的灰度关系，即图像中灰度的空间相关特性。灰度共生矩阵就是一种通过研究灰度的空间相关特性来描述纹理的方法。灰度共生矩阵是涉及像素距离和角度的矩阵函数，它通过计算图像中一定距离和一定方向的两点灰度之间的相关性，来反映图像在方向、间隔、变化幅度及快慢上的综合信息。灰度直方图是对图像上单个像素具有某个灰度进行统计的结果，而灰度共生矩阵是对图像上保持某距离的两像素分别具有某灰度的状况进行统计得到的。假设原图灰度级为L，则灰度共生矩阵是一个LxL大小的计数矩阵，每个点的值为满足位置条件的“灰度对”个数，设这个矩阵中坐标为(p,q)的点处值为num，则表示原图中满足条件的灰度对(p,q)有num个，即位置关系满足“条件”的灰度为p和q的像素对个数为num。

矢量图像，也称为面向对象的图像或绘图图像。矢量文件中的图形元素称为对象。每个对象都是一个自成一体的实体，它具有颜色、形状、轮廓、大小和屏幕位置等属性。既然每个对象都是一个自成一体的实体，就可以在维持它原有清晰度和弯曲度的同时，多次移动和改变它的属性，而不会影响图例中的其它对象。这些特征使基于矢量的程序特别适用于图例和三维建模，因为它们通常要求能创建和操作单个对象。基于矢量的绘图同分辨率无关。这意味着它们可以按最高分辨率显示到输出设备上，因此在印刷时可以任意放大或缩小图形而不会影响出图的清晰度。

请参阅图1，图1为本发明实施例所提供的未利用地上土壤污染风险区域的识别***的场景示意图，该识别***可以包括识别装置，该识别装置具体可以集成在平板电脑、手机、笔记本电脑、台式电脑等具备储存单元并安装有微处理器而具有运算能力的终端中，主要用于根据监测区域在两个时相的遥感图像中每个像素点的灰度，获得反映各遥感图像的纹理特征的纹理特征影像；根据两时相遥感图像的灰度差异以及纹理特征影像的纹理差异，获得监测区域内灰度和纹理特征变化均大于预设阈值的区域，作为目标区域；根据监测区域内的未利用地的范围，利用高分辨率影像对目标区域进行识别，获得未利用地上的土壤污染风险区域。

此外，该识别***还可以包括服务器，识别装置可以将识别出未利用地上的土壤污染风险区域的影像文件上传至服务器，由服务器存储该影像文件，当识别装置需要获取该影像文件时，可以从服务器获取等。

以下分别进行详细说明。

在本实施例中，将从识别方法的角度进行描述。一种未利用地上土壤污染风险区域的识别方法，包括：根据监测区域在两个时相的遥感图像中每个像素点的灰度，获得反映各遥感图像的纹理特征的纹理特征影像；根据两时相遥感图像的灰度差异以及纹理特征影像的纹理差异，获得监测区域内灰度和纹理特征变化均大于预设阈值的区域，作为目标区域；根据监测区域内的未利用地的范围，利用高分辨率影像对所述目标区域进行识别，获得未利用地上的土壤污染风险区域。

请参阅图2，图2是本发明实施例的未利用地上土壤污染风险区域的识别方法的流程示意图，如图所示，该识别方法包括：

S101、根据监测区域在两个时相的遥感图像中每个像素点的灰度，获得反映各遥感图像的纹理特征的纹理特征影像。

本发明实施例中的识别方法可以应用在识别未利用地上的具有较高可能被污染的区域(即土壤污染风险区域)的场景。因此，监测区域显然包括了未利用地，之所以需要获取监测区域在两个时相的遥感图像，目的是通过对比两个时间段的监测区域的遥感图像，获得监测区域里发生了变化的区域，而发生了变化的区域就可能是因为受到污染而从遥感影像中可以观察出区别的区域。由于Landsat系列卫星具有较长的时间序列，为利用近四十年代遥感影像监测地球资源及环境提供了很好的数据支持。作为本领域技术人员可以理解的是，由灰度共生矩阵获得的纹理特征影像可以反应出灰度的空间相关特性，在本发明实施例中，遥感图像中每个像素点的灰度值能够反应区域对应位置的地表情况，通过获取两个时相的遥感图像的纹理特征影像，也就相当于获得了监测区域在两个时相上的地物分布情况。

S102、根据两时相遥感图像的灰度差异以及两时相纹理特征影像的纹理差异，获得监测区域内灰度和纹理特征变化均大于预设阈值的区域，作为目标区域。

由上述内容可知，遥感图像的灰度可以在一定程度上反应真实地表的情况，单纯利用两时相的遥感图像的灰度差异可以直接识别出监测区域在两时相上的改变，只用灰度进行识别就是简单对比灰度变化，地表的细节也不能得到有效的凸显，因此这种方式是不准确的。而纹理特征影像则是对灰度信息的再加工，能够体现出图像的纹理。

由于遥感图像的纹理是对像元之间关系的一种客观描述，若将一幅遥感图像划分为多幅子图像，在每幅子图像都会有相似的像元关系出现，也就是在每幅图像中出现相同的纹理，从整幅遥感图像上能够发现有相同纹理出现多次，由此表明纹理具有自相似性；另一方面，同一地区可以有不同分辨率的遥感图像，分辨率会影响纹理的粗糙程度，对于每一种纹理，粗糙程度只是反映局部具有相同关系的像元对个数的变化，而像元间的关系没有改变，比如像元间的灰度差，相对位置等，所以纹理具有标度不变性。结合到本发明实施例的应用场景中，可以理解的是，出于获取未利用地的土壤污染风险区域的考虑，监测区域本身的地理位置会比较偏僻，并且地理环境偏向于自然、原始、少被人类意志所改造，因此这种区域的环境中地物分布一般是均匀的，纹理特征影像的纹理会比较自然，若根据两个时相的灰度共生矩阵获知纹理特征变化大，说明对应的区域很可能发生了不寻常的变化，采用上述方法能够快速地识别出这些变化大的地方，同时基于纹理特征影像的性质，即使两个时相的遥感图像分辨率不同，也不会影响后续获取目标区域的结果，进一步增强了本发明实施例的识别方法的适应性。

S103、根据监测区域内的未利用地的范围，利用高分辨率影像对目标区域进行识别，获得未利用地上的土壤污染风险区域。

需要说明的是，本发明实施例在获得了监测区域内变化较大的区域，即目标区域后，结合监测区域内的未利用地获得目标区域中的未利用地，即未利用地中变化较大的区域，再利用高分辨率影像对未利用地中变化较大的区域进行目视解译和空间关系分析，最终获得未利用地上的土壤污染风险区域，本发明实施例克服对监测区域进行逐像元查找、缺乏针对性的弊端，极大地减少了工作量，提高了监测效率。并且，由于未利用地具有环境单一的特征，因此将未利用地反映至图像中会出现整体的颜色变化较小的特点，本发明实施例通过利用高分辨率影像对未利用地进行目视解释和空间关系分析时，可以快速地从图像中颜色变化大的区域中获得未利用地上的土壤污染风险区域，能够大大提高监测效率。

在上述实施例的基础上，步骤S102具体包括：

S1021、根据两时相的遥感图像中对应位置像素点的灰度，获得监测区域内灰度差值绝对值大于预设阈值的区域，作为第一区域；

具体地，两个时相的遥感图像的大小均为M×N，将前一个时相的遥感图像中第m行第n列的像素点的灰度与后一个时相的遥感图像中的第m行第n列的像素点的灰度求差，将差值的绝对值作为对应位置的像素点的灰度差异，将该灰度差异与预设阈值进行比对，如果该灰度差异大于了预设阈值，就将该像素点对应的监测区域内的真实区域作为第一区域的一部分，这样通过遍历所有像素点，就获得了第一区域。显然，M、N、m和n均为正整数，并且m∈M，n∈N。

S1022、根据两时相的纹理特征影像中对应位置像素点的纹理特征值，获得监测区域内纹理特征值差值的绝对值大于预设阈值的区域，作为第二区域；

具体地，两个时相的纹理特征影像的大小均为M×N，将前一个时相的纹理特征影像中第m行第n列的像素点的纹理特征值与后一个时相的纹理特征影像中的第m行第n列的像素点的纹理特征值求差，将差值的绝对值作为对应位置的像素点的纹理差异，将该纹理差异与预设阈值进行比对，如果该纹理差异大于了预设阈值，就将该像素点对应的监测区域内的真实区域作为第二区域的一部分，这样通过遍历所有像素点，就获得了第二区域。

S1023、将第一区域与第二区域的重合区域作为所述目标区域。

需要说明的是，由于未利用地具有环境单一的特征，因此将未利用地中未变化区域反映至图像中会出现整体的纹理及灰度变化较小的特点，通过对比两个时相的监测区域遥感图像，获得监测区域内发生变化的区域，即疑似受污染区域。

本发明实施例通过对比遥感图像的灰度差异以及纹理特征影像的纹理差异，实现了获取目标区域的作用。并且，由于目标区域是指灰度异常区域与纹理异常区域的重合区域，因此在实际应用时，获取目标区域的方法还可以采用以下形式：

根据两时相的遥感图像中对应位置像素点的灰度，获得监测区域内灰度差值绝对值大于预设阈值的区域，作为第三区域；

根据两时相的纹理特征影像中对应第三区域的像素点的纹理特征值，获得第三区域中纹理特征值差值的绝对值大于预设阈值的区域，作为第四区域；

将所述第四区域作为目标区域。

本发明实施例区别于上一个实施例的地方在于，首先通过灰度差异对监测区域进行筛选，获得一个范围较小的区域，即第三区域，再根据两时相的纹理特征影像中对应第三区域的像素点的纹理特征值，进一步对第三区域进行筛选，获得第四区域，并将第四区域作为目标区域，相比上一个实施例缩小了识别面积，提供了识别效率。

在上述实施例的基础上，本发明实施例主要描述获得纹理特征影像的具体步骤，参见图3，图3为本发明实施例通过遥感影像获得纹理特征影像的流程示意图，如图所示，该流程具体包括：

S201、确定两时相遥感图像的灰度等级，选定若干个特征量作为纹理统计量。

需要说明的是，为了减少计算量和计算时间，首先需要将两时相的遥感图像的灰度级统一压缩至一定的级数。由于灰度共生矩阵会随着选择的两个点的取值变化而变化，因此，在具体实施时，本发明实施例会预先确定两个点的角度和距离，以使得两个时相的灰度共生矩阵是相互对应的。

一般来说，如果图像的是由具有相似灰度值的像素块构成，则灰度共生矩阵的对角元素会有比较大的值；如果图像像素灰度值在局部有变化，那么偏离对角线的元素会有比较大的值。通常可以用一些标量来表征灰度共生矩阵的特征，即特征量，常用的特征量有：

ASM(能量，angular second moment)，也即每个矩阵元素的平方和，反映了图像灰度分布均匀程度和纹理粗细度。如果灰度共生矩阵的所有值均相等，则ASM值小；相反，如果其中一些值大而其它值小，则ASM值大。当共生矩阵中元素集中分布时，此时ASM值大。ASM值大表明一种较均一和规则变化的纹理模式。如果灰度共生矩阵中的值集中在某一块(比如对连续灰度值图像，值集中在对角线；对结构化的图像，值集中在偏离对角线的位置)，则ASM有较大值，若灰度共生矩阵中的值分布较均匀(如噪声严重的图像)，则ASM有较小的值。

对比度(contrast，CON)，直接反映了某个像素值及其领域像素值的亮度的对比情况。如果偏离对角线的元素有较大值，即图像亮度值变化很快，则CON会有较大取值，这也符合对比度的定义。对比度也反映了图像的清晰度和纹理沟纹深浅的程度。纹理沟纹越深，其对比度越大，视觉效果越清晰；反之，对比度小，则沟纹浅，效果模糊。灰度差即对比度大的象素对越多，这个值越大。灰度共生矩阵中远离对角线的元素值越大，CON越大。

逆差矩(inverse different moment，IDM)，反映图像纹理的同质性，度量图像纹理局部变化的多少。其值大则说明图像纹理的不同区域间缺少变化，局部非常均匀。如果灰度共生矩阵对角元素有较大值，IDM就会取较大的值。因此连续灰度的图像会有较大IDM值。

熵(entropy，ENT)，熵是图像所具有的信息量的度量，纹理信息也属于图像的信息，是一个随机性的度量，它表示了图像中纹理的非均匀程度或复杂程度。若灰度共生矩阵值分布均匀，也即图像近于随机或噪声很大，熵会有较大值。

自相关(correlation，COR)，反应了图像纹理的一致性。如果图像中有水平方向纹理，则水平方向矩阵的COR大于其余矩阵的COR值。它度量空间灰度共生矩阵元素在行或列方向上的相似程度，因此，相关值大小反映了图像中局部灰度相关性。当矩阵元素值均匀相等时，相关值就大；相反，如果矩阵像元值相差很大则相关值小。

根据灰度共生矩阵的若干个特征量，能够更准确地表达地表地物的分布，进一步地，根据两时相对应的灰度共生矩阵的特征量，也就可以准确地获得地标地物的变化情况，同时结合灰度变化程度，获得监测区域内灰度和纹理特征变化均大于预设阈值的区域，作为目标区域。通过上述步骤，本发明实施例可以在监测区域内初步找出可能出现监测区域内出现污染的区域，这种将灰度差异与灰度共生矩阵相结合的方式能够更准确地找出监测区域内环境发生变化的区域。

S202、利用滑动窗口对所述遥感图像截取子图像，对每次截取的子图像利用灰度共生矩阵法计算所述纹理统计量，并将所述纹理统计量作为对应子图像中心像素的纹理特征值。

具体地，设遥感图像的大小为M×N，滑动窗口的大小为L×L，其中L＜M，L＜N，利用滑动窗口由遥感图像的一个端点开始滑动，每次滑动将截取一个大小为L×L的子图像，然后在该子图像中计算0°、45°、90°和135°方向上的特征量，例如ASM、CON、IDM、ENT，对于每一种特征量，将该特征量4个方向的平均值作为结果，显然，每一个子图像的中心像素的纹理特征值的个数与特征量的种类一致。

S203、获取所有像素点的纹理特征值，对任意一个像素点，将该像素点的纹理特征值中的各特征量进行加权求和，作为该像素点的新灰度值，以获得纹理特征影像。

例如，在获得每个像素点的纹理特征值后，若纹理特征值包含ENT、COR、IDM以及CON，分别对ENT、COR、IDM以及CON赋予权重值，然后将每个特征量乘以权重后做总和，将求和后的值作为像素点的新灰度值，通过统计每个像素点的新灰度值，就可以获得纹理特征影像。具体权重值的设置可以通过预先测试一定数量的图像，判断各特征量对识别的影响大小，根据各特征量的影响设置权重值。

由于两个时相的遥感影像或多或少会存在拍摄环境不同的情况，如果直接以从***等单位直接获取的遥感影像进行识别，两个时相的遥感影像的灰度并不是基于同一个标准建立的，进而会在后续步骤中产生误差，因此，本发明实施例在步骤S101之前，还包括对遥感图像进行预处理的步骤。参见图4，图4为对遥感图像进行预处理的流程示意图，如图所示，该步骤包括：

S301、获取覆盖监测区域的云量小于预设门限的两时相的Landsat影像；

需要说明的是，由于目前卫星的透云透雾能力还不完善，在选择遥感影像时需要尽量避免有云天气，可以设置门限值为10％，即只选取覆盖监测区域的云量小于10％的遥感图像，具体地界定监测区域的云量的方法已被本领域技术人员所熟知，在此不再赘述。Landsat影像是指美国NASA陆地探测卫星***获取的遥感图像。

S302、对所述两时相的Landsat影像分别进行辐射定标、大气校正以及几何校正处理，并叠加矢量图像，裁剪出监测区域在两个时相的原始遥感影像。

需要说明的是，辐射定标是将传感器所得的测量值转化为绝对亮度或变换为与地表反射率、表面温度等物理量有关的相对值的处理过程，能够消除传感器本身的误差，确定传感器入口处的准确辐射值。辐射定标主要有三种定标方法，分别为实验室定标、机上/星上定标以及场地定标。大气校正将辐射亮度或者表面反射率转换为地表实际反射率，能够消除大气散射、吸收、反射引起的误差。现有的成熟技术主要有三种，分别为基于辐射传输模型的校正方法、基于简化辐射传输模型的黑暗像元法以及基于统计学模型的反射率反演。几何校正利用GCP纠正各种因素引起的几何变形，对影像进行地理坐标定位，获得真实坐标信息。主要包括：空间像元位置的变换和变换后像元亮度值的计算。几何校正用于纠正***和非***因素引起的几何畸变。对两时相的Landsat影像进行图像校正的流程顺序为：辐射定标——大气校正——几何校正，经过校正后的影像降低了提取纹理、灰度特征变化区域时的误差，使相对稳定的同名地物的辐射值在不同时相的遥感图像上保持一致。

S303、对两个时相的原始遥感影像分别进行主成分分析，提取前一定个数的主成分，构成对应时段的所述遥感影像。

需要说明的是，为减少后续处理计算量，对步骤S302处理后的图像进行主成分分析，主成分分析用于将原来影像中的大部分信息集中到尽可能少的几个新主成分影像中，并使各个主成分包含的信息彼此不重叠，从而减少总数据量，消除冗余信息，突出有效信息。经验证，第一主成分对原始数据的贡献率均在85％以上，因此在一个可选实施例中选择第一主成分作为对应时段的遥感影像。

在上述实施例的基础上，监测区域内未利用地可以通过目视解释获得，具体地，参见图5，图5为通过目视解释获得监测区域内未利用地的流程示意图，如图所示，该步骤包括：

S401、根据遥感影像的光谱特征，结合野外实测资料和高分辨率影像确定真实地物类别；

需要说明的是，不同波段对地物各有不同的响应，例如MSS₄波段对水体有较强的透射能力，可用于识别水下地物，因MSS₅波段对叶绿素的吸收作用，植被在图像中具有较暗的色调，叶绿素在MSS₆波段具有较强的反射能量，适合反映地形，对土地资源解译极为有利，因此，可以根据遥感影像的光谱特征对遥感影像中的地物进行初步分类，再根据野外实测的资料以及高分辨率影像，确定遥感影像中地物的准确类别。

S402、根据真实地物在遥感影像上的色调、形状、纹理及空间位置特征，建立未利用地解译标志。一般来说，解译标志有形状(Shape)、大小(Size)、色调(Tone)和色彩(Color)、阴影(Shadow)、水系(River System)、地貌形态标志、纹理(Texture)和位置(Location)。

S403、根据未利用地解译标志进行目视解译，确定监测区域内的未利用地的矢量边界。

需要说明的是，本发明实施例实现了准确获得监测区域内未利用地的目的，为后续识别未利用地中的土壤污染风险区域确定了识别范围，并且在建立未利用地解译标志的过程中，也对未利用地上的污染情况进行了初步的了解，使得后续在识别土壤污染风险区域时的效率更高。

通过上述实施例获得未利用地上的土壤污染风险区域后，为了进一步判断疑似污染区域，分析出相关的污染源，本发明实施例的识别方法还包括：利用POI数据对土壤污染风险区域进行空间验证，判别污染风险点的污染源。

具体地，利用POI数据对土壤污染风险区域进行空间验证分析，结合Arcgis软件中的近邻分析工具，通过设置搜寻半径来查找最邻近的POI点位，确定污染区域是否为附近企业造成，并分析造成污染的相应污染源。

POI数据是真实地物的代表，POI数据包含了地物的名称、经纬度、类别信息，它甚至可以代表在空间中占一定面积的地理存在，它的主要内容包括：名称、类别、经纬度、地址名称，获取POI数据的手段多样，例如高德地图、报读地图的POI数据均可以开源获得，POI数据分类较细，并且可以根据实际情况定时更新。

在本实施例中，将从识别装置的角度进行描述，该识别装置具体可以集成在电脑、平板电脑、手机等具有存储单元并安装有处理器或微处理器而具有运算能力的终端中。

请参阅图6，图6是本发明实施例的未利用地上土壤污染风险区域的识别装置的功能框图，如图所示，该识别装置包括：

纹理特征影像获取模块501，用于根据监测区域在两个时相的遥感图像中每个像素点的灰度，获得反映各遥感图像的纹理特征的纹理特征影像。

本发明实施例中的识别方法可以应用在识别未利用地上的具有较高可能被污染的区域(即土壤污染风险区域)的场景。因此，监测区域显然包括了未利用地，之所以需要获取监测区域在两个时相的遥感图像，目的是通过对比两个时间段的监测区域的遥感图像，获得监测区域里发生了变化的区域，而发生了变化的区域就可能是因为受到污染而从遥感影像中可以观察出区别的区域。作为本领域技术人员可以理解的是，由灰度共生矩阵获得的纹理特征影像可以反应出灰度的空间相关特性，在本发明实施例中，遥感图像中每个像素点的灰度值能够反应区域对应位置的地表情况，通过获取两个时相的遥感图像的纹理特征影像，也就相当于获得了监测区域在两个时相上的地物分布情况。

目标区域获取模块502，用于根据两时相遥感图像的灰度差异以及纹理特征影像的纹理差异，获得监测区域内灰度和纹理特征变化均大于预设阈值的区域，作为目标区域。

由上述内容可知，遥感图像的灰度可以在一定程度上反应真实地表的情况，单纯利用两时相的遥感图像的灰度差异可以直接识别出监测区域在两时相上的改变，只用灰度进行识别就是简单对比灰度变化，地表的细节也不能得到有效的凸显，因此这种方式是不准确的。而纹理特征影像则是对灰度信息的再加工，能够体现出图像的纹理。

由于遥感图像的纹理是对像元之间关系的一种客观描述，若将一幅遥感图像划分为多幅子图像，在每幅子图像都会有相似的像元关系出现，也就是在每幅图像中出现相同的纹理，从整幅遥感图像上能够发现有相同纹理出现多次，由此表明纹理具有自相似性；另一方面，同一地区可以有不同分辨率的遥感图像，分辨率会影响纹理的粗糙程度，对于每一种纹理，粗糙程度只是反映局部具有相同关系的像元对个数的变化，而像元间的关系没有改变，比如像元间的灰度差，相对位置等，所以纹理具有标度不变性。结合到本发明实施例的应用场景中，可以理解的是，出于获取未利用地的土壤污染风险区域的考虑，监测区域本身的地理位置会比较偏僻，并且地理环境偏向于自然、原始、少被人类意志所改造，因此这种区域的环境中地物分布一般是均匀的，纹理特征影像的纹理会比较自然，若根据两个时相的灰度共生矩阵获知纹理特征变化大，说明对应的区域很可能发生了不寻常的变化，采用上述方法能够快速地识别出这些变化大的地方，同时基于纹理特征影像的性质，即使两个时相的遥感图像分辨率不同，也不会影响后续获取目标区域的结果，进一步增强了本发明实施例的识别方法的适应性。

识别模块503，用于根据监测区域内的未利用地的范围，利用高分辨率影像对所述目标区域进行识别，获得未利用地上的土壤污染风险区域。

需要说明的是，本发明实施例在获得了监测区域内变化较大的区域，即目标区域后，结合监测区域内的未利用地获得目标区域中的未利用地，即未利用地中变化较大的区域，再利用高分辨率影像对未利用地中变化较大的区域进行目视解译和空间关系分析，最终获得未利用地上的土壤污染风险区域，本发明实施例克服对监测区域进行逐像元查找、缺乏针对性的弊端，极大地减少了工作量，提高了监测效率。并且，由于未利用地具有环境单一的特征，因此将未利用地反映至图像中会出现整体的颜色变化较小的特点，本发明实施例通过利用高分辨率影像对未利用地进行目视解释和空间关系分析时，可以快速地从图像中颜色变化大的区域中获得未利用地上的土壤污染风险区域，能够大大提高监测效率。

在上述实施例的基础上，目标区域获取模块具有包括：

灰度区域获取单元，用于根据两时相的遥感图像中对应位置像素点的灰度，获得监测区域内灰度差值绝对值大于预设阈值的区域，作为第一区域。

具体地，两个时相的遥感图像的大小均为M×N，将前一个时相的遥感图像中第m行第n列的像素点的灰度与后一个时相的遥感图像中的第m行第n列的像素点的灰度求差，将差值的绝对值作为对应位置的像素点的灰度差异，将该灰度差异与预设阈值进行比对，如果该灰度差异大于了预设阈值，就将该像素点对应的监测区域内的真实区域作为第一区域的一部分，这样通过遍历所有像素点，就获得了第一区域。显然，M、N、m和n均为正整数，并且m∈M，n∈N。

纹理区域获取单元，用于根据两时相的纹理特征影像中对应位置像素点的纹理特征值，获得监测区域内纹理特征值差值的绝对值大于预设阈值的区域，作为第二区域。

具体地，两个时相的纹理特征影像的大小均为M×N，将前一个时相的纹理特征影像中第m行第n列的像素点的纹理特征值与后一个时相的纹理特征影像中的第m行第n列的像素点的纹理特征值求差，将差值的绝对值作为对应位置的像素点的纹理差异，将该纹理差异与预设阈值进行比对，如果该纹理差异大于了预设阈值，就将该像素点对应的监测区域内的真实区域作为第二区域的一部分，这样通过遍历所有像素点，就获得了第二区域。

目标区域获取单元，用于将第一区域与第二区域的重合区域作为所述目标区域。

本发明实施例的目标区域获取模块通过对比遥感图像的灰度差异以及纹理特征影像的纹理差异，实现了获取目标区域的作用。并且，由于目标区域是指灰度异常区域与纹理异常区域的重合区域，因此在实际应用时，获取目标区域的方法还可以采用以下形式：

根据两时相的遥感图像中对应位置像素点的灰度，获得监测区域内灰度差值绝对值大于预设阈值的区域，作为第三区域；

根据两时相的纹理特征影像中对应第三区域的像素点的纹理特征值，获得第三区域中纹理特征值差值的绝对值大于预设阈值的区域，作为第四区域；

将所述第四区域作为目标区域。

本发明实施例区别于上一个实施例的地方在于，首先通过灰度差异对监测区域进行筛选，获得一个范围较小的区域，即第三区域，再根据两时相的纹理特征影像中对应第三区域的像素点的纹理特征值，进一步对第三区域进行筛选，获得第四区域，并将第四区域作为目标区域，相比上一个实施例缩小了识别面积，提供了识别效率。

在上述实施例的基础上，本发明实施例主要描述纹理特征影像获取模块，该模块具体包括：

预处理单元，用于确定两时相遥感图像的灰度等级，选定若干个特征量作为纹理统计量。

需要说明的是，为了减少计算量和计算时间，首先需要将两时相的遥感图像的灰度级统一压缩至一定的级数。由于灰度共生矩阵会随着选择的两个点的取值变化而变化，因此，在具体实施时，本发明实施例会预先确定两个点的角度和距离，以使得两个时相的灰度共生矩阵是相互对应的。

一般来说，如果图像的是由具有相似灰度值的像素块构成，则灰度共生矩阵的对角元素会有比较大的值；如果图像像素灰度值在局部有变化，那么偏离对角线的元素会有比较大的值。通常可以用一些标量来表征灰度共生矩阵的特征，即特征量，常用的特征量有：

ASM(能量，angular second moment)，也即每个矩阵元素的平方和，反映了图像灰度分布均匀程度和纹理粗细度。如果灰度共生矩阵的所有值均相等，则ASM值小；相反，如果其中一些值大而其它值小，则ASM值大。当共生矩阵中元素集中分布时，此时ASM值大。ASM值大表明一种较均一和规则变化的纹理模式。如果灰度共生矩阵中的值集中在某一块(比如对连续灰度值图像，值集中在对角线；对结构化的图像，值集中在偏离对角线的位置)，则ASM有较大值，若灰度共生矩阵中的值分布较均匀(如噪声严重的图像)，则ASM有较小的值。

对比度(contrast，CON)，直接反映了某个像素值及其领域像素值的亮度的对比情况。如果偏离对角线的元素有较大值，即图像亮度值变化很快，则CON会有较大取值，这也符合对比度的定义。对比度也反映了图像的清晰度和纹理沟纹深浅的程度。纹理沟纹越深，其对比度越大，视觉效果越清晰；反之，对比度小，则沟纹浅，效果模糊。灰度差即对比度大的象素对越多，这个值越大。灰度共生矩阵中远离对角线的元素值越大，CON越大。

逆差矩(inverse different moment，IDM)，反映图像纹理的同质性，度量图像纹理局部变化的多少。其值大则说明图像纹理的不同区域间缺少变化，局部非常均匀。如果灰度共生矩阵对角元素有较大值，IDM就会取较大的值。因此连续灰度的图像会有较大IDM值。

熵(entropy，ENT)，熵是图像所具有的信息量的度量，纹理信息也属于图像的信息，是一个随机性的度量，它表示了图像中纹理的非均匀程度或复杂程度。若灰度共生矩阵值分布均匀，也即图像近于随机或噪声很大，熵会有较大值。

自相关(correlation，COR)，反应了图像纹理的一致性。如果图像中有水平方向纹理，则水平方向矩阵的COR大于其余矩阵的COR值。它度量空间灰度共生矩阵元素在行或列方向上的相似程度，因此，相关值大小反映了图像中局部灰度相关性。当矩阵元素值均匀相等时，相关值就大；相反，如果矩阵像元值相差很大则相关值小。

根据灰度共生矩阵的若干个特征量，能够更准确地表达地表地物的分布，进一步地，根据两时相对应的灰度共生矩阵的特征量，也就可以准确地获得地标地物的变化情况，同时结合灰度变化程度，获得监测区域内灰度和纹理特征变化均大于预设阈值的区域，作为目标区域。通过上述步骤，本发明实施例可以初步找出可能出现监测区域内出现污染的区域，这种将灰度差异与灰度共生矩阵相结合的方式能够更准确地找出监测区域内环境发生变化的区域。

纹理特征值获取单元，用于利用滑动窗口对所述遥感图像截取子图像，对每次截取的子图像利用灰度共生矩阵法计算所述纹理统计量，并将所述纹理统计量作为对应子图像中心像素的纹理特征值。

具体地，设遥感图像的大小为M×N，滑动窗口的大小为L×L，其中L＜M，L＜N，利用滑动窗口由遥感图像的一个端点开始滑动，每次滑动将截取一个大小为L×L的子图像，然后在该子图像中计算0°、45°、90°和135°方向上的特征量，例如ASM、CON、IDM、ENT，对于每一种特征量，将该特征量4个方向的平均值作为结果，显然，每一个子图像的中心像素的纹理特征值的个数与特征量的种类一致。

特征影像获取单元，用于获取所有像素点的纹理特征值，对任意一个像素点，将该像素点的纹理特征值中的各特征量进行加权求和，作为该像素点的新灰度值，以获得纹理特征影像。

例如，在获得每个像素点的纹理特征值后，若纹理特征值包含ENT、COR、IDM以及CON，分别对ENT、COR、IDM以及CON赋予权重值，然后将每个特征量乘以权重后做总和，将求和后的值作为像素点的新灰度值，通过统计每个像素点的新灰度值，就可以获得纹理特征影像。具体权重值的设置可以通过预先测试一定数量的图像，判断各特征量对识别的影响大小，根据各特征量的影响设置权重值。

由于两个时相的遥感影像或多或少会存在拍摄环境不同的情况，如果直接以从地质局等单位直接获取的遥感影像进行识别，两个时相的遥感影像的灰度并不是基于同一个标准建立的，进而会在后续步骤中产生误差，因此，本发明实施例的识别装置还包括：预处理模块，用于对遥感图像进行预处理。具体地，该预处理模块包括：

Landast影像单元，用于获取覆盖监测区域的云量小于预设门限的两时相的Landsat影像。

需要说明的是，由于目前卫星的透云透雾能力还不完善，在选择遥感影像时需要尽量避免有云天气，可以预设门限值为10％，即只选取覆盖监测区域的云量小于10％的遥感图像，具体地界定监测区域的云量的方法已被本领域技术人员所熟知，在此不再赘述。Landsat影像是指美国NASA陆地探测卫星***获取的遥感图像。

校正单元，用于对两时相的Landsat影像分别进行辐射定标、大气校正以及几何校正处理，并叠加矢量图像，裁剪出监测区域在两个时相的原始遥感影像

需要说明的是，辐射定标是将传感器所得的测量值变化为绝对亮度或变换为与地表反射率、表面温度等物理量有关的相对值的处理过程，能够消除传感器本身的误差，确定传感器入口处的准确辐射值。辐射定标主要有三种定标方法，分别为实验室定标、机上/星上定标以及场地定标。大气校正将辐射亮度或者表面反射率转换为地表实际反射率，能够消除大气散射、吸收、反射引起的误差。现有的成熟技术主要有三种，分别为基于辐射传输模型的校正方法、基于简化辐射传输模型的黑暗像元法以及基于统计学模型的反射率反演。几何校正利用GCP纠正各种因素引起的几何变形，对影像进行地理坐标定位，获得真实坐标信息。主要包括：空间像元位置的变换和变换后像元亮度值的计算。几何校正用于纠正***和非***因素引起的几何畸变。对两时相的Landsat影像进行图像校正的流程顺序为：辐射定标——大气校正——几何校正，经过图像校正后的影像降低提取纹理、灰度特征变化区域时的误差，使相对稳定的同名地物的辐射值在不同时相的遥感图像上保持一致。

主成分分析单元，用于对两个时相的原始遥感影像分别进行主成分分析，提取前一定个数的主成分，构成对应时段的所述遥感影像。

需要说明的是，为减少后续处理计算量，对校正单元处理后的图像进行主成分分析，主成分分析用于将原来影像中的大部分信息集中到尽可能少的几个新主成分影像中，并使各个主成分包含的信息彼此不重叠，从而减少总数据量，消除冗余信息，突出有效信息。经验证，第一主成分对原始数据的贡献率均在85％以上，因此在一个可选实施例中选择第一主成分作为对应时段的遥感影像。

通过上述实施例获得未利用地上的土壤污染风险区域后，为了进一步判断疑似污染区域，分析出相关的污染源，本发明实施例的识别装置，还包括：

污染源识别模块，用于利用POI数据对所述土壤污染风险区域进行空间验证，判别污染风险点的污染源。

具体地，利用POI数据对土壤污染风险区域进行空间验证分析，结合Arcgis软件中的近邻分析工具，通过设置搜寻半径来查找最邻近的POI点位，确定污染区域是否为附近企业造成，并分析造成污染的相应污染源。

POI数据可以是对一个垃圾堆、一个商业住宅、一个工厂企业等信息的描述，它的主要内容包括：名称、类别、经纬度、地址名称，获取POI数据的手段多样，例如高德地图、报读地图的POI数据均可以开源获得，POI数据分类较细，并且可以根据实际情况定时更新。

以上所描述的装置实施例仅仅是示意性的，其中所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部模块来实现本实施例方案的目的。本领域普通技术人员在不付出创造性的劳动的情况下，即可以理解并实施。

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到各实施方式可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件。基于这样的理解，上述技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在计算机可读存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行各个实施例或者实施例的某些部分所述的方法。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (10)

1.一种未利用地上土壤污染风险区域的识别方法，其特征在于，包括：

根据监测区域在两个时相的遥感图像中每个像素点的灰度，获得反映各遥感图像的纹理特征的纹理特征影像；

根据两时相遥感图像的灰度差异以及纹理特征影像的纹理差异，获得监测区域内灰度和纹理特征变化均大于预设阈值的区域，作为目标区域；

根据监测区域内的未利用地的范围，利用高分辨率影像对所述目标区域进行识别，获得未利用地上的土壤污染风险区域。

2.如权利要求1所述的识别方法，其特征在于，所述根据两时相遥感图像的灰度差异以及两时相纹理特征影像的纹理差异，获得监测区域内灰度和纹理特征变化均大于预设阈值的区域，作为目标区域的步骤，具体包括：

根据两时相的遥感图像中对应位置像素点的灰度，获得监测区域内灰度差值绝对值大于预设阈值的区域，作为第一区域；

根据两时相的纹理特征影像中对应位置像素点的纹理特征值，获得监测区域内纹理特征值差值的绝对值大于预设阈值的区域，作为第二区域；

将第一区域与第二区域的重合区域作为所述目标区域。

3.如权利要求1所述的识别方法，其特征在于，所述获得反映各遥感图像的纹理特征的纹理特征影像的步骤之前，还包括获得所述监测区域在两个时相的遥感图像的步骤：

获取覆盖监测区域的云量小于预设门限的两时相的Landsat影像；

对所述两时相的Landsat影像分别进行辐射定标、大气校正以及几何校正处理，并叠加矢量图像，裁剪出监测区域在两个时相的原始遥感影像；

对两个时相的原始遥感影像分别进行主成分分析，提取前一定个数的主成分，构成对应时段的所述遥感影像。

4.如权利要求1所述的识别方法，其特征在于，所述根据监测区域内的未利用地的范围，利用高分辨率影像对所述目标区域进行识别之前还包括：获取监测区域内的未利用地的步骤：

根据遥感影像的光谱特征，结合野外实测资料和高分辨率影像确定真实地物类别；

根据真实地物在遥感影像上的色调、形状、纹理及空间位置特征，建立未利用地解译标志；

根据所述未利用地解译标志进行目视解译，确定监测区域内的未利用地的矢量边界。

5.如权利要求1所述的识别方法，其特征在于，所述获得未利用图上的土壤污染风险区域的步骤之后，还包括：

利用POI数据对所述土壤污染风险区域进行空间验证，判别污染风险点的污染源。

6.如权利要求1-5任意一项所述的识别方法，其特征在于，所述根据监测区域在两个时相的遥感图像中每个像素点的灰度，获得反映各遥感图像的纹理特征的纹理特征影像的步骤，具体包括：

确定所述两时相的遥感图像的灰度等级，选定若干个特征量作为纹理统计量；

利用滑动窗口对所述遥感图像截取子图像，对每次截取的子图像利用灰度共生矩阵法计算所述纹理统计量，并将所述纹理统计量作为对应子图像中心像素的纹理特征值；

获取所有像素点的纹理特征值，对任意一个像素点，将该像素点的纹理特征值中的各特征量进行加权求和，作为该像素点的新灰度值，以获得纹理特征影像。

7.如权利要求6所述的识别方法，其特征在于，所述纹理统计量包括以下特征量：能量、对比度、逆差矩以及熵。

8.一种未利用地上土壤污染风险区域的识别装置，其特征在于，包括：

纹理特征影像获取模块，用于根据监测区域在两个时相的遥感图像中每个像素点的灰度，获得反映各遥感图像的纹理特征的纹理特征影像；

目标区域获取模块，用于根据两时相遥感图像的灰度差异以及纹理特征影像的纹理差异，获得监测区域内灰度和纹理特征变化均大于预设阈值的区域，作为目标区域；

识别模块，用于根据监测区域内的未利用地的范围，利用高分辨率影像对所述目标区域进行识别，获得未利用地上的土壤污染风险区域。

9.如权利要求8所述的识别装置，其特征在于，所述目标区域获取模块具体用于：

根据两时相的遥感图像中对应位置像素点的灰度，获得监测区域内灰度差值绝对值大于预设阈值的区域，作为第一区域；

根据两时相的纹理特征影像中对应位置像素点的纹理特征值，获得监测区域内纹理特征值差值的绝对值大于预设阈值的区域，作为第二区域；

将第一区域与第二区域的重合区域作为所述目标区域。

10.如权利要求8所述的识别装置，其特征在于，还包括：

污染源识别模块，用于利用POI数据对所述土壤污染风险区域进行空间验证，判别污染风险点的污染源。