一种污染场地遥感探测与空间热区识别系统和方法与流程

本发明涉及一种污染场地遥感探测与空间热区识别系统和方法，属于城市规划及环境保护行业应用的遥感和地理信息技术应用领域，重点应用在污染场地管理工作中，服务于宏观监测污染场地空间分布，探测污染场地集中分布区域，在当前污染场地管理、测绘与地理信息和生态环境保护领域具有广泛而急迫的应用需求。

背景技术：

近年来污染场地管理行业在场地信息收集、污染调查以及法规建设方面已经有了长足的进步，但当前大量污染场地处于高风险状态，在土地用途转换以及二次开发利用过程中尚缺乏有效监管，严重威胁着人居环境安全和公众健康，环境隐患尤为突出。初步估计，我国各类工业污染场地至少以数十万计，多数分布在经济发达地区和老工业基地。据《中国环境年鉴》统计，仅1998～2009年我国关停并转迁企业总数累积就达到15.5万之多。自国家2010年启动实施高分专项工程以来，高分辨率对地观测系统技术迅速发展，国产高分影像被广泛应用于生态环境监测和评估等领域，利用高分卫星影像实时、高效和准确地探测污染场地的空间位置和集中分布区域，加强污染场地监测，在污染场地管理和生态环境保护领域具有广泛而急迫的应用需求。

中国专利公报公开显示，“污染场地”正成为专利授权的热点领域和新兴方向。如“一种石油类污染场地土壤的原位联合修复方法”(公开号cn105772499b、公开日2019年4月23日)，“一种污染土壤修复范围及边界的优化方法”(公开号cn106446435b、公开日20198年4月5日)，“一种重金属污染场地注浆修复剂及其制备方法”(公开号cn105602568b、公开日2018年12月25日)，“一种用于污染场地修复的太阳能通风系统装置”(公开号cn106001076b、公开日2018年5月11日)等一批专利取得授权。但是，这些发明专利主要集中于针对某种具体类型的污染场地修复工作或发明某种污染场地修复装置，较少有针对污染场地遥感识别和热区探测方面的相关专利，不能满足污染场地规划设计、现状评估的应用需求。因此，针对污染场地遥感探测方面的现有技术相对较少。

当前在污染场地遥感探测和热区识别中存在以下困难：(1)基于科学合理的分类系统，确定不同类型污染物边界；(2)通过遥感影像快速、准确和实时地区分场地污染物，获取污染物位置、面积、形状等信息(3)污染场地热区识别，分析污染场地热区空间分布；

为解决上述问题，本方法提出一种新的污染场地遥感图像探测方法，实现了针对污染场地集中分布区域的识别方法，且具有实用性强的特点，在当前污染场地现状监测、污染场地规划、生态环境修复等相关领域具有广阔的应用前景。

技术实现要素：

本发明技术解决问题：克服现有技术的不足，提供一种污染场地遥感识别和空间热区探测的系统和方法。该系统方法在污染场地监测、测绘与地理信息和生态环境保护领域具有急迫的应用需求和重大的应用价值，为识别污染场地空间分布提供科学依据，实现在大数据支持下的以污染场地管理为目标的可视化空间决策支持功能。

本发明所采用的技术方案为：一种污染场地遥感探测与空间热区识别方法，包括如下步骤：

步骤(1)，空间数据采集与预处理；

首先，获取gf-1号高分辨率遥感影像并对遥感影像进行预处理，包括辐射校正、大气校正、正射校正操作，作为污染场地遥感分类的数据基础；然后，针对研究区建立地物分类系统，作为污染场地遥感分类的分类依据；

步骤(2)，污染场地遥感探测，基于所述的分类系统与预处理后的遥感影像，采用面向对象方法进行污染场地遥感探测；

首先，建立遥感解译标志库，为遥感解译提供样本参考和先验知识，解译标志包括光谱、形状、大小、纹理；

然后，基于解译标志，确定分类特征，构建分类规则；

最后，基于分类规则，使用模糊函数确定地物类型的隶属度，进行地物分类，提取不同类型污染场地，从而获得区域内污染场地分类结果；

步骤(3)，污染场地热区探测；

首先，基于所述的污染场地分类结果，利用gis空间分析方法，计算污染场地空间分布密度；

然后，构建污染场地空间分布热度模型，并计算污染场地空间分布热度指数；

最后，对污染场地空间分布热度指数等级进行划分，探寻污染场地空间分布热区。

所述空间数据采集与预处理包括：

