CN114943194A - 一种基于地质统计学的河流污染溯源方法 - Google Patents

一种基于地质统计学的河流污染溯源方法 Download PDF

Info

Publication number
CN114943194A
CN114943194A CN202210527841.2A CN202210527841A CN114943194A CN 114943194 A CN114943194 A CN 114943194A CN 202210527841 A CN202210527841 A CN 202210527841A CN 114943194 A CN114943194 A CN 114943194A
Authority
CN
China
Prior art keywords
formula
river
pollution source
water quality
theta
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202210527841.2A
Other languages
English (en)
Other versions
CN114943194B (zh
Inventor
崔婷婷
高岐涛
龙玉桥
盖永伟
宇润练
李伟
王小红
陈建东
赵士文
刘勇
刘克琳
陈韬
杨杰
卜建东
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Taiyuan Institute Of Water Resources And Soil Conservation
Nanjing Hydraulic Research Institute of National Energy Administration Ministry of Transport Ministry of Water Resources
Original Assignee
Taiyuan Institute Of Water Resources And Soil Conservation
Nanjing Hydraulic Research Institute of National Energy Administration Ministry of Transport Ministry of Water Resources
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Taiyuan Institute Of Water Resources And Soil Conservation, Nanjing Hydraulic Research Institute of National Energy Administration Ministry of Transport Ministry of Water Resources filed Critical Taiyuan Institute Of Water Resources And Soil Conservation
Priority to CN202210527841.2A priority Critical patent/CN114943194B/zh
Publication of CN114943194A publication Critical patent/CN114943194A/zh
Application granted granted Critical
Publication of CN114943194B publication Critical patent/CN114943194B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/12—Computing arrangements based on biological models using genetic models
    • G06N3/126—Evolutionary algorithms, e.g. genetic algorithms or genetic programming
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q30/00—Commerce
    • G06Q30/018—Certifying business or products
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/10—Services
    • G06Q50/26—Government or public services
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A20/00—Water conservation; Efficient water supply; Efficient water use
    • Y02A20/152—Water filtration

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Business, Economics & Management (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Optimization (AREA)
  • Data Mining & Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Health & Medical Sciences (AREA)
  • Evolutionary Computation (AREA)
  • Software Systems (AREA)
  • Computational Mathematics (AREA)
  • Biophysics (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Marketing (AREA)
  • Economics (AREA)
  • Databases & Information Systems (AREA)
  • General Business, Economics & Management (AREA)
  • Artificial Intelligence (AREA)
  • Strategic Management (AREA)
  • Computer Hardware Design (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Geometry (AREA)
  • Development Economics (AREA)
  • Tourism & Hospitality (AREA)
  • Medical Informatics (AREA)
  • Genetics & Genomics (AREA)
  • Human Resources & Organizations (AREA)
  • Computational Linguistics (AREA)
  • Educational Administration (AREA)
  • Biomedical Technology (AREA)
  • Molecular Biology (AREA)

Abstract

本发明公开了一种基于地质统计学的河流污染溯源方法，应用于环境污染技术领域，具体步骤如下：获得区面上的多时段水质与流速数据；利用地质统计学方法构建污染源溯源的反演数学模型；根据多时段水质与流速数据初步判断污染源所处的区域，并假定可能的污染源位置；利用流速数据建立分析河段的水动力及水质模拟模型，并通过水动力及水质模拟模型计算传递函数；采用GPU并行加速方法改进原最小化问题的计算方法求解污染源溯源的反演数学模型，通过比较模型计算结果与实测的水质与流速数据拟合程度，确定污染源的位置与排放过程。该方法为了提高计算效率采用GPU并行加速方法改进原最小化问题的计算方法求解污染源溯源的反演数学模型。

Description

一种基于地质统计学的河流污染溯源方法
技术领域
本发明涉及环境保护技术领域，更具体的说是涉及一种基于地质统计学的河流污染溯源方法。
背景技术
由水污染导致的水质型缺水已经成为制约国民经济可持续发展的重要因素。水污染的来源主要有城市生活污水、农业面源污染、工业污水等。目前，我国不断加大对环境保护的力度，在水污染治理方面已经开展了大量工作，水环境质量在整体上有较大的提升。然而，仍存在部分偷排或污染物意外泄漏的情况。这些污染排放存在隐蔽、随机的特点。一旦发生，会对水环境产生较大的影响，但是由于其排放不规律等原因，难以快速确定污染源，为水环境的管理与治理带来了一定困难。
传统的污染源调查工作方法主要依据水质监测数据来判定污染源的位置与强度。水质监测数据主要来自于固定的监测站点与人工取样，这些方法虽然精确度高，但仅能提供有限点位的水质数据，耗时耗力，调查效率低，无法反应区面上的水质变化，也难以查找隐蔽、随机的污染源。基于无人机平台的高光谱或多光谱水质分析技术，侦测范围广、使用条件灵活、数据分辨率高，近年来已出现了大量的相关研究，为高效地获取较大区面上水质数据提供了技术基础。根据多时段、大区面上的水质数据，利用地质统计学方法反演污染源位置与排放过程，能够极大的节省甄别隐蔽污染源所需的人力物力，有利于提高污染源排查和相关环境保护工作的效率和准确性。
进一步，在溯源问题较为复杂时，由于计算负荷较大，需要的计算时间较长，对于计算溯源问题时，计算量比较大的问题是本领域技术人员亟需解决的问题。
发明内容
有鉴于此，本发明提供了一种针对偷排或污染物意外泄漏时污染源具有隐蔽、随机、难以确定的特点，提出一种新的污染源溯源方法。该方法先利用无人机平台的高光谱或多光谱水质分析技术为该方法获得区面上的多时段水质数据，再利用地质统计学方法反演出隐蔽污染源的位置与排放过程，同时为了提高计算效率采用GPU并行加速方法改进原最小化问题的计算方法求解污染源溯源的反演数学模型。
为了实现上述目的，本发明提供如下技术方案：
一种基于地质统计学的河流污染溯源方法，具体步骤如下：
获得区面上的多时段水质与流速数据；
利用地质统计学方法构建污染源溯源的反演数学模型；
根据多时段水质与流速数据初步判断污染源所处的区域，并假定可能的污染源位置；
利用流速数据建立分析河段的水动力及水质模拟模型，并通过水动力及水质模拟模型计算传递函数；
采用GPU并行加速方法改进原最小化问题的计算方法求解污染源溯源的反演数学模型，通过比较模型计算结果与实测的水质与流速数据拟合程度，确定污染源的位置与排放过程。
可选的，在上述的一种基于地质统计学的河流污染溯源方法中，构建污染源溯源的反演数学模型具体步骤如下：
假定河流中污染物未发生化学反应或发生一级化学反应，将污染物排放过程与河流中污染物浓度观测值间的关系概化为式1：
z＝Hs+v 式1；
式中：z为一个由观测数据组成的m×1维向量；s为待求未知函数经离散后形成的一个n×1维向量，代表污染物排放过程；v为测量误差组成的m×1维向量；H为由传递函数构成的灵敏度矩阵，代表了污染物在河水中的迁移过程；
假定v的均值为0，其协方差矩阵为R，s为一个随机向量，其数学期望和协方差分别由式2和式3表示：
E[s]＝Xβ 式2；
Q(θ)＝E[(s-Xβ)(s-Xβ)T] 式3；
式中：X为一个已知的n×p的矩阵；β为p个漂移系数组成的向量；Q(θ)为关于未知参数θ的已知函数，通常取高斯分布函数。
可选的，在上述的一种基于地质统计学的河流污染溯源方法中，根据污染物浓度的非负性，对未知函数s进行如下变换：
Figure BDA0003645299820000036
将式4代入式1得到
Figure BDA0003645299820000031
变换后，
Figure BDA0003645299820000032
是非线性的，需要通过迭代方法求得θ和
Figure BDA0003645299820000033
最终利用式6求出未知函数s的最佳估计
Figure BDA0003645299820000034
Figure BDA0003645299820000035
可选的，在上述的一种基于地质统计学的河流污染溯源方法中，通过水动力及水质模拟模型计算传递函数具体步骤如下：
由于假定河流中污染物未发生化学反应或发生一级化学反应，污染物迁移过程用线性偏微分方程表示：
Figure BDA0003645299820000041
式中c为污染物浓度，u为流速，D为弥散系数，t为时间，x为迁移距离，r(c)为反应项；
此时，污染物浓度由传递函数方程表示：
Figure BDA0003645299820000042
其中s(t)为实际排放过程，f(·)为传递函数。
可选的，在上述的一种基于地质统计学的河流污染溯源方法中，通过差分方法计算传递函数，从而得到含有一级化学反应过程的传递函数；
通过变量替换将式8变为式9：
Figure BDA0003645299820000043
当式9中t>0时，s(T-t)＝1，t≤0时，s(T-t)＝0，得到穿透曲线方程：
Figure BDA0003645299820000044
使式10对时间t求导，即可得到传递函数在坐标x时刻t的值；
Figure BDA0003645299820000045
根据假定的污染源位置，利用正向的河流污染物迁移模型计算s(t)＝1条件下各观测点的浓度曲线，然后采用向后差分方法，式12计算各观测点在各个时刻的传递函数值，从而得到由传递函数所组成的灵敏度矩阵H，最后利用地质统计学方法反演出相应污染源污染物的排放过程；
Figure BDA0003645299820000046
可选的，在上述的一种基于地质统计学的河流污染溯源方法中，待求未知函数离散形式s的求解过程可分成两个步骤，首先是确定未知参数θ，而后估计未知函数s；
当式13表示的概率取最大值时，θ的取值即为待求的参数值；
Figure BDA0003645299820000051
Σ＝HQHT+R 式14；
Ξ＝Σ-1-Σ-1HX(XTHTΣ-1HX)-1XTHTΣ-1 式15；
式13描述的极值问题，等同于求式16的最小值；
Figure BDA0003645299820000052
L(θ)对θ的导数取0值时，L(θ)取最小值；当迭代计算收敛后，再利用θ求Q，并代入式17，解得到n×m的系数矩阵Λ和p×n的拉格朗日乘子阵M；
Figure BDA0003645299820000053
对未知函数s的最佳估计
Figure BDA0003645299820000054
及其协方差分别如式18和式19；
Figure BDA0003645299820000055
V＝-XM+Q-QHTΛT 式19；
得到假定污染源对应的排放过程。
可选的，在上述的一种基于地质统计学的河流污染溯源方法中，采用GPU并行加速方法改进原最小化问题的计算方法具体步骤如下：
要求L(θj)<L(θj-1)，j代表第j次迭代计算的次数；使用系数λ控制目标函数的搜索方向，通过反复调整λ直到搜索结果满足上述要求，要求第j次搜索使用的系数θj对应的搜索结果应满足L(θj)<L(θj-1)；当计算结果满足该要求后，根据θ的收敛标准，确定第j次的搜索结果θj是否为最终结果，或是否进入第j+1次运算，会通过大量的循环迭代计算找到合适的λ和θ。
可选的，在上述的一种基于地质统计学的河流污染溯源方法中，将GPU中Grid中的若干个Block作为一个Block组，每个Block组执行一个λj所对应的计算工作，组内多个Block执行最速下降法中的矩阵运算；各个Block组的计算结果中，满足L(θj)＝min(L(θ1)，L(θ2)，L(θ3)，…，L(θp))的λj为所求，并根据θ的收敛标准，决定第j个的搜索结果θj是否为最终结果，或是否进入下次运算。
可选的，在上述的一种基于地质统计学的河流污染溯源方法中，将假定污染源与排放过程输入正向模型中，计算各离散网格的水质与流速，根据计算结果与测量结果的拟合效果确定污染源的位置与排放过程；采用下式描述拟合效果：
Figure BDA0003645299820000061
式中
Figure BDA0003645299820000062
为离散网格i的测量值，ci为离散网格i的计算值。
经由上述的技术方案可知，与现有技术相比，本发明公开提供了一种基于地质统计学的河流污染溯源方法，充分利用了无人机平台的灵活性，以及光谱分析技术便于取得区面数据的优势，克服了区域水质取样密度低及取样频次低的不足，在水质与流速数据的基础上，利用地质统计学方法反演污染物在河水中的迁移过程，通过对比分析计算结果与测量数据的拟合程度，确定隐蔽污染源的位置与排放过程。与传统水质分析方法或人工探查方法相比，能够反映区面上的水质信息，便于查找隐蔽污染源，同时提高污染源调查工作的效率，有利于提高河流污染治理及管理的技术水平。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。
图1为本发明的方法流程图；
图2为本发明的并行加速示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明实施例公开了一种基于地质统计学的河流污染溯源方法，该方法先利用无人机平台的高光谱或多光谱水质分析技术为该方法获得区面上的多时段水质数据，再利用地质统计学方法反演出隐蔽污染源的位置与排放过程。具体步骤如下：
(1)利用无人机平台的高光谱或多光谱水质分析技术为该方法获得区面上的多时段水质与流速数据。
(2)利用地质统计学方法构建污染源溯源的反演数学模型。
(3)根据已知数据与相关调查初步判断污染源所处的区域，并假定可能的污染源位置。
(4)利用流速数据建立分析河段的水动力及水质模拟模型，并得通过模型计算传递函数。
(5)求解污染源溯源的反演数学模型，通过比较模型计算结果与实测的水质与流速数据拟合程度，分析确定污染源的位置与排放过程。
本发明实施例的的主要内容如下：
1.获取区面上的水质与流速数据
无人机平台和高光谱或多光谱水质分析技术的作用是为本发明提供区面上的水质参数与流速数据。目前，已有大量利用无人机平台和高光谱或多光谱手段分析水质与流速的研究，本发明此处仅简述水质与流速数据获取的过程。
(1)水质参数
首先，使用搭载有高光谱或多光谱成像设备的无人机，在拟分析河段上按一定航迹飞行，并获取飞行路线上的高光谱或多光谱低空遥感影像。其次，利用实测水质数据和遥感影像，确定水质参数与光谱曲线间的定量关系。再次，利用机器学习技术构建水质参数反演模型。最终，利用水质参数反演模型分析不同时段的影像数据，从而获得分析河段上相应时段的水质参数空间分布数据。
(2)流速参数
利用获得的影像数据，分析河水表面水流特征和漂浮物交错轨迹的流动时间计算流动速度，然后分析其分布并间接地计算水流速度，最终得到不同时段的流速空间分布数据。
(3)水质与流速参数的空间分配
根据拟使用的河流水质模拟模型，采用结构性或非结构性网格将分析河段离散化。根据离散后的空间信息，将水质与流速参数分配到各网格中。
2.建立溯源数学模型
(1)地质统计学溯源模型
假定河流中污染物未发生化学反应或发生一级化学反应，将污染物排放过程与河流中污染物浓度观测值间的关系概化为式1：
z＝Hs+v 式1
式中：z为一个由观测数据组成的m×1维向量；s为待求未知函数经离散后形成的一个n×1维向量，代表污染物排放过程；v为测量误差组成的m×1维向量；H为由传递函数构成的灵敏度矩阵，代表了污染物在河水中的迁移过程。
假定v的均值为0，其协方差矩阵为R，s为一个随机向量，其数学期望和协方差分别由式2和式3表示。
E[s]＝Xβ 式2
Q(θ)＝E[(s-Xβ)(s-Xβ)T] 式3
式中：X为一个已知的n×p的矩阵；β为p个漂移系数组成的向量；Q(θ)为关于未知参数θ的已知函数，通常取高斯分布函数。
(2)污染物浓度非负性约束
考虑污染物浓度的非负性，对未知函数s进行如下变换：
Figure BDA0003645299820000091
将式4代入式1得到
Figure BDA0003645299820000092
实施上述变换后，
Figure BDA0003645299820000093
是非线性的，需要通过迭代方法求得θ和
Figure BDA0003645299820000094
最终利用式6求出未知函数s的最佳估计
Figure BDA0003645299820000095
Figure BDA0003645299820000096
3.计算传递函数
为便于说明问题，此处用一维问题说明传递函数的计算方法。由于假定河流中污染物未发生化学反应或发生一级化学反应，污染物迁移过程可用线性偏微分方程表示。
Figure BDA0003645299820000097
式中c为污染物浓度，u为流速，D为弥散系数，t为时间，x为迁移距离，r(c)为反应项。
此时，污染物浓度可由传递函数方程表示：
Figure BDA0003645299820000098
其中s(t)为实际排放过程，f(·)为传递函数。在迁移问题较复杂的情况下，难以得到传递函数的解析形式。因此，本发明通过差分方法计算传递函数，从而得到含有一级化学反应过程的传递函数。
首先，通过变量替换将式8变为式9：
Figure BDA0003645299820000101
当式9中t>0时，s(T-t)＝1，t≤0时，s(T-t)＝0，得到穿透曲线方程：
Figure BDA0003645299820000102
使式10对时间t求导，即可得到传递函数在坐标x时刻t的值。
Figure BDA0003645299820000103
因此，可先根据假定的污染源位置，利用正向的河流污染物迁移模型计算s(t)＝1条件下各观测点的浓度曲线，然后采用向后差分方法(式12)计算各观测点在各个时刻的传递函数值，从而得到由传递函数所组成的灵敏度矩阵H，最后利用地质统计学方法反演出相应污染源污染物的排放过程。
Figure BDA0003645299820000104
4.求解溯源数学模型
待求未知函数离散形式s的求解过程可分成两个步骤，首先是确定未知参数θ，而后估计未知函数s。当式13表示的概率取最大值时，θ的取值即为待求的参数值。
Figure BDA0003645299820000105
Σ＝HQHT+R 式14
Ξ＝Σ-1-Σ-1HX(XTHTΣ-1HX)-1XTHTΣ-1 式15
式13描述的极值问题，等同于求式16的最小值。
Figure BDA0003645299820000111
L(θ)对θ的导数取0值时，L(θ)取最小值。可采用Gauss-Newton迭代法、最速下降法等方法求解此极值问题，也可以结合遗传算法等智能优化方法求解此问题。当迭代计算收敛后，再利用θ求Q，并代入式17，解得到n×m的系数矩阵Λ和p×n的拉格朗日乘子阵M。
Figure BDA0003645299820000112
对未知函数s的最佳估计
Figure BDA0003645299820000113
及其协方差分别如式18和式19所示。
Figure BDA0003645299820000114
V＝-XM+Q-QHTΛT 式19
此时，得到假定污染源对应的排放过程。
5.确定污染源的位置与排放过程
将假定污染源与排放过程输入正向模型中，计算各离散网格的水质与流速，根据计算结果与测量结果的拟合效果确定污染源的位置与排放过程。本发明采用下式描述拟合效果：
Figure BDA0003645299820000115
式中
Figure BDA0003645299820000116
为离散网格i的测量值，ci为离散网格i的计算值。
6.并行加速
采用地质统计学溯源时，求解式16的最小化问题需要进行大量迭代计算。在溯源问题较为复杂时，由于计算负荷较大，需要的计算时间较长。为提高计算效率，采用GPU并行加速方法改进原最小化问题的计算方法。
以最速下降法为例，该方法法要求每次迭代后的目标函数值应是下降的，即要求L(θj)<L(θj-1)，j代表第j次迭代计算的次数。该方法使用系数λ控制了目标函数的搜索方向，即通过反复调整λ直到搜索结果满足上述要求，即要求第j次搜索使用的系数θj对应的搜索结果应满足L(θj)<L(θj-1)。当计算结果满足该要求后，根据θ的收敛标准，确定第j次的搜索结果θj是否为最终结果，或是否进入第j+1次运算。因此，在串行的计算方法中，会通过大量的循环迭代计算找到合适的λ和θ。
针对这一特点，本文将GPU中Grid(网格)中的几个Block(线程块)作为一个Block组，每个Block组执行一个λj所对应的计算工作，组内多个Block执行最速下降法中的矩阵运算。各个Block组的计算结果中，满足L(θj)＝min(L(θ1)，L(θ2)，L(θ3)，…，L(θn))的λj为所求，并根据θ的收敛标准，决定第j个的搜索结果θj是否为最终结果，或是否进入下次运算。
由于遗传算法等智能优化算法也可以采用GPU并行的方式来提高计算效率。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的装置而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims (9)

1.一种基于地质统计学的河流污染溯源方法，其特征在于，具体步骤如下：
获得区面上的多时段水质与流速数据；
利用地质统计学方法构建污染源溯源的反演数学模型；
根据多时段水质与流速数据初步判断污染源所处的区域，并假定可能的污染源位置；
利用流速数据建立分析河段的水动力及水质模拟模型，并通过水动力及水质模拟模型计算传递函数；
采用GPU并行加速方法改进原最小化问题的计算方法求解污染源溯源的反演数学模型，通过比较模型计算结果与实测的水质与流速数据拟合程度，确定污染源的位置与排放过程。
2.根据权利要求1所述的一种基于地质统计学的河流污染溯源方法，其特征在于，构建污染源溯源的反演数学模型具体步骤如下：
假定河流中污染物未发生化学反应或发生一级化学反应，将污染物排放过程与河流中污染物浓度观测值间的关系概化为式1：
z＝Hs+v 式1；
式中：z为一个由观测数据组成的m×1维向量；s为待求未知函数经离散后形成的一个n×1维向量，代表污染物排放过程；v为测量误差组成的m×1维向量；H为由传递函数构成的灵敏度矩阵，代表了污染物在河水中的迁移过程；
假定v的均值为0，其协方差矩阵为R，s为一个随机向量，其数学期望和协方差分别由式2和式3表示：
E[s]＝Xβ 式2；
Q(θ)＝E[(s-Xβ)(s-Xβ）T] 式3；
式中：X为一个已知的n×p的矩阵；β为p个漂移系数组成的向量；Q(θ)为关于未知参数θ的已知函数，通常取高斯分布函数。
3.根据权利要求2所述的一种基于地质统计学的河流污染溯源方法，其特征在于，根据污染物浓度的非负性，对未知函数s进行如下变换：
Figure FDA0003645299810000021
将式4代入式1得到
Figure FDA0003645299810000022
变换后，
Figure FDA0003645299810000023
是非线性的，需要通过迭代方法求得θ和
Figure FDA0003645299810000024
最终利用式6求出未知函数s的最佳估计
Figure FDA0003645299810000025
Figure FDA0003645299810000026
4.根据权利要求1所述的一种基于地质统计学的河流污染溯源方法，其特征在于，通过水动力及水质模拟模型计算传递函数具体步骤如下：
由于假定河流中污染物未发生化学反应或发生一级化学反应，污染物迁移过程用线性偏微分方程表示：
Figure FDA0003645299810000027
式中c为污染物浓度，u为流速，D为弥散系数，t为时间，x为迁移距离，r(c)为反应项；
此时，污染物浓度由传递函数方程表示：
Figure FDA0003645299810000028
其中s(t)为实际排放过程，f(·)为传递函数。
5.根据权利要求4所述的一种基于地质统计学的河流污染溯源方法，其特征在于，通过差分方法计算传递函数，从而得到含有一级化学反应过程的传递函数；
通过变量替换将式8变为式9：
Figure FDA0003645299810000031
当式9中t>0时，s(T-t)＝1，t≤0时，s(T-t)＝0，得到穿透曲线方程：
Figure FDA0003645299810000032
使式10对时间t求导，即可得到传递函数在坐标x时刻t的值；
Figure FDA0003645299810000033
根据假定的污染源位置，利用正向的河流污染物迁移模型计算s(t)＝1条件下各观测点的浓度曲线，然后采用向后差分方法，式12计算各观测点在各个时刻的传递函数值，从而得到由传递函数所组成的灵敏度矩阵H，最后利用地质统计学方法反演出相应污染源污染物的排放过程；
Figure FDA0003645299810000034
6.根据权利要求1所述的一种基于地质统计学的河流污染溯源方法，其特征在于，待求未知函数离散形式s的求解过程可分成两个步骤，首先是确定未知参数θ，而后估计未知函数s；
当式13表示的概率取最大值时，θ的取值即为待求的参数值；
Figure FDA0003645299810000035
Σ＝HQHT+R 式14；
Ξ＝Σ-1-Σ-1HX(XTHTΣ-1HX)-1XTHTΣ-1 式15；
式13描述的极值问题，等同于求式16的最小值；
Figure FDA0003645299810000036
L(θ)对θ的导数取0值时，L(θ)取最小值；当迭代计算收敛后，再利用θ求Q，并代入式17，解得到n×m的系数矩阵Λ和p×n的拉格朗日乘子阵M；
Figure FDA0003645299810000037
对未知函数s的最佳估计
Figure FDA0003645299810000038
及其协方差分别如式18和式19；
Figure FDA0003645299810000041
V＝-XM+Q-QHTΛT 式19；
得到假定污染源对应的排放过程。
7.根据权利要求6所述的一种基于地质统计学的河流污染溯源方法，其特征在于，采用GPU并行加速方法改进原最小化问题的计算方法具体步骤如下：
要求L(θj)<L(θj-1)，j代表第j次迭代计算的次数；使用系数λ控制目标函数的搜索方向，通过反复调整λ直到搜索结果满足上述要求，要求第j次搜索使用的系数θj对应的搜索结果应满足L(θj)<L(θj-1)；当计算结果满足该要求后，根据θ的收敛标准，确定第j次的搜索结果θj是否为最终结果，或是否进入第j+1次运算，会通过大量的循环迭代计算找到合适的λ和θ。
8.根据权利要求7所述的一种基于地质统计学的河流污染溯源方法，其特征在于，将GPU中Grid中的若干个Block作为一个Block组，每个Block组执行一个λj所对应的计算工作，组内多个Block执行最速下降法中的矩阵运算；各个Block组的计算结果中，满足L(θj)＝min(L(θ1)，L(θ2)，L(θ3)，…，L(θn))的λj为所求，并根据θ的收敛标准，决定第j个的搜索结果θj是否为最终结果，或是否进入下次运算。
9.根据权利要求1所述的一种基于地质统计学的河流污染溯源方法，其特征在于，将假定污染源与排放过程输入正向模型中，计算各离散网格的水质与流速，根据计算结果与测量结果的拟合效果确定污染源的位置与排放过程；采用下式描述拟合效果：
Figure FDA0003645299810000042
式中
Figure FDA0003645299810000043
为离散网格i的测量值，ci为离散网格i的计算值。
CN202210527841.2A 2022-05-16 2022-05-16 一种基于地质统计学的河流污染溯源方法 Active CN114943194B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210527841.2A CN114943194B (zh) 2022-05-16 2022-05-16 一种基于地质统计学的河流污染溯源方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210527841.2A CN114943194B (zh) 2022-05-16 2022-05-16 一种基于地质统计学的河流污染溯源方法

Publications (2)

Publication Number Publication Date
CN114943194A true CN114943194A (zh) 2022-08-26
CN114943194B CN114943194B (zh) 2023-04-28

Family

ID=82906336

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210527841.2A Active CN114943194B (zh) 2022-05-16 2022-05-16 一种基于地质统计学的河流污染溯源方法

Country Status (1)

Country Link
CN (1) CN114943194B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117871171A (zh) * 2024-03-13 2024-04-12 太原市水利勘测设计院 一种河道淤泥取样装置

Citations (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104103005A (zh) * 2014-08-05 2014-10-15 环境保护部华南环境科学研究所 一种有限条件下突发性水环境事件污染源的溯源方法
CN105739951A (zh) * 2016-03-01 2016-07-06 浙江工业大学 一种基于gpu的l1最小化问题快速求解方法
CN105956664A (zh) * 2016-04-27 2016-09-21 浙江大学 一种河流点源突发污染事故溯源方法
CN106228007A (zh) * 2016-07-19 2016-12-14 武汉大学 突发事件污染源追溯方法
CN107895081A (zh) * 2017-11-14 2018-04-10 武汉大学 多个瞬时污染源的溯源方法
CN109084840A (zh) * 2018-08-16 2018-12-25 天狼联盟材料科技研究（广东）有限公司 一种基于物联网的河涌水域污染监控和分段管理方法
CN109190280A (zh) * 2018-09-18 2019-01-11 东北农业大学 一种基于核极限学习机替代模型的地下水污染源反演识别方法
CN109886830A (zh) * 2019-01-02 2019-06-14 同济大学 一种基于用户投诉信息的供水管网污染源追踪定位方法
CN111091082A (zh) * 2019-12-09 2020-05-01 北京空间机电研究所 一种基于高分辨率遥感数据的流域污染溯源方法
CN111272960A (zh) * 2020-02-20 2020-06-12 中国环境科学研究院 一种同位素和测年相结合的浅层地下水硝酸盐源解析方法
CN111667116A (zh) * 2020-06-08 2020-09-15 哈尔滨工业大学(深圳)(哈尔滨工业大学深圳科技创新研究院) 一种地下水污染源辨识方法及系统
CN112233247A (zh) * 2020-09-29 2021-01-15 广东新禾道信息科技有限公司 大气污染防治网格化监控设备及其使用方法
CN112307602A (zh) * 2020-10-13 2021-02-02 河海大学 一种地下水污染源信息和水力渗透系数场联合反演的方法
CN113128129A (zh) * 2021-05-07 2021-07-16 大连理工大学 一种突发水污染正逆耦合溯源方法及系统
CN113267607A (zh) * 2021-05-11 2021-08-17 浙江大学 一种场地有机污染物迁移过程的特征参数识别系统

Patent Citations (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104103005A (zh) * 2014-08-05 2014-10-15 环境保护部华南环境科学研究所 一种有限条件下突发性水环境事件污染源的溯源方法
CN105739951A (zh) * 2016-03-01 2016-07-06 浙江工业大学 一种基于gpu的l1最小化问题快速求解方法
CN105956664A (zh) * 2016-04-27 2016-09-21 浙江大学 一种河流点源突发污染事故溯源方法
CN106228007A (zh) * 2016-07-19 2016-12-14 武汉大学 突发事件污染源追溯方法
CN107895081A (zh) * 2017-11-14 2018-04-10 武汉大学 多个瞬时污染源的溯源方法
CN109084840A (zh) * 2018-08-16 2018-12-25 天狼联盟材料科技研究（广东）有限公司 一种基于物联网的河涌水域污染监控和分段管理方法
CN109190280A (zh) * 2018-09-18 2019-01-11 东北农业大学 一种基于核极限学习机替代模型的地下水污染源反演识别方法
CN109886830A (zh) * 2019-01-02 2019-06-14 同济大学 一种基于用户投诉信息的供水管网污染源追踪定位方法
CN111091082A (zh) * 2019-12-09 2020-05-01 北京空间机电研究所 一种基于高分辨率遥感数据的流域污染溯源方法
CN111272960A (zh) * 2020-02-20 2020-06-12 中国环境科学研究院 一种同位素和测年相结合的浅层地下水硝酸盐源解析方法
CN111667116A (zh) * 2020-06-08 2020-09-15 哈尔滨工业大学(深圳)(哈尔滨工业大学深圳科技创新研究院) 一种地下水污染源辨识方法及系统
CN112233247A (zh) * 2020-09-29 2021-01-15 广东新禾道信息科技有限公司 大气污染防治网格化监控设备及其使用方法
CN112307602A (zh) * 2020-10-13 2021-02-02 河海大学 一种地下水污染源信息和水力渗透系数场联合反演的方法
CN113128129A (zh) * 2021-05-07 2021-07-16 大连理工大学 一种突发水污染正逆耦合溯源方法及系统
CN113267607A (zh) * 2021-05-11 2021-08-17 浙江大学 一种场地有机污染物迁移过程的特征参数识别系统

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
姜继平，董芙嘉，刘仁涛， 袁一星: "基于河流示踪实验的Bayes污染溯源:算法参数、影响因素及频率法对比", 《中国环境科学》 *
石宝山 等: "基于GPU 加速的污染物输移高分辨率数值模型", 《水动力学研究与进展》 *
龙玉桥，崔婷婷，李伟，吴春勇，李砚阁: ""地质统计学法在地下水污染溯源中的应用及参数敏感性分析"", 《水利学报》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117871171A (zh) * 2024-03-13 2024-04-12 太原市水利勘测设计院 一种河道淤泥取样装置
CN117871171B (zh) * 2024-03-13 2024-05-24 太原市水利勘测设计院 一种河道淤泥取样装置

Also Published As

Publication number Publication date
CN114943194B (zh) 2023-04-28

Similar Documents

Publication Publication Date Title
Liao et al. Deep learning for air quality forecasts: a review
CN106650825B (zh) 一种机动车尾气排放数据融合系统
Shi et al. Surface modelling of soil pH
CN112785450B (zh) 一种土壤环境质量分区方法及系统
CN106372402A (zh) 一种大数据环境下模糊区域卷积神经网络的并行化方法
CN107273995A (zh) 空气质量预报方法
CN109541172A (zh) 土壤属性值的计算方法及装置
Barros et al. Machine learning for whole-building life cycle assessment: A systematic literature review
CN115330153A (zh) 一种重金属污染土壤治理修复决策方法
CN105787184A (zh) 一种基于pm2.5的大气气溶胶光学厚度估计方法
Gao et al. A two-point machine learning method for the spatial prediction of soil pollution
CN108764527B (zh) 一种土壤有机碳库时空动态预测最优环境变量筛选方法
Lark Multi-objective optimization of spatial sampling
CN114943194A (zh) 一种基于地质统计学的河流污染溯源方法
CN110321528B (zh) 一种基于半监督地理空间回归分析的高光谱影像土壤重金属浓度评估方法
Bai et al. Groundwater contamination source identification using improved differential evolution Markov chain algorithm
Haas Statistical assessment of spatio-temporal pollutant trends and meteorological transport models
CN117669270A (zh) 一种对微站组网数据进行时空一致性订正的方法
CN117953398A (zh) 一种基于无人机巡河技术的水域状况智能监测方法及系统
CN1327376C (zh) 基于支持向量机的软测量仪表建模方法
Hou et al. Bidirectional machine learning–assisted sensitivity-based stochastic searching approach for groundwater DNAPL source characterization
CN116842385B (zh) 一种基于履带车振动特征的lstm路面不平度识别方法
CN116645001B (zh) 一种基于多维度数据分析的金属矿山环境评价方法及装置
CN117973125B (zh) 基于人工智能和大数据的污染源清单反演方法、系统及应用
Guo et al. Research on Optimum Algorithm of Charging Pile Location for New Energy Electric Vehicle

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant