本发明涉及水源污染的监测领域，尤其涉及一种基于水基因的污染溯源方法及系统。
背景技术：
溯源，追本溯源，探寻事物的根本、源头；最早是1997年欧盟为应对“疯牛病”问题而逐步建立并完善起来的食品安全管理制度。这套食品安全管理制度由政府进行推动，覆盖食品生产基地、食品加工企业、食品终端销售等整个食品产业链条的上下游，通过类似银行取款机系统的专用硬件设备进行信息共享，服务于最终消费者。一旦食品质量在消费者端出现问题，可以通过食品标签上的溯源码进行联网查询，查出该食品的生产企业、食品的产地、具体农户等全部流通信息，明确事故方相应的法律责任。此项制度对食品安全与食品行业自我约束具有相当重要的意义。现在该技术得到广泛应用，除食品以外，在药品、服饰、电子、渔船等各行各业都能见到溯源技术的影子。
流域是一类复杂的自然地理区域，它是以地表水和地下水为主要纽带，密切连接特定区域水循环、土地覆被、生态系统等自然支撑系统和人口、经济与社会等社会经济系统的综合生态地域系统。随着工业化进程的日益加快和城市化水平的不断提高，各流域的环境治理和监测任务越发重大，因此，如何保护水资源并对污染水排放进行有效的监控已成为国家可持续法阵的重要国策之一。但是，现有对不同流域的水资源监控仅仅是依靠设定在固定区域的水质自动检测站，并不能快速的监测到水源污染的来源。
技术实现要素：
有鉴于此，本发明实施例提供一种基于水基因的污染溯源方法，至少部分解决现有技术中存在的问题。
本发明的目的在于提供一种基于水基因的污染溯源方法，其特征在于，所述方法包括如下步骤：
s1：将溯源点上游3公里范围确定为溯源范围；
s2：结合各溯源点的静态特征数据，动态特征数据，筛选溯源范围内的污染源；
s3：根据监测数据判断溯源点是否被污染，若被污染，则将溯源点的污染水样与该点溯源范围内的污染源进行匹配，确定污染源行业；
s4:污染源定位。
本发明所提供的基于水基因的污染溯源方法，具有如下技术特征，所述s2包括如下步骤：
s2.1：计算不同的静态特征的污染贡献率；
s2.2：计算不同的动态特征的污染贡献率；
s2.3：根据贡献率筛选溯源范围内的污染源。
本发明所提供的基于水基因的污染溯源方法，具有如下技术特征，所述s2.1包括如下步骤：
s2.1.1：行业及排口排污量确定：
a：按行业类别划分不同污染源；
b：计算不同行业的排放量：
x1_1＝[x1_1_1,x1_1_2,…,x1_1_i,…,x1_1_h]，
其中：x1_1_i为行业i的排放量，
行业i的排放量x1_1_i＝sum(x1_1_i_m)，
其中，x1_1_i_m为行业i中企业m的排放量，
静态特征数据x1包括：二次污染普查数据、排口数据、人口数据、poi和卫星识别数据；
s2.1.2：根据环保部的监察数据计算行业和排口贡献率：
c：各行业初始贡献率x1_1＝x1_1/(sum(x1_1))，
其中：x1_1为步骤b中计算所得的排放量；
d：分别统计历史已排查行业与历史未排查行业的各行业初始贡献率，历史已排查行业的初始贡献率为x1_1_p，历史未排查行业的初始贡献率为x1_1_np；
e：分别计算历史已排查行业的行业贡献率和历史未排查行业的行业贡献率，
历史已排查行业的行业贡献率y1_1＝1/(1+e^(-5*(mu/m-0.4)))*x1_1_p，
其中：m为某行业排查发现问题的次数，mu为该行业排查得总次数；
历史未排查行业贡献率
y1_2＝(mean(y1_1)+1/(1+e^(-5*(nu/n-0.4))))*0.5*x1_1_np
其中：n为排查的总次数，nu为排查中发现问题得次数；
f：通过上述步骤得到的历史已排查行业的行业贡献率和历史未排查行业的行业贡献率计算行业及排口贡献率y1，
y1＝y1_1andy1_2。
本发明所提供的基于水基因的污染溯源方法，具有如下技术特征，所述s2.2包括如下步骤：
s2.2.1：计算各类面源污染贡献量：
y2_1＝x2_q*x2_l*x2_v，
其中：x2_q为点位流量，x2_l为降雨量数据，x2_v为土地利用数据；
s2.2.2：计算各类面源污染贡献率：
y2＝(y2_1-min(y2_1))/(max(y2_1)-min(y2_1))。
本发明所提供的基于水基因的污染溯源方法，具有如下技术特征，所述s2.3包括如下步骤：
s2.2.1:融合静态特征数据与动态特征数据计算综合贡献率：
y3＝y1andy2；
s2.2.2：根据贡献率筛选污染源：
y4＝[y3>u]，其中，u为贡献率阈值。
本发明所提供的基于水基因的污染溯源方法，具有如下技术特征，所述s3包括如下步骤：
s3.1：污染时刻的光谱谱图绘制；
s3.2：将污染源光谱谱图与污染源的光谱谱图进行特征相似度对比；
s3.3：确定污染行业。
本发明所提供的基于水基因的污染溯源方法，具有如下技术特征，所述s3.1包括如下步骤：
s3.1.1：对污染时刻的光谱数据进行去噪处理，
利用小波技术对光谱数据进行去噪处理，去除噪声后的光谱数据z1为：
z1＝f(z)
其中：z为从光学仪器中获取的200-712mn吸光度数据，光谱数据每2nm一个数据，即z的维度为256；
s3.1.2：波段筛选，
波段重要性分析
band_importance＝f(z1,label)
其中：band_importance为各波处重要性系数，label为光谱数据的行业标签，
计算累积重要性，并根据帕累托法则选取重要的波段，最终筛选204-324nm波段的数据s，
s＝z2[204:324]；
s3.1.3：光谱特征提取：
s3.1.3.1：光谱数据归一化：
s＝s/s[u]
其中：u为共性波长；
s3.1.3.2：光谱特征提取：
t1＝ds(m)/dm
t2＝d(ds(m)/dm)/dm
其中：m为波长，t1为吸光度的一阶导数，t2为吸光度的二阶导数；
s3.1.4：加强特征形成谱图矩阵，所形成的谱图矩阵h：
h＝h1+h2[::-1]，
其中h1＝f(t1,r)，
h2＝f(t2,r)，
r为小波函数。
本发明所提供的基于水基因的污染溯源方法，具有如下技术特征，所述s3.2包括如下步骤：
s3.2.1谱图矩阵对比：
p1＝pnsr(hi,h)/40，
其中：hi为污染源谱图，h为污染时刻谱图，pnsr为峰值信噪比；
s3.2.2光谱数据对比：
s3.2.2.1光谱数据去噪
s3.2.2.2归一化
s＝(s-min(s))/(max(s)-min(s))；
s3.2.3相似度对比：
p2＝(1+abs(s’))/(1+abs(s’)+2*abs(s’-si))，其中：s’＝s/s[1]，s为污染时刻光谱数据，si为污染源数据。
本发明所提供的基于水基因的污染溯源方法，具有如下技术特征，所述s3.3包括如下步骤：
s3.3.1：计算权重：
v1＝sta(p1)
v2＝sta(p2)
sta()为方差公式
w1＝v1/(v1+v2)
w2＝v2/(v1+v2)；
s3.3.2：污染行业确定：
p＝max(w1*s1+w2*s2)。
本发明所提供的基于水基因的污染溯源方法，具有如下技术特征，所述步骤s4包括如下步骤：
s4.1：对水质超标位置的上游区域进行网络划分；
s4.2：对上述划分的网络进行评分；
s4.3：选取上述评分最高的网络进行排查。
本发明所提供的基于水基因的污染溯源方法，具有如下技术特征，所述步骤s4.2包括如下步骤：
s4.2.1：对网络的距离进行评分：
超标点位与网格中心距离
d1＝((dlat-xlat)^2+(dlng-xlng)^2)^0.5
其中dlat，dlng为网格中心坐标，xlat，xlng为超标点位坐标网格中心与河岸距离
d2＝min(((hlat-xlat)^2+(hlng-xlat)^2)^0.5)
其中，hlat，hlng为网格中心做河岸线垂线与河岸线相交的坐标则网络的距离评分d为：d＝d1+d2
d＝(d-min(d))/(max(d)-min(d))；
s4.2.2：对网络的密度进行评分：
p＝n/s
其中n为网格内企业数量，s为网格面积，
p＝(p-min(p))/(max(p)-min(p))；
s4.2.3：对网络的排查情况进行评分：
历史已排查网格评分
r1＝2*(m/mu)/((0.5+0.5*(m/(n/nu)))
其中，m为某点位排查次数，mu为该点位排查得总次数，n为排查的总次数，nu为排查中发现问题得次数
历史未排查网格评分
r2＝(mean(r1)+1/(1+e^(-5*(nu/n-0.4))))*0.5
其中n为某排查总次数，nu为排查中发现问题得次数
r＝r1andr2；
s4.2.4：根据距离评分、密度评分和排查情况评分，相加得到该网络评分。
有益效果
本发明所提供的基于水基因的污染溯源方法，将溯源这一制度引入到城市河道、流域等污染的检测领域中，通过污染行业和污染来源的数据进行前期统计，后续分析对比污染的特征，能更清晰明确的反应出水污染的污染来源和污染行业。
具体实施方式
以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需说明的是，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
要说明的是，下文描述在所附权利要求书的范围内的实施例的各种方面。应显而易见，本文中所描述的方面可体现于广泛多种形式中，且本文中所描述的任何特定结构及/或功能仅为说明性的。基于本发明，所属领域的技术人员应了解，本文中所描述的一个方面可与任何其它方面独立地实施，且可以各种方式组合这些方面中的两者或两者以上。举例来说，可使用本文中所阐述的任何数目个方面来实施设备及/或实践方法。另外，可使用除了本文中所阐述的方面中的一或多者之外的其它结构及/或功能性实施此设备及/或实践此方法。
另外，在以下描述中，提供具体细节是为了便于透彻理解实例。然而，所属领域的技术人员将理解，可在没有这些特定细节的情况下实践所述方面。
本发明实施例提供一种基于水基因的污染溯源方法，该方法包括如下步骤：
s1：将溯源点上游3公里范围确定为溯源范围；
s2：结合各溯源点的静态特征数据，动态特征数据，筛选溯源范围内的污染源；
s3：根据监测数据判断溯源点是否被污染，若被污染，则将溯源点的污染水样与该点溯源范围内的污染源进行匹配，确定污染源行业；
s4:污染源定位。
在本发明的部分实施例中，所述s2包括如下步骤：
s2.1：计算不同的静态特征的污染贡献率；
s2.2：计算不同的动态特征的污染贡献率；
s2.3：根据贡献率筛选溯源范围污染源。
在本发明的部分实施例中，输入静态数据，再对静态数据进行分类统计，即将静态数据中的企业按照行业划分，统计该区域内各行业的排放量及排口的排放量。再输入各行业排水量及环保部门的监察数据，计算各类污染来源的贡献率，即企业问题大的贡献率大。具体步骤如下：
s2.1.1：行业及排口排污量确定：
a：按行业类别划分不同污染源：行业类别的划分依据为《国民经济行业分类》，根据关键字初步分类后进行人工核查是否有误。
b：计算不同行业的排放量：
x1_1＝[x1_1_1,x1_1_2,…,x1_1_i,…,x1_1_h]，
其中：x1_1_i为行业i的排放量，
行业i的排放量x1_1_i＝sum(x1_1_i_m)，
其中，x1_1_i_m为行业i中企业m的排放量，
静态特征数据x1包括：二次污染普查数据、排口数据、人口数据、poi和卫星识别数据；二次污染普查数据指工业、农业、生活污染的信息，静态数据的获得通过网络爬取，客户提供以及公司摸排等方式。
s2.1.2：根据环保部的监察数据计算行业和排口贡献率：
c：各行业初始贡献率x1_1＝x1_1/(sum(x1_1))，
其中：x1_1为步骤b中计算所得的排放量；
d：分别统计历史已排查行业与历史未排查行业的各行业初始贡献率，历史已排查行业的初始贡献率为x1_1_p，历史未排查行业的初始贡献率为x1_1_np；
e：分别计算历史已排查行业的行业贡献率和历史未排查行业的行业贡献率，
历史已排查行业的行业贡献率y1_1＝1/(1+e^(-5*(mu/m-0.4)))*x1_1_p，其中：m为某行业排查发现问题的次数，mu为该行业排查得总次数；
历史未排查行业贡献率
y1_2＝(mean(y1_1)+1/(1+e^(-5*(nu/n-0.4))))*0.5*x1_1_np
其中：n为排查的总次数，nu为排查中发现问题得次数；
f：通过上述步骤得到的历史已排查行业的行业贡献率和历史未排查行业的行业贡献率计算行业及排口贡献率y1，
y1＝y1_1andy1_2。
在本发明的部分实施例中，输入流量数据、雨情数据以及土地利用数据，分析面源冲刷对水质的影响。面源污染又称非点源污染，主要由土壤泥沙颗粒、氮磷等营养物质、农药、各种大气颗粒物等组成，通过地表径流、土壤侵蚀、农田排水等方式进入水、土壤或大气环境。s2.2的具体步骤如下：
s2.2.1：计算各类面源污染贡献量：
y2_1＝x2_q*x2_l*x2_v，
其中：x2_q为点位流量，x2_l为降雨量数据，x2_v为土地利用数据；
s2.2.2：计算各类面源污染贡献率：
y2＝(y2_1-min(y2_1))/(max(y2_1)-min(y2_1))。
在上述实施例中，考虑了面源污染对水源环境污染的影响，使得污染的监测、治理和管理更加精确。
在本发明的部分实施例中，输入静态数据与动态数据贡献率值，确定该区域所有的污染源，并根据贡献率大小，剔除个别贡献率小的污染源。s2.3具体步骤如下：
s2.2.1:融合静态特征数据与动态特征数据计算综合贡献率：
y3＝y1andy2；
s2.2.2：根据贡献率筛选污染源：
y4＝[y3>u]，其中，u为贡献率阈值。
在本发明的部分实施例中，所述s3包括如下步骤：
s3.1：污染时刻的光谱谱图绘制；
s3.2：将污染源光谱谱图与污染源的光谱谱图进行特征相似度对比；
s3.3：确定污染行业。
在本发明的部分实施例中，所述s3.1包括如下步骤：
s3.1.1：对污染时刻的光谱数据进行去噪处理，
利用小波技术对光谱数据进行去噪处理，去除噪声后的光谱数据z1为：
z1＝f(z)
其中：z为从光学仪器中获取的200-712mn吸光度数据，光谱数据每2nm一个数据，即z的维度为256；
s3.1.2：波段筛选，
波段重要性分析
band_importance＝f(z1,label)
其中：band_importance为各波处重要性系数，label为光谱数据的行业标签，
计算累积重要性，即选取波段重要性系数的总和，并根据帕累托法则选取重要的波段，最终筛选204-324nm波段的数据s，
s＝z2[204:324]；
s3.1.3：光谱特征提取：
s3.1.3.1：光谱数据归一化：
s＝s/s[u]
其中：u为共性波长；
s3.1.3.2：光谱特征提取：
t1＝ds(m)/dm
t2＝d(ds(m)/dm)/dm
其中：m为波长，t1为吸光度的一阶导数，t2为吸光度的二阶导数；
s3.1.4：加强特征形成谱图矩阵，所形成的谱图矩阵h：
h＝h1+h2[::-1]，
其中h1＝f(t1,r)，
h2＝f(t2,r)，
r为小波函数。
在本发明的部分实施例中，所述s3.2包括如下步骤：
s3.2.1谱图矩阵对比：
p1＝pnsr(hi,h)/40，
其中：hi为污染源谱图，h为污染时刻谱图，pnsr为峰值信噪比；
s3.2.2光谱数据对比：
s3.2.2.1光谱数据去噪
s3.2.2.2归一化
s＝(s-min(s))/(max(s)-min(s))；
s3.2.3相似度对比：
p2＝(1+abs(s’))/(1+abs(s’)+2*abs(s’-si))，其中：s’＝s/s[1]，s为污染时刻光谱数据，si为污染源数据。
在本发明的部分实施例中，所述s3.3包括如下步骤：
s3.3.1：计算权重：
v1＝sta(p1)
v2＝sta(p2)
sta()为方差公式
w1＝v1/(v1+v2)
w2＝v2/(v1+v2)；
s3.3.2：污染行业确定：
p＝max(w1*s1+w2*s2)。
在本发明的部分实施例中，所述步骤s4包括如下步骤：
s4.1：对水质超标位置的上游区域进行网络划分，划分的网络面积为300*300米；
s4.2：对上述划分的网络进行评分；
s4.3：选取上述评分最高的网络进行排查。
在本发明的部分实施例中，所述步骤s4.2包括如下步骤：
s4.2.1：对网络的距离进行评分：
超标点位与网格中心距离
d1＝((dlat-xlat)^2+(dlng-xlng)^2)^0.5
其中dlat，dlng为网格中心坐标，xlat，xlng为超标点位坐标，该坐标为经纬度坐标。
网格中心与河岸距离
d2＝min(((hlat-xlat)^2+(hlng-xlat)^2)^0.5)
其中，hlat，hlng为网格中心做河岸线垂线与河岸线相交的坐标，该坐标为经纬度坐标。
则网络的距离评分d为：d＝d1+d2
d＝(d-min(d))/(max(d)-min(d))；
s4.2.2：对网络的密度进行评分：
p＝n/s
其中n为网格内企业数量，s为网格面积，
p＝(p-min(p))/(max(p)-min(p))；
s4.2.3：对网络的排查情况进行评分：
历史已排查网格评分
r1＝2*(m/mu)/((0.5+0.5*(m/(n/nu)))
其中，m为某点位排查次数，mu为该点位排查得总次数，n为排查的总次数，nu为排查中发现问题得次数
历史未排查网格评分
r2＝(mean(r1)+1/(1+e^(-5*(nu/n-0.4))))*0.5
其中n为某排查总次数，nu为排查中发现问题得次数
r＝r1andr2；
s4.2.4：根据距离评分、密度评分和排查情况评分，相加得到该网络评分。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。