(1.1)所述的gf-1号高分辨率遥感影像，其多光谱空间分辨率为16m，全色波段的空间分辨率最高为2m，gf-1号遥感影像下载自地理国情监测云平台；

(1.2)遥感影像的预处理，基于envi5.1软件进行；首先，对多光谱波段，进行辐射定标、大气校正、正射校正预处理；对全色波段，进行辐射定标和大气校正；然后，将多光谱与全色波段进行图像融合，获得处理后的遥感图像，便于模块二进行地物识别；

(1.3)所述的地物分类系统，包括7个大类：植被、水体、裸地、建筑、道路、工矿区污染场地和其他，其中工矿区污染场地可进一步细分为4类：排土排矸场、尾矿库、煤矸石堆场、采煤坑。

所述污染场地遥感探测步骤具体如下：

(2.1)所述的各地物类型遥感解译类型，包括裸地、建筑、道路、排土排矸场、尾矿库、煤矸石堆放场、采煤坑、植被、水体；

(2.2)所述的分类特征，包括亮度、比率、各方向灰度均值、面积、形状指数、长宽比、距离、矩形相似度、归一化水体指数、归一化建筑指数和归一化植被指数；

(2.3)所述的具体分类规则，见如下规则(1)：

其中，b为亮度、r为比率、g为各方向灰度均值、a为面积、s为形状指数、lw为长宽比、rf为矩形相似度、ndwi为归一化水体指数、ndbi为归一化建筑指数、ndvi为归一化植被指数。

所述污染场地热区识别具体如下：

(3.1)根据所述的污染场地遥感分类数据，在arcgis中，利用arctoolbox里的featuretopoint工具，将污染场地矢量面状要素转为点状要素，得到污染场地的点位空间分布数据；然后，基于密度分析方法，分析污染场地分布特征，使用gis空间分析中的基于点要素的核密度分析方法，计算污染场地的周围邻域中的空间分布密度，从而获得污染场地分布密度指数数据；

(3.2)基于污染场地面积指数数据与污染场地分布密度指数数据，构建污染场地空间分布热度模型，计算污染场地空间分布热度指数，污染场地空间分布热度指数计算公式(2)如下：

式中，di为第i个栅格单元的污染场地分布密度值；dmin为所有di的最小值，dmax为所有di的最大值；ri为第i个栅格单元的污染场地空间分布热度指数；

(3.3)基于污染场地空间分布热度指数数据，在gis中，利用基于自然断点的数据标准分类方法，识别污染场地空间分布热区并划分等级，在arcgis中，对污染场地空间分布热度指数数据进行自然断点分级，共分为4类，见如下规则(2)：

其中，ri为第i个栅格单元的污染场地空间分布热度指数；

从而形成污染场地空间分布高聚集区、中聚集区、低聚集区和无分布区，完成污染场地空间分布热区的识别。

另一方面，本发明还提出一种污染场地遥感探测与空间热区识别系统，包括：空间数据采集与预处理模块、污染场地遥感识别模块和污染场地热区探测模块，其中：

空间数据采集与预处理模块，首先，获取gf-1号高分辨率遥感影像并对遥感影像进行预处理，包括辐射校正、大气校正、正射校正操作，作为污染场地遥感分类的数据基础；然后，针对研究区建立地物分类系统，作为污染场地遥感分类的分类依据；

污染场地遥感探测模块，基于所述的分类系统与预处理后的遥感影像，采用面向对象方法进行污染场地遥感探测；首先，建立遥感解译标志库，为遥感解译提供样本参考和先验知识，解译标志包括光谱、形状、大小、纹理；然后，基于解译标志，确定分类特征，构建分类规则；最后，基于分类规则，使用模糊函数确定地物类型的隶属度，进行地物分类，提取不同类型污染场地，从而获得区域内污染场地分类结果；

污染场地热区探测模块，首先，基于所述的污染场地分类结果，利用gis空间分析方法，计算污染场地空间分布密度；然后，构建污染场地空间分布热度模型，并计算污染场地空间分布热度指数；最后，对污染场地空间分布热度指数等级进行划分，探寻污染场地空间分布热区。

进一步的，所述空间数据采集与预处理模块具体实现如下：

(1)所述的gf-1号高分辨率遥感影像，其多光谱空间分辨率为16m，全色波段的空间分辨率最高为2m，gf-1号遥感影像下载自地理国情监测云平台；

(2)遥感影像的预处理，基于envi5.1软件进行；首先，对多光谱波段，进行辐射定标、大气校正、正射校正预处理；对全色波段，进行辐射定标和大气校正；然后，将多光谱与全色波段进行图像融合，获得处理后的遥感图像，便于模块二进行地物识别；

(3)所述的地物分类系统，包括7个大类：植被、水体、裸地、建筑、道路、工矿区污染场地和其他，其中工矿区污染场地可进一步细分为4类：排土排矸场、尾矿库、煤矸石堆场、采煤坑。

进一步的，污染场地遥感探测模块具体实现如下：

(1)所述的各地物类型遥感解译标志，包括裸地、建筑、道路、排土排矸场、尾矿库、煤矸石堆放场、采煤坑、植被、水体；

(2)所述的分类特征，包括亮度、比率、各方向灰度均值、面积、形状指数、长宽比、距离、矩形相似度、归一化水体指数、归一化建筑指数和归一化植被指数；

(3)所述的具体分类规则，见规则(1)：

其中，b为亮度、r为比率、g为各方向灰度均值、a为面积、s为形状指数、lw为长宽比、rf为矩形相似度、ndwi为归一化水体指数、ndbi为归一化建筑指数、ndvi为归一化植被指数。

进一步的，所述污染场地热区识别模块具体实现如下：

(1)根据所述的污染场地遥感分类数据，在arcgis中，利用arctoolbox里的featuretopoint工具，将污染场地矢量面状要素转为点状要素，得到污染场地的点位空间分布数据；然后，基于密度分析方法，分析污染场地分布特征，使用gis空间分析中的基于点要素的核密度分析方法，计算污染场地的周围邻域中的空间分布密度，从而获得污染场地分布密度指数数据；

(2)基于污染场地面积指数数据与污染场地分布密度指数数据，构建污染场地空间分布热度模型，计算污染场地空间分布热度指数，污染场地空间分布热度指数计算公式如下：

式中，di为第i个栅格单元的污染场地分布密度值；dmin为所有di的最小值，dmax为所有di的最大值；ri为第i个栅格单元的污染场地空间分布热度指数；

(3)基于污染场地空间分布热度指数数据，在gis中，利用基于自然断点的数据标准分类方法，识别污染场地空间分布热区并划分等级，在arcgis中，对污染场地空间分布热度指数数据进行自然断点分级，共分为4类，具体划分规则(2)如下：

其中，ri为第i个栅格单元的污染场地空间分布热度指数；

从而形成污染场地空间分布高聚集区、中聚集区、低聚集区和无分布区，完成污染场地空间分布热区的识别。

本发明与现有技术相比的优点在于：

(1)本发明针对污染场地监测遇到的突出问题，提出了一种污染场地遥感探测与空间热区识别方法。针对中国快速的工业化和城市化进程导致污染场地数量和规模显著增加等现象，提出了服务于基于遥感图像的污染场地识别完整的方法体系，研发的模型系统具有直观明确、可操作性强、具有较强的可视化程度和辅助决策支持功能，对于环境监管部门具有重要的意义。

(2)本发明提出了一种污染场地遥感探测方法，使用高分卫星遥感影像，应用计算机自动分类方式，区分不同类型污染物，实现快速和准确的获取污染场地位置、面积、形状等信息，为污染场地管控方面的规划者或决策者提供规划参考依据，具有广泛的行业应用需求。

(3)本发明提出了一种利用污染场地空间分布信息、基础地理信息，通过核密度分析方法，探测污染场地分布热区的方法，可以快速获取和检测哪些区域污染场地分布较多，哪些区域污染场地分布较为聚集。解决了污染场地空间分布情况与变化情况等一系列污染场地检测与分析识别问题。

附图说明

图1为本发明污染场地遥感探测与空间热区识别系统结构图；

图2为高分一号多光谱遥感影像预处理技术流程图；

图3为污染场地遥感探测模块技术流程图；

图4为污染场地热区探测模块技术流程图；

图5为污染场地遥感探测结果示意图；

图6为污染场地热区分布示意图。

具体实施方式

下面结合附图及实施例对本发明进行详细说明。

如图1所示，本发明所述的系统包括空间数据采集与预处理模块、污染场地遥感探测模块、污染场地热区识别模块。

如图2所示，本发明中空间数据采集与预处理模块用于进行空间数据采集与预处理，包括如下步骤：

步骤1-1，对下载的高分一号原始影像进行辐射校正操作；首先，借助envi5.1软件，打开原始遥感影像；然后，借助radiometriccalibration工具，选择需要进行辐射校正的影像，设置定标形式(calibrationtype)为默认的radiance，将outputinterleave改为bil顺序，输出类型默认浮点型，将scalefactor改为0.10，选择输出路径，运行，得到辐射定标结果，辐射定标结果以dat文件格式存储，以此作为步骤1-2进行大气校正的输入数据；

