一种基于时空交互关系的土壤镉风险预测方法

1.本发明属于土壤镉(cd)风险预测领域，具体涉及一种基于时空交互关系的土壤cd风险预测方法。


背景技术：

2.由于污染企业与农用地土壤污染之间的源汇理论关系，越来越多的研究者开始研究土壤重金属与污染企业之间的关系，进而为土壤重金属污染溯源提供证据链。然而，传统的源解析方法，如多元和地统计分析，主成分分析、正定矩阵分解等，在一定程度上都忽略了土壤重金属的时空污染特征，部分地区有可能得出相反的结果，难以有效和精准的指导土壤重金属管理。此外，土壤重金属污染往往与人为活动密切相关，尤其是污染企业，使用传统的源解析方法没有提供源汇之间可靠的时空耦合关系，对于污染的时空变异性也不能很好的解决，这使得目前由污染企业变迁引起的土壤重金属污染防控治理工作变得相当艰难。
3.目前，相继出台了一些关于土壤重金属污染分区管理的技术规定或标准，但这些分区管理办法仅仅考虑到当前土壤重金属污染水平，而忽略了土壤重金属污染的未来风险，比如由于污染企业引起的土壤污染风险。无论如何，基于未来风险的土壤污染管理可以作为未来土壤管理的重要参考，模拟土壤cd未来风险概率区的动态变化特征对于当前和未来土壤cd污染管理的有效决策具有重要意义。但是，如何实现土壤重金属污染的未来风险预测，是目前亟待解决的技术问题。


技术实现要素：

4.本发明的目的在于解决现有技术中存在的问题，并提供一种基于时空交互关系的土壤cd风险预测方法。
5.本发明所采用的具体技术方案如下：
6.一种基于时空交互关系的土壤镉风险预测方法，其步骤如下：
7.s1、获取目标预测区域内多个不同历史时期与土壤cd污染相关的数据集，并通过数据的提取和交叉验证，得到不同历史时期目标预测区域内的土壤cd污染企业空间点位分布数据，以及不同历史时期目标预测区域内各采样点的土壤cd浓度和土壤ph；
8.s2、针对每一个历史时期，使用核密度法对该历史时期目标预测区域内的土壤cd污染企业空间点位分布数据进行核密度分析，得到污染企业密度值，同时基于该历史时期目标预测区域内各采样点的土壤cd浓度和土壤ph，通过插值方法得到目标预测区域内的土壤cd浓度空间分布和土壤ph空间分布，结合不同ph范围对应的土壤cd浓度风险筛选值，确定目标预测区域内不同位置的土壤cd风险级别；
9.s3、针对每一个历史时期，对目标预测区域进行网格化后，统计各网格内的土壤cd风险级别与污染企业密度值，再对该历史时期土壤cd风险级别与污染企业密度值进行双变量局部空间自相关分析，得到土壤cd与污染企业的时空交互关系；
10.s4、针对每一个历史时期，根据对应的所述土壤cd风险级别、受体脆弱性和所述时空交互关系，按照预设的风险区分级标准进行目标预测区域内不同等级土壤cd风险区的识别；
11.s5、获取目标预测区域内每一个历史时期对应的土壤cd风险演变影响因素，结合各历史时期的土壤cd风险区识别结果，运用斑块生成土地利用模拟模型(patch-generating land use simulation model,plus)实现对目标预测区域内未来土壤cd风险区的预测。
12.作为优选，所述s1中，所述与土壤cd污染相关的数据集包括与土壤cd污染相关的涉土企业数据、poi数据、历史工业企业数据和多时期高分辨率遥感影像，以及不同时期土壤采样点数据；其中，涉土企业数据包括与土壤cd污染相关的企业的名称及其所属地市，poi数据包括企业名称及其经纬度信息；历史工业企业数据包括企业名称和生产活动年份；土壤采样点数据包括不同采样点处的土壤cd浓度和土壤ph数据。
13.进一步的，所述s1中，所述土壤cd污染企业空间点位分布数据的获取方法为：
14.s11、将所述与土壤cd污染相关的涉土企业数据中的非结构化数据转换为结构化数据，从而得到与土壤cd污染相关的涉土企业名称；对所述涉土企业名称、以及历史工业企业数据和poi数据中的企业名称分别进行分词处理，提取企业名称中不同层级的实体，随后进行相似度匹配，将企业名称中各层级实体完全匹配的涉土企业作为污染企业，将历史工业企业数据中该污染企业的生产活动年份和poi数据中的地理位置信息关联至该污染企业中，得到含地理位置和生产活动年份的污染企业数据。
15.s12、从多时期高分辨率遥感影像信息中提取每个污染企业及其周边场地在各生产活动年份的高分辨率遥感影像，然后使用基于深度学习的图像语义分割模型，对每个污染企业对应的高分辨率遥感影像数据进行建筑物特征提取，以判断各生产活动年份中影像区域内是否存在建筑或者企业工厂，实现不同年份污染企业空间点位分布的遥感验证，剔除未通过遥感验证的污染企业数据，按生产年份对剩余污染企业数据进行划分，得到不同历史时期目标预测区域内的土壤cd污染企业空间点位分布数据。
16.进一步的，所述s11中，企业名称使用分词引擎jieba进行分词处理，被分为行政区划、字号、行业、组织形式四个层级的实体，随后通过计算编辑距离进行相似度匹配，若两个企业名称中四个层级的实体完全匹配，则视为这两个企业名称匹配。
17.进一步的，所述s12中，所述基于深度学习的图像语义分割模型为u-net卷积神经网络模型；在对高分辨率遥感影像进行建筑物特征提取时，根据不同历史时期污染企业的经纬度及其对应遥感影像上的信息获取该点位周边相应的遥感影像；对获取的遥感影像进行数据预处理，生成相应的标注数据并运用u-net卷积神经网络模型进行训练建模；再利用训练好的u-net卷积神经网络模型对高分辨率遥感影像进行图像分割，返还的结果图中的每个像元有建筑与非建筑两类；将属于同一生产年份的返还的结果图与污染企业数据进行比较，若污染企业所在位置在返还的结果图中存在建筑，则说明该点位的污染企业真实存在；若污染企业所在位置在返还的结果图中只存在非建筑，则说明该点位的污染企业并不存在，将该污染企业信息进行二次审查以判断该点位是否存在污染企业。
18.作为优选，所述s3的具体方法如下：针对每一个历史时期，对所述土壤cd风险级别进行划分后，根据目标预测区域的边界，得到其最小外接矩形所围成的范围；然后以该最小
外接矩形的某个顶点开始，生成格网和对应的网格点；运用网格点提取不同历史时期污染企业核密度和土壤cd风险分级值，分别代表边界中不同历史时期污染企业的聚集程度和土壤cd风险水平；将带有土壤cd风险级别与污染企业密度值的网格点进行双变量局部莫兰分析，得到二者的时空交互关系，具体公式如下：
[0019][0020]
其中和分别是变量a和b在网格i和网格j上的不同历史时期土壤cd风险级别和污染企业密度，变量a和b分别代表网格内不同历史时期土壤cd风险级别与污染企业密度；w
ij
为网格i和网格j的空间权重矩阵，根据网格i和网格j之间的欧式距离权重所得；i
ab
表示网格i处的a属性与网格j处b属性的局部莫兰指数；当i
ab
显著为正时，则网格i处的土壤cd风险级别与网格j处的污染企业密度具有显著的局部正相关关系；当i
ab
显著为负时，则认为网格i处的土壤cd风险级别与网格j处的污染企业密度具有显著的局部负相关关系；当i
ab
不显著时，则认为网格i处的土壤cd风险级别与网格j处的污染企业密度无明显的时空交互关系。
[0021]
作为优选，所述s4中，在预设的风险区分级标准中，土壤cd风险区被分为高、中、低三种风险控制区和风险不确定区。
[0022]
作为优选，所述受体脆弱性为人口密度。
[0023]
作为优选，所述s5中，所述斑块生成土地利用模拟模型通过提取历史时期之间各等级土壤cd风险区的变化特征，并采用随机森林算法逐一对各土壤cd风险演变影响因素进行挖掘，获取各等级的土壤cd风险区的变化发生概率以及各土壤cd风险演变影响因素对变化特征的贡献，最终在变化发生概率的约束下，基于最新的土壤cd风险演变影响因素，实现对目标预测区域内未来土壤cd风险区的预测。
[0024]
作为优选，所述斑块生成土地利用模拟模型实现对目标预测区域内未来土壤cd风险区的预测的具体步骤如下：首先，将两个历史时期的各等级土壤cd风险区进行叠加分析，从中提取每种等级土壤cd风险区的变化特征；其次，基于随机森林模型建立每种等级土壤cd风险区发生变化的发生概率，并计算各土壤cd风险演变影响因素对每种土地利用变化的贡献率，发生概率用于估计网格单元中每种等级土壤cd风险区的动态变化趋势；最后，通过整合时间序列的土壤cd风险区数据和模型产生的组合概率来预测未来等级土壤cd风险区，其中通过以下公式估计网格单元p未来被k等级的土壤cd风险区占用的组合概率：
[0025][0026]
其中表示网格单元p在迭代时间t从原始等级的土壤cd风险区转换到目标k等级的土壤cd风险区的组合概率；p
p,k
表示网格单元p上目标k等级的土壤cd风险区出现的概率；表示在迭代时间t对应网格单元p上目标k等级的土壤cd风险区的邻域效应；表示在迭代时间t上目标k等级的土壤cd风险区的惯性系数；scc→k表示从原始c等级的土壤cd风险区到目标k等级的土壤cd风险区的转换成本；斑块生成土地利用模拟模型使用轮盘赌选择机制来确定哪种土地利用类型将占据网格单元。
[0027]
本发明相对于现有技术而言，具有以下有益效果：
[0028]
本发明基于公开的涉土企业信息，同时匹配当地政府调查并公开的涉污企业数据，获取与污染相关的涉土企业数据，后续无需对企业数据进行分类以及判别是否该企业是污染企业，可直接根据目标预测区域所在地政府公开的污染企业数据进行数据清洗与建模；此外，本发明根据不同时期污染企业的核密度分析得到密度值，并与对应时期的土壤cd风险级别进行双变量局部莫兰分析，分析二者的时空交互关系，从而使得离散的cd污染点数据与污染企业点数据能够较为准确的探索时空交互关系；最后通过不同时期土壤cd风险分区特征及其风险分区演变的影响因素，运用斑块生成土地利用模拟模型对当前不同时期的土壤cd风险分区进行建模，并预测未来土壤cd风险分区特征，扩展了现有土壤cd分区的方法和思路，对进一步管理与预防土壤cd污染具有重要的理论与实际意义，并存在推广应用的价值。
附图说明
[0029]
图1是实施例中的部分结果展示，其中(a)2002年土壤cd风险水平分区，(b)2002年人口密度的分区，(c)2002年污染企业密度分布，(d)2002年土壤cd与污染企业的局部莫兰空间交互关系；
[0030]
图2是实施例中的部分结果展示，其中(a)2012年土壤cd风险水平分区，(b)2012年人口密度的分区，(c)2012年污染企业密度分布，(d)2012年土壤cd与污染企业的局部莫兰空间交互关系；
[0031]
图3是实施例中的部分结果展示，其中(a)2002年土壤cd风险管理分区示意图，(b)2012年土壤cd风险管理分区示意图，(c)斑块生成土地利用模拟模型预测的2012年土壤cd风险管理分区示意图，(d)斑块生成土地利用模拟模型预测的2022年土壤cd风险管理分区示意图；
[0032]
图4是实施例中斑块生成土地利用模拟模型基于2002年风险分区预测2012年风险分区的模型精度。
具体实施方式
[0033]
下面结合附图和具体实施方式对本发明做进一步阐述和说明。本发明中各个实施方式的技术特征在没有相互冲突的前提下，均可进行相应组合。
[0034]
本发明是通过一种基于多源数据的相似度匹配算法进行污染企业空间识别，主要基于多源企业数据，建立相似度计算方法进行数据匹配，同时采用u-net模型建立污染企业识别的遥感验证，并通过得到的不同时期的污染企业数据与土壤cd污染数据进行双变量局部莫兰指数分析来构建污染企业与土壤cd风险分级之间的时空交互关系，同时运用源汇理论和时空交互关系来确定不同时期土壤cd风险分区特征，最终运用斑块生成土地利用模拟模型对未来土壤cd风险分区进行预测，进而对土壤cd污染管理分区研究工作起到一个指向性的作用。
[0035]
在本发明的一个较佳实施例中，提供了一种基于时空交互关系的土壤cd风险预测方法，其具体实现步骤如下：
[0036]
s1、获取目标预测区域内多个不同历史时期与土壤cd污染相关的数据集，并通过数据的提取和交叉验证，得到不同历史时期目标预测区域内的土壤cd污染企业空间点位分
布数据，以及不同历史时期目标预测区域内各采样点的土壤cd浓度和土壤ph。
[0037]
在本实施例中，上述s1步骤的具体实现方式如下：
[0038]
1)数据获取：
[0039]
获取与土壤cd污染相关的数据集，包括：与土壤cd污染相关的涉土企业数据、poi数据、历史工业企业数据和多时期高分辨率遥感影像，以及不同时期土壤采样点数据；其中，涉土企业数据包括与土壤cd污染相关的企业的名称及其所属地市，高德poi数据包括企业名称及其经纬度信息；历史工业企业数据包括企业名称和生产活动年份；土壤采样点数据包括不同采样点处的土壤cd浓度和土壤ph数据。
[0040]
本实施例中，与污染相关的涉土企业数据、poi数据、历史工业企业数据和多时期高分辨率遥感影像可使用基于http协议与公开网页端进行获取，poi数据采用高德poi数据；土壤采样点由田间采样获取，可从土壤普查数据中查询或者其他途径获取，每个采样点应当含有该采样点位置的土壤cd浓度和土壤ph。引入土壤ph的原因是，土壤中的cd风险是与ph密切相关的，因此需要结合土壤cd浓度和土壤ph来综合判断土壤cd风险区。另外，高德poi数据和历史工业企业数据应当尽量覆盖目标预测区域内的企业，以利于后续尽可能匹配到更多的污染企业数据。
[0041]
2)数据预处理：
[0042]
将与土壤cd污染相关的涉土企业数据中的非结构化数据转换为结构化数据，从而得到与土壤cd污染相关的涉土企业名称；对得到的涉土企业名称、以及历史工业企业数据和poi数据中的企业名称分别进行分词处理，提取企业名称中不同层级的实体，随后进行相似度匹配，将企业名称中各层级实体完全匹配的涉土企业作为污染企业，将历史工业企业数据中该污染企业的生产活动年份和poi数据中的地理位置信息关联至该污染企业中，得到含地理位置和生产活动年份的污染企业数据。
[0043]
在本实施例中，与污染相关的涉土企业中以非结构化文档数据pdf、word或图片格式存在信息，使用光学字符识别(optical character recognition，ocr)技术将其识别为基于office open xml标准的压缩文件格式，从而转换为结构化数据。
[0044]
在本实施例中，企业名称可使用分词引擎jieba进行分词处理，被分为行政区划、字号、行业、组织形式四个层级的实体，随后通过计算编辑距离进行相似度匹配，若两个企业名称中四个层级的实体完全匹配，则视为这两个企业名称匹配。jieba分词是基于词典的分词方法，为了得到较好的切分效果，需要使用自定义词典的方式来对未登录词进行更好的识别，通过加入了行政区划名和通名特征词到自定义词典，使得对名称的切分更为准确。本实施例中相似度匹配通过计算两个字符的编辑距离，从而得到字符串之间的差异性，差异越小其相似度越高。编辑距离是是指两个字符串之间，由一个转成另一个所需的最少编辑操作次数。如果它们的编辑距离越小，说明它们差异性越小，相似度越高。编辑距离允许的字符操作包括:a.删除一个字符，b.插入一个字符，c.修改一个字符。其核心思想是对一个字符串内部的单个字符进行最少的编辑操作后得到另一个字符串，算法数学定义为公式2.2。
[0045][0046]
其中lev
a,b
(i,j)表示字符串a前i个字符和字符串b前j个字符之间的编辑距离。当min(i,j)＝0，a和b之间的编辑距离为max(i,j)；当min(i,j)≠0，为如下三种情况的最小值：1.lev
a,b
(i-1,j) 1,2.lev
a,b
(i,j-1) 1,和3.他们分别表示删除ai，插入bj和替换bj，为一个指示函数，表示当ai＝bj的时候取0；当ai≠bj的时候，其值为1。
[0047]
基于编辑距离，两个字符串之间的相似度计算公式如下：
[0048][0049]
式中，|a|和|b|为字符串a和b的字符串长度，max(|a|,|b|)是字符串a和b中较大的字符串长度，lev
a,b
(i,j)为由a转换为b所需的最小编辑操作字数。
[0050]
3)遥感验证：
[0051]
从多时期高分辨率遥感影像信息中提取每个污染企业及其周边场地在各生产活动年份的高分辨率遥感影像，然后使用基于深度学习的图像语义分割模型，对每个污染企业对应的高分辨率遥感影像数据进行建筑物特征提取，以判断各生产活动年份中影像区域内是否存在建筑或者企业工厂，实现不同年份污染企业空间点位分布的遥感验证，剔除未通过遥感验证的污染企业数据，按生产年份对剩余污染企业数据进行划分，得到不同历史时期目标预测区域内的土壤cd污染企业空间点位分布数据。
[0052]
在本实施例中，基于深度学习的图像语义分割模型可采用u-net卷积神经网络模型。在对高分辨率遥感影像进行建筑物特征提取时，根据不同历史时期污染企业坐落点位的经纬度及其对应遥感影像上的信息获取该点位周边相应的遥感影像；然后对获取的遥感影像进行数据预处理，生成相应的标注数据并运用u-net卷积神经网络模型进行训练建模；再利用训练好的u-net卷积神经网络模型对高分辨率遥感影像进行图像分割，返还的结果图中的每个像元有建筑与非建筑两类；最后将属于同一生产年份的返还的结果图与污染企业数据进行比较，若污染企业所在位置在返还的结果图中存在建筑，则说明该点位的污染企业真实存在；若污染企业所在位置在返还的结果图中只存在非建筑，则说明该点位的污染企业并不存在，将该污染企业信息进行二次审查以判断该点位是否存在污染企业。此处的二次审查可通过人工审核或其他形式实现。
[0053]
u-net卷积神经网络模型属于现有技术，模型主要通过卷积与池化的步骤对图像进行特征提取，并通过反卷积的步骤将特征图恢复至原始图像的尺寸大小，卷积与池化的具体计算如下公式所示：
[0054][0055]
其中假设k为模型的卷积层，表示第k层中产生的第j个特征图。将x
ik-1
的特征图
与相应卷积核中的像元进行卷积计算并求和，b代表相应的偏置。f(x)则是该层的激励函数，此处选用relu作为激活函数。
[0056][0057]
其中down()作为下采样的函数，在平均池化的步骤下，将固定大小的像素区域里的所有像素相加，最后得到的特征图的大小变为原来的1/n。u-net卷积神经网络模型的具体结构以及训练过程不再赘述。
[0058]
s2、针对每一个历史时期，使用核密度法对该历史时期目标预测区域内的土壤cd污染企业空间点位分布数据进行核密度分析，得到污染企业密度值，同时基于该历史时期目标预测区域内各采样点的土壤cd浓度和土壤ph，通过插值方法得到目标预测区域内的土壤cd浓度空间分布和土壤ph空间分布，结合不同ph范围对应的土壤cd浓度风险筛选值，确定目标预测区域内不同位置的土壤cd风险级别。
[0059]
不同ph范围对应的土壤cd浓度风险筛选值可根据相关的标准规范或者专家经验进行确定。在本实施例中，以运用《土壤环境质量—农用地土壤污染风险管控标准》gb 15618—2018来确定不同ph范围对应的土壤cd浓度风险筛选值，进而以每个点位的土壤cd浓度风险筛选值为基准判断每个点位是否超过cd风险筛选值，从而以判断结果实现对土壤cd风险级别的划分。
[0060]
s3、针对每一个历史时期，对目标预测区域进行网格化后，统计各网格内的土壤cd风险级别与污染企业密度值，再对该历史时期土壤cd风险级别与污染企业密度值进行双变量局部空间自相关分析，得到土壤cd与污染企业的时空交互关系。
[0061]
在本实施例中，步骤s3的具体实现方法如下：针对每一个历史时期，完成土壤cd风险级别的划分后，根据目标预测区域的边界，得到其最小外接矩形所围成的范围；然后以该最小外接矩形的某个顶点开始，生成格网和对应的网格点；运用网格点提取不同历史时期污染企业核密度和土壤cd风险分级值，分别代表边界中不同历史时期污染企业的聚集程度和土壤cd风险水平；将带有土壤cd风险级别与污染企业密度值的网格点进行双变量局部莫兰分析，得到二者的时空交互关系，具体公式如下：
[0062][0063]
其中和分别是变量a和b在网格i和网格j上的不同历史时期土壤cd风险级别和污染企业密度，变量a和b分别代表网格内不同历史时期土壤cd风险级别与污染企业密度；w
ij
为网格i和网格j的空间权重矩阵，根据网格i和网格j之间的欧式距离权重所得；i
ab
表示网格i处的a属性与网格j处b属性的局部莫兰指数；当i
ab
显著为正时，则网格i处的土壤cd风险级别与网格j处的污染企业密度具有显著的局部正相关关系；当i
ab
显著为负时，则认为网格i处的土壤cd风险级别与网格j处的污染企业密度具有显著的局部负相关关系；当i
ab
不显著时，则认为网格i处的土壤cd风险级别与网格j处的污染企业密度无明显的时空交互关系。
[0064]
s4、针对每一个历史时期，根据对应的所述土壤cd风险级别、受体脆弱性和所述时空交互关系，按照预设的风险区分级标准进行目标预测区域内不同等级土壤cd风险区的识
别。
[0065]
本实施例中，将目标预测区域内不同历史时期的人口密度作为受体脆弱性指标，当人口密度高出当年我国平均人口密度10倍则认为是高脆弱性，否则为低脆弱性。土壤cd风险区的具体等级可根据实际需要进行设置，而每一个分级的标准则可以根据相关指南或者规范或者专家经验进行调整。在本实施例中，在预设的风险区分级标准中，土壤cd风险区被分为高、中、低三种风险控制区和风险不确定区。在本实施例中，运用源汇理论和前一步确定的时空交互关系来确定不同时期土壤cd风险区的等级，其采用的风险区分级标准如下：
[0066]
表1风险区分级标准
[0067][0068]
s5、获取目标预测区域内每一个历史时期对应的土壤cd风险演变影响因素，结合各历史时期的土壤cd风险区识别结果，运用斑块生成土地利用模拟模型(plus)实现对目标预测区域内未来土壤cd风险区的预测。
[0069]
本发明中，斑块生成土地利用模拟模型plus属于现有模型软件，下载链接为：
https://github.com/hpscil/patch-generating_land_use_simulation_model模型软件可以独立运行在windows vista/7/8/x64位环境，plus模型软件可以独立运行在windows vista/7/8/x64位环境。plus模型通过提取历史时期之间各等级土壤cd风险区的变化特征，并采用随机森林算法逐一对各土壤cd风险演变影响因素进行挖掘，获取各等级的土壤cd风险区的变化发生概率以及各土壤cd风险演变影响因素对变化特征的贡献，最终在变化发生概率的约束下，基于最新的土壤cd风险演变影响因素，实现对目标预测区域内未来土壤cd风险区的预测。作为一种土地利用变化模拟模型，plus模型集成了土地扩展分析策略和基于多类型随机斑块种子的元胞自动机模型，是模拟土地利用类型动态变化的有效方法。因此，基于源汇理论关系，考虑土壤cd污染水平和受体脆弱性水平，运用双变量局部莫兰指数建立土壤cd与污染企业之间的时空交互关系，同时通过斑块生成土地利用模拟模型对未来土壤cd风险进行预测，对于科学合理的土壤cd污染风险分区管理，以及指导政策制定者和利益相关者管控土壤cd污染具有迫切的现实意义。
[0070]
在本实施例中，斑块生成土地利用模拟模型实现对目标预测区域内未来土壤cd风险区的预测的具体步骤如下：首先，将两个历史时期的各等级土壤cd风险区进行叠加分析，从中提取每种等级土壤cd风险区的变化特征；其次，基于随机森林模型建立每种等级土壤cd风险区发生变化的发生概率，并计算各土壤cd风险演变影响因素对每种土地利用变化的贡献率，发生概率用于估计网格单元中每种等级土壤cd风险区的动态变化趋势；最后，通过整合时间序列的土壤cd风险区数据和模型产生的组合概率来预测未来等级土壤cd风险区，其中通过以下公式估计网格单元p未来被k等级的土壤cd风险区占用的组合概率：
[0071][0072]
其中表示网格单元p在迭代时间t从原始等级的土壤cd风险区转换到目标k等级的土壤cd风险区的组合概率；p
p,k
表示网格单元p上目标k等级的土壤cd风险区出现的概率；表示在迭代时间t对应网格单元p上目标k等级的土壤cd风险区的邻域效应；表示在迭代时间t上目标k等级的土壤cd风险区的惯性系数；scc→k表示从原始c等级的土壤cd风险区到目标k等级的土壤cd风险区的转换成本；斑块生成土地利用模拟模型使用轮盘赌选择机制来确定哪种土地利用类型将占据网格单元。
[0073]
毫无疑问，网格单元上组合概率最高的主导的土壤cd风险区等级优先被分配转换，但其他组合概率相对较低的土壤cd风险区等级仍然有机会被分配，即使机会很小。为了实现这一点，斑块生成土地利用模拟模型使用轮盘赌选择机制来确定哪种土壤cd风险区等级将占据网格单元。轮盘选择机制是指不同类型轮盘所占据的区域代表分配概率。一种土壤cd风险区等级占的区域越多，分配概率就越大，且不会剥夺其他区域被分配的机会和可能性。
[0074]
下面将上述实施例中的基于时空交互关系的土壤镉风险预测方法应用于一个具体实例中，以进一步展示本发明所能够实现的技术效果，以便本领域技术人员更好地理解本发明的实质。下述实例中的具体实现步骤框架如前述s1～s5，不再完全重复赘述，下面主要展示其具体实现细节以及技术效果。
[0075]
实施例
[0076]
选取中国东南沿海某区域作为研究区，使用本发明前述实施例s1～s5所示的方法
进行土壤镉风险预测，具体步骤如下：
[0077]
首先，基于http协议与公开网页端获得研究区与污染相关的涉土企业数据、高德poi数据、历史工业企业数据和多时期高分辨率遥感影像；其中，涉土企业数据有污染企业信息，高德poi数据有gcj-02坐标系的经纬度信息，由高德地图api下载而来，历史工业企业数据来源与中国工业企业历史数据库，高分辨率遥感影像数据来源于中国科学院资源环境科学与数据中心(https://www.resdc.cn/)。土壤cd浓度和ph数据由研究区田间采样点获取。
[0078]
将与污染相关的涉土企业中以非结构化文档数据pdf、word或图片格式存在信息，使用ocr技术将其识别为基于office open xml标准的压缩文件格式，转换为结构化数据。根据污染企业数据与历史工业企业数据和高德poi数据，分别对以上数据中的企业名称使用分词引擎jieba进行分词处理，将其分词成“城市” “企业名” “行业名” “后缀”四项；若两个数据中的企业名称经分词处理后的四项完全匹配，则将该信息返回至污染企业数据中，得到含地理位置和开始生产活动年份的污染企业数据；同时，将poi数据的gcj-02坐标系转换为wgs-84地理坐标系。
[0079]
其次，根据每个污染企业及其对应点位上的遥感影像信息，提取该点位周边相应的遥感影像。通过获取的遥感影像进行数据的预处理，并生成相应的标注数据，对u-net卷积神经网络模型进行训练建模。具体如下：
[0080]
将获得的相应点位的遥感影像与相应标注图进行切块处理，分成256*256大小的图像，并进行翻转与裁剪等数据增强的操作提高数据的多样性，以训练样本大小为10作为模型的训练输入，初始学习率为5e-5，损失函数选用二元交叉熵。在达到200k次的迭代后，模型基本达到局部最优点并收敛。在后续的模型验证过程中，对模型分割的结果以总体精度、召回率、f1-score作为模型性能的精度评价。
[0081][0082][0083][0084]
其中tp代表着是分类正确的并且类别为建筑的像元个数；tn表示分类正确的并且类别为背景的像元个数；fn表示本为建筑类别的像元被分为背景的像元个数；fp表示类别为背景的像元被分为建筑类别的像元个数。
[0085]
然后，将训练好的u-net卷积神经网络模型在研究区进行测试验证，模型对相应点位的影像进行图像分割，返还的结果图中的每个像元只有建筑与非建筑两类。若返还的结果图中存在建筑，则说明该点位的企业真实存在；若返还的结果图中只存在非建筑，则说明该点位的企业并不存在；将该结果与污染企业数据进行比较，若两者相同，则自动筛选出有效污染企业数据；若两者不同，对该企业信息进行二次审查以判断该点位是否存在企业。
[0086]
再后，使用核密度法对前述获取的不同时期污染企业点位数据进行空间插值，设定输出的分辨率为100m，得到研究区不同时期污染企业的空间分布密度图(图1c和图2c)，即完成核密度分析。运用反距离加权法预测研究区不同时期土壤cd浓度空间分布，并运用
《土壤环境质量—农用地土壤污染风险管控标准》对土壤cd风险级别进行划分(图1a和图2a)。根据研究区的研究边界生成特定100m分辨率的网格以及对应的网格点，基于网格点统计网格内不同时期土壤cd风险级别与污染企业密度值。对不同时期土壤cd风险级别与污染企业密度值进行双变量局部空间自相关分析，选取的空间相邻关系为k-nearest neighbors，得到不同时期土壤cd与污染企业的时空交互关系(图1d和图2d)。
[0087]
最后，如图2所示，根据研究区不同时期土壤cd风险级别、受体脆弱性(以人口密度为例(图1b和图2b))和时空交互关系，参考表1对研究区不同时期土壤cd进行风险识别(图3a和图3b)。随后通过收集研究区不同时期土壤cd风险演变的影响因素(包括数字高程模型、污染企业密度、pm2.5、年平均降水、年平均气温、人口密度和夜间灯光)，运用2002风险管理分区和影响因素预测2012年风险管理分区并构建模型(图3c)，以总体精度、kappa系数以及不同等级风险控制区的生产者精度和使用者精度作为评价指标，从而得到模型的精度评价结果(图4)；最后运用斑块生成土地利用模拟模型对研究区未来土壤cd风险类别进行预测(图3d)，以便后续进行监测与管理。
[0088]
以上所述的实施例只是本发明的一种较佳的方案，然其并非用以限制本发明。有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型。因此凡采取等同替换或等效变换的方式所获得的技术方案，均落在本发明的保护范围内。