步骤1-2，对经辐射定标后的高分一号影像进行大气校正；首先，获取影像对应区域平均海拔高程，通过elevation(gmted2010)，打开envi自带全球900米分辨率的dem数据，打开computestatistics输入文件对话框，选择gmted2010.jp2数据。单击statssubset按钮，打开selectstatisticssubset对话框，选择统计区域对应的高分一号图像。在统计面板中，选择basicstats基本统计，得到平均海拔高程(单位：“米”)；然后，进行大气校正操作；借助flaashatmosphericcorrection工具，设置相关参数，选择上一步生成的辐射定标结果作为输入辐射定标数据，选择将单一比例参数应用到所有波段，将singlescalefactor默认设置为1，sensortype设置为unknown-msi，sensoraltitude设置为传感器高度645km，groundelevation设置为计算的影像对应区域平均高程，需要将米转为千米，pixelsize设置为16m，flightdate设置为影像获取时间，可通过查看原始影像元数据viewmetadata获取，atmosphericmodel大气模型设置根据影像实际情况选择，aerosolretrieval气溶胶反演算法，设置为none，设置大气校正结果输出位置；最后，在multispectralsettings面板进行多光谱参数设置，filterfunctionfile设置为对应的波谱响应函数，其他参数保持默认，运行大气校正模块，获得大气校正结果，作为步骤1-3的输入数据；

步骤1-3，对高分一号影像进行正射校正；首先，利用rpcorthorectificationworkflow，打开正射校正工作流，选择大气校正后的影像；然后，进入orthorectification参数设置，图像像元大小设置为16，重采样方式设置为三次卷积，其他默认；最后，设置正射校正后的影像输出路径，运行正射校正模块，获得正射校正后的影像；

步骤1-4，对高分一号影像进行几何校正；首先，利用imageregistrationworkflow，选择基准影像为完全覆盖待处理区域的landsat8影像，待校正影像为正射校正后的影像；然后，设置影像配准主要参数；最后，设置输出路径；运行几何校正模块，获得几何校正输出结果，完成原始高分一号影像的预处理，经预处理后的影像作为步骤2的输入数据，进行污染场地的遥感探测；

步骤1-5，建立进行污染场地遥感探测的地物分类系统，将典型地物分为7个大类：植被、水体、裸地、建筑、道路、工矿区污染场地和其他，其中工矿区污染场地可进一步细分为5类：排土排矸场(采矿中剥离的表土和围岩)、尾矿(矿物含量较低而无法用于生产的废弃物堆积)、煤矸石堆场、采煤坑等，作为步骤2污染场地遥感探测的分类基础。

如图3所示，本发明中污染场地遥感探测模块用于进行污染场地遥感探测，具体实现如下：

步骤2-1，建立遥感解译标志，为遥感解译提供样本参考和先验知识，解译标志包括光谱、形状、大小、纹理等，根据区域地表状况，分为植被、水体、裸地、建筑、道路、排土排矸场、尾矿库、煤矸石堆场和采煤坑等地物类型。

步骤2-2，基于步骤2-1确定的解释标志，确定分类特征，构建分类规则，分类特征具体包括亮度(brightness)、比率(ratio)、各方向灰度均值(glcmmeanalldir)、面积(area)、形状指数(shapeindex)、长宽比(length/width)、距离(distanceto)、矩形相似度(rectangularfit)、归一化水体指数(ndwi)、归一化建筑指数(ndbi)和归一化植被指数(ndvi)等，具体见规则(1)，如下：

其中，裸地分类特征为亮度大于2000，比率小于0.3，各方向灰度均值介于4.8-9.3之间；建筑分类特征为归一化建筑指数小于0.18；道路分类特征为长度大于4.6，长宽比大于2.5；排土排矸场分类特征为面积大于82，形状指数大于0.9；尾矿库分类特征为矩形相似度介于0.34-0.67之间；煤矸石堆放场分类特征为长度大于2031，面积大于82；采煤坑分类特征为长度介于756-898之间；植被分类特征为归一化植被指数大于0.2；水体特征为归一化水体指数小于-0.045；

步骤2-3，基于步骤2-2获得的预处理后的gf-1号遥感影像，对影像进行面向对象分割；

首先，基于易康软件，选择实验影像，选择加载多光谱的蓝、绿、红以及近红外四个波段信息，完成图像加载；

然后，利用processtree工具，打开进程树窗口，添加新的分类树节点，将name修改为segmentation，其他参数默认；

最后，点击segmentation，新建rule，建立insertchild，name设置为“分割”，算法选择设置为多尺度分割算法multiresolutionsegmentation，imageobjectdomain设置为像素水平，根据分类区域地物种类及分布情况，将scaleparameter设置为50，levelname设置为level1，形状因子设置为0.5，紧致度设置为0.5；

步骤2-4，基于步骤2-2确定的分类规则和步骤2-3获得的分割后的影像，对典型地物进行面向对象分类；

首先，基于易康软件，打开classhierarchy窗体，调出类层次窗体，添加地物类型，设置类别名称，基于规则2的分类规则，设置contained，设置分类特征和模糊函数；

然后，通过processtree中的“分割”节点，添加新的分类树节点，name设置为地物类型名称，algorithm设置为classification，imageobjectdomain设置为imageobjectlevel；最后，将activeclasses设置为在classhierarchy窗体已经添加过的相应地物类型；

最后，执行分类操作，得到分类结果，将工矿业污染场地进行提取，即可获得污染场地空间分布区域，作为污染场地热区识别模块的输入数据；

如图4所示，本发明中污染场地热区识别模块具体实现如下：

步骤3-1，基于步骤2-4得到的地物分类图，提取污染场地类型地物要素，并计算污染场地分布密度；

首先，基于arcgis，利用featuretopoint工具将污染场地面状要素转换成点要素，形成污染场地点位矢量数据图层；

然后，基于核密度分析方法，对污染场地点要素空间分布密度进行分析。在spatialanalyst支持下，选择kerneraldensity工具，依次输入各个参数：其中，输入要素设置为前一步骤生成的污染场地点位矢量数据图层；populationfield设置为“none”，即每个污染场地点状要素只计数一次；输出栅格单元尺寸设置为10，搜索半径设置为300，面积单位设置为“square_kilometers”，即平方公里。

最后，运行工具，得到栅格格式的污染场地密度分布数据，作为步骤3-2的输入数据；

步骤3-2，使用步骤3-1获得的污染场地分布密度栅格，基用公式(1)计算污染场地空间分布热度指数；

式中，di为第i个栅格单元的污染场地分布密度值；dmin为所有di的最小值，dmax为所有di的最大值；ri为第i个栅格单元的污染场地空间分布热度指数；

基于步骤3-1得到的污染场地分布密度栅格，打开数据属性信息，并在属性窗口中的classificationstatistics信息中查得栅格单元值的最大值dmax和最小值dmin；然后，在mapalgebra工具集下，选择并打开rastercalculator工具，基于公式(1)，输入栅格计算表达式“(d–dmin)/(dmax–dmin)”，其中，d为污染场地分布密度栅格，输出的结果栅格即为污染场地空间分布热度指数栅格数据，作为步骤3-3的输入数据；

步骤3-3，进行污染场地空间分布热度等级的划分。基于步骤3-3生成的污染场地空间分布热度指数栅格数据，利用规则(2)如下：

即将栅格单元的污染场地空间分布热度指数介于0.47-1之间栅格单元作为污染场地空间分布高聚集区，将栅格单元的污染场地空间分布热度指数介于0.24-0.47之间栅格单元作为污染场地空间分布中聚集区，将栅格单元的污染场地空间分布热度指数介于0.11-0.24之间作为污染场地空间分布低聚集区，将栅格单元的污染场地空间分布热度指数介于0-0.11之间作为污染场地空间分布无分布区，在arcgis软件中，基于自然断点分级的数据标准分类方法对前述的污染场地空间分布热度指数进行分级，共分为4类，形成污染场地的高聚集区、中聚集区、低聚集区和无分布区，完成污染场地空间分布热区的识别。

如图5所示，为污染场地遥感探测结果示意图，即步骤2-4的结果。

基于步骤2-2确定的分类规则和步骤2-3获得的分割后的影像，对典型地物进行面向对象分类，典型地物类型包括7个大类：植被、水体、裸地、建筑、道路、工矿区污染场地和其他，其中工矿区污染场地可进一步细分为4类：排土排矸场、尾矿库、煤矸石堆场、采煤坑等。

如图6所示，为污染场地热区分布示意图，即步骤3-3的结果。

基于步骤3-3生成的污染场地空间分布热度指数栅格数据，利用规则(2)，在arcgis软件中，基于自然断点分级的数据标准分类方法对前述的污染场地空间分布热度指数进行分级，共分为4类，图中颜色绿色、黄色和红色表示的区域依次为污染场地的低聚集区、中聚集区和高聚集区。

以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本申请专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。