CN114970302B - 一种基于地下水监测***的区域地下水情预测方法 - Google Patents
一种基于地下水监测***的区域地下水情预测方法 Download PDFInfo
- Publication number
- CN114970302B CN114970302B CN202210500900.7A CN202210500900A CN114970302B CN 114970302 B CN114970302 B CN 114970302B CN 202210500900 A CN202210500900 A CN 202210500900A CN 114970302 B CN114970302 B CN 114970302B
- Authority
- CN
- China
- Prior art keywords
- groundwater
- prediction
- observation
- water
- prediction result
- 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.)
- Active
Links
- 239000003673 groundwater Substances 0.000 title claims abstract description 154
- 238000012544 monitoring process Methods 0.000 title claims abstract description 48
- 238000000034 method Methods 0.000 title claims abstract description 18
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims abstract description 134
- 238000005086 pumping Methods 0.000 claims abstract description 42
- 238000009826 distribution Methods 0.000 claims description 49
- 238000003809 water extraction Methods 0.000 claims description 26
- 239000002245 particle Substances 0.000 claims description 18
- 238000009499 grossing Methods 0.000 claims description 17
- 238000004364 calculation method Methods 0.000 claims description 16
- 230000009466 transformation Effects 0.000 claims description 16
- 238000004422 calculation algorithm Methods 0.000 claims description 12
- 238000000342 Monte Carlo simulation Methods 0.000 claims description 10
- 238000012417 linear regression Methods 0.000 claims description 10
- 238000004458 analytical method Methods 0.000 claims description 9
- 238000012549 training Methods 0.000 claims description 9
- 230000008878 coupling Effects 0.000 claims description 6
- 238000010168 coupling process Methods 0.000 claims description 6
- 238000005859 coupling reaction Methods 0.000 claims description 6
- 238000001514 detection method Methods 0.000 claims description 5
- 239000002349 well water Substances 0.000 claims description 5
- 235000020681 well water Nutrition 0.000 claims description 5
- 238000010276 construction Methods 0.000 claims description 4
- 239000011159 matrix material Substances 0.000 claims description 4
- 238000005457 optimization Methods 0.000 claims description 4
- 238000001556 precipitation Methods 0.000 claims description 4
- 238000010835 comparative analysis Methods 0.000 claims description 3
- 238000003708 edge detection Methods 0.000 claims description 3
- 238000007619 statistical method Methods 0.000 claims description 3
- 238000012216 screening Methods 0.000 claims description 2
- 238000011002 quantification Methods 0.000 abstract description 2
- 230000006870 function Effects 0.000 description 34
- 238000004519 manufacturing process Methods 0.000 description 4
- 238000004088 simulation Methods 0.000 description 3
- 238000010801 machine learning Methods 0.000 description 2
- 238000004321 preservation Methods 0.000 description 2
- 238000013398 bayesian method Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 239000003086 colorant Substances 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000002950 deficient Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 230000002262 irrigation Effects 0.000 description 1
- 238000003973 irrigation Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000003911 water pollution Methods 0.000 description 1
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/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
- 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/06—Energy or water supply
-
- 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/02—Reliability analysis or reliability optimisation; Failure analysis, e.g. worst case scenario performance, failure mode and effects analysis [FMEA]
-
- 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)
- Business, Economics & Management (AREA)
- Theoretical Computer Science (AREA)
- Evolutionary Computation (AREA)
- General Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Economics (AREA)
- Strategic Management (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Water Supply & Treatment (AREA)
- Tourism & Hospitality (AREA)
- Public Health (AREA)
- General Business, Economics & Management (AREA)
- Marketing (AREA)
- Human Resources & Organizations (AREA)
- Artificial Intelligence (AREA)
- Primary Health Care (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Software Systems (AREA)
- Computer Hardware Design (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明提供了一种基于地下水监测***的区域地下水情预测方法，所述方法包括如下步骤：通过地下水观测***中的观测数据，结合实时监测数据，基于支持向量回归SVR，建立各个观测井抽水量预测函数和抽水量‑地下水水位预测函数，对地下水水位进行实时预测，解决了传统地下水水情的估计建模耗时长且需要大量水文地质参数，受限于复杂的地质条件、水情预测参数的不确定性与数学建模负担等因素，尚未形成较精确的区域地下水水情量化预测工具，无法细化到区域各部分地下水水情的精准预测，不能很好地服务于城市、镇、乡村区域地下水管理的问题。
Description
技术领域
本发明属于计算机机器学习技术领域，具体涉及一种基于地下水监测***的区域地下水情预测方法。
背景技术
我国大部分地区在气候变迁与水资源分布不均的影响下，水资源的运用越显严峻，由于水文地质条件复杂，有关地下水水文资料较为匮乏。由于地下水具有水量稳定、水质好等特点，是农业灌溉、企业生产和城市居民生活的重要水源之一，尤其是在地表缺水的干旱和半干旱地区，地下水常常被作为当地的主要用水源。地下水资源在区域供水问题中，起到了至关重要的作用，若没有对下水资源进行妥善管理，地下水资源的过度使用容易造成地层下陷、塌陷等问题，造成局部地区水资源衰减，甚至地下水污染问题，带来严重的环境问题。2020年，整个华北平原都由于近年来降雨量普遍较少，都出现了地下水位下降的趋势。因此有必要加强地下水资源的合理使用，实施科学有效的监管，保护地下水的安全，促进地下水资源的可持续开发利用。
传统地下水水情的估计都是以各种水力试验与数值模式推算地下水***可用水量，建模耗时长且需要大量水文地质参数，受限于复杂的地质条件、水情预测参数的不确定性与数学建模负担等因素，尚未形成较精确的区域地下水水情量化预测工具，无法细化到区域各部分地下水水情的精准预测，不能很好地服务于区域供水管理，尤其是根据单一化的预警参数模型对整个区域进行停水、供水、限制供水量、供水时间做出部署，会为居民的生活造成不便、甚至会影响企业生产或工程进度安排。
发明内容
因此，本发明要解决的技术问题在于提供一种基于地下水监测***的区域地下水情预测方法，能够解决了传统地下水水情的估计建模耗时长且需要大量水文地质参数，受限于复杂的地质条件、水情预测参数的不确定性与数学建模负担等因素，尚未形成较精确的区域地下水水情量化预测工具，无法细化到区域各部分地下水水情的精准预测，不能很好地服务于城市、镇、乡村区域地下水管理的问题。
为了解决上述问题，本发明提供了一种基于地下水监测***的区域地下水情预测方法，所述方法包括如下步骤：
S1：通过地下水观测***中的观测数据，结合实时监测数据，基于支持向量回归SVR，建立各个观测井抽水量预测函数和抽水量-地下水水位预测函数，对地下水水位进行实时预测；
S2：采用马尔可夫-蒙特卡洛法计算抽水井与其对应的观测井的抽水-水位泄降关系，利用观测网实时监测的观测数据，同时计算多个地下水水情参数，基于贝叶斯理论建立的概率分布，通过对样本进行迭代筛选，获取粒子来近似后验概率分布，同时计算各个观测井水位预测结果和观测井抽水量预测结果，避免各个预测参数之间存在耦合；
S3：根据观测井的坐标和分布情况，划定覆盖全部观测井的预测区域，预测区域覆盖整个监测网区域，对整个预测区域进行网格划分，生成n×n个网格点；
S4：构造监测网区域的泰森多边形，结合欧几里德距离变换对应坐标集合Ue并入到观测井坐标集合Um当中，形成二维数组Um+e再合并泰森多边形顶点坐标集合UV，构成二维数组表示集合观测井坐标、欧几里德距离变换坐标和泰森多边形顶点坐标的坐标点集合Um+e+V；
利用观测井水位预测结果，根据观测区域构造的泰森多边形对覆盖整个监测网的预测区域进行模拟，对生成的n×n个网格点赋值，设定泰森多边形中观测井i的泰森多边形Vi区域网格点赋值，为观测井i的地下水位预测结果，同时将观测井抽水量预测结果赋值给Vi区域网格点，生成该区域抽水量的初步预测结果，同时生成该区域地下水水位的初步预测结果；
S5：通过观测井坐标、欧几里德距离变换坐标以及泰森多边形顶点坐标构建包含所有散点的凸多边形点集Um+e+V，将凸多边形分割成互不相交的三角形，对相邻三角形进行调整获得带权三角剖分，通过三角形内插值计算，根据坐标点集中各点对应的初步预测结果进行修正，对生成的地下水水位-抽水量区域分布初步预测结果进行实时修正；
S6：采用高斯平滑函数对修正后的地下水水位-抽水量分布区域水情预测结果进行平滑，得到最终预测结果，采用若梯度双阈值检测进一步修正预测结果，对采用支持向量回归算法的地下水水情预测结果和采用马尔可夫-蒙特卡洛法算法的地下水水情预测结果进行统计分析，根据对比残差统计学分布比较分析预测结果的精确度，创建一个保存结果的数组，确定保存优化预测模型作为稳定的地下水水位预测***模型。
可选的，S1具体包括如下步骤：
S11：对应不同各观测井坐标按照时间序列建立一个训练样本集表示抽水量的历史观测数据，对应不同各观测井坐标按照时间序列建立一个训练样本集/>表示地下水水位的历史观测数据；
S12：在样本空间中训练模型函数，在如下约束条件yi∈{-1,+1}，松弛变量ξi≥0，松弛变量，ηi≥0，i＝0,1,2…t，i表示第i个时间点；t表示观测时间点总数；损失函数为：
其中，C为>0的常数；ε给定决策面位置；μ表示决策面的法向量；ξi和ηi为松弛变量；满足优化目标，得到降水量预测模型函数f1(xa，ti)，其中，ti表示时间，xa表示抽水量；
在如下约束条件得到地下水水位预测模型函数f2(xa，xb)，其中，xa表示抽水量；xb表示地下水水位；即对未来时刻的抽水量和地下水水位进行实时预测；
S13：对抽水量观测数据和地下水水位的历史观测数据/>进行线性回归，构建多种的抽水量-地下水水位线性回归模型，计算得到残差，进行对比分析，将线性回归模型函数的预测结果与支持向量回归预测的地下水水位结果进行对比，分析预测的精度，得到优化的模型。
可选的，步骤S2中，基于马尔可夫-蒙特卡洛法计算抽水井与其对应的观测井的抽水-水位泄降关系，避免各个预测参数之间存在耦合，具体包括如下步骤：
S21：基于贝叶斯理论建立水情参数连续计算的概率模型，水情观测历史数据可以同时计算包括抽水量和地下水水位的多个水情参数的概率分布；
设定t0时刻及t0时刻前各个观测井及收集的水情观测数据为Xt，后验概率P(xa,xb,t|Xt)可以表示为：
其中，xa表示抽水量；xb表示地下水水位；t0预测起始时间；t
为预测时间；
S22：取样本简化模型φt(xa，xb，t|Xt)符合q分布，q(·)为先验分布P(φt|φt-1)的转移矩阵，i表示符合q分布的第i个粒子，记作t表示时刻，由i＝1，…，N循环迭代；
后验概率分布可以近似归纳为：
其中，t时刻第i个粒子的第一权系数ⅴt i，总数可以归一化为1，N为粒子总数，δ表示等比例概率；
基于更新后的观测数据Xt，参与到水情参数的计算中，所有粒子权重重新标准化为：
其中，L(Xt|φt)＝L(Xt|xa，xb，t)作为第二权系数，通过若干带有权重的粒子，计算抽水井与其对应的观测井的抽水-水位泄降关系，构建区域地下水情预测模型。
可选的，S3具体包括如下步骤：
S31：根据观测井的坐标和分布情况划定覆盖全部观测井的预测区域，预测区域覆盖整个监测网区域，基于由监测网提供的实时监测数据，对预测区域进行网格划分和观测井坐标数组的建立；
S32：根据预测区域范围划分的网格，创建网格点区域中心的二维数组以横纵坐标(xn，yn)为中心的n×n区域Un＝{(x1,y1),(x2,y2),…,(xn-1,yn-1),(xn,yn)}，将数组中所有点初始化为单位矩阵函数，其中l和h分别为x向与y方向相邻网格点之间的间距，xn＝l×n，yn＝h×n；
S33：构建一个对应观测井位置坐标的二维数组Um＝{(x1,y1),(x2,y2),…,(xm-1,ym-1),(xm,ym)}，(xm，ym)为Um中第m个观测井对应的水情预测结果。
可选的，S4具体包括如下步骤：
S41：构造监测网区域的泰森多边形：根据N个观测井坐标对预测区域构造泰森多边形，以相邻观测井作中垂线生成泰森多边形，即VN图；
利用观测井水位预测结果，根据观测区域构造的VN图对覆盖整个监测网的预测区域进行模拟，对生成的n×n个网格点赋值，设定泰森多边形当中观测井i的泰森多边形Vi区域网格点赋值，为观测井i的地下水位预测结果，同时将观测井抽水量预测结果赋值给Vi区域网格点；
S42：以Df(xi，yi)表示坐标(xi，yi)的欧几里德距离变换计算，以f(xj，yj)表示坐标(xj，yj)的消耗函数，(xi,yi)∈Un，(xj,yj)∈Un，Df(xi,yi)＝min(((xi,yi)-(xj,yj))2+f(xj,yj))遍历n×n区域Un＝{(x1,y1),(x2,y2),…,(xn-1,yn-1),(xn,yn)}，计算得到欧几里德距离变换对应坐标集合Ue并入到观测井坐标集合Um当中，形成二维数组Um+e，再合并泰森多边形顶点坐标集合UV，构成二维数组表示集合观测井坐标、欧几里德距离变换对应坐标集合Ue坐标和泰森多边形顶点坐标的坐标点集合Um+e+V。
可选的，S5具体包括如下步骤：
S51：将生成的二维数组，对应观测井位置坐标、欧几里德距离变换对应坐标集合Ue坐标以及泰森多边形顶点坐标的二维数组Um+e+V，每个坐标点对应Um+e+V中的观测井相应的水情预测结果，Um+e+V当中每个坐标点都对应一个时间序列的地下水水情预测结果，由Um+e+V中的点作为端点构成三角剖分网络的封闭线段；
S52：Um+e+V中的点作为端点构成三角剖分网络凸多边形，在网格点Un中包络的点集为Uk＝{(x1,y1),(x2,y2),…,(xk-1,yk-1),(xk,yk)}，其任意两点形权重为ωij，将凸多边形分割成若干个互不相交的三角形，生成带权三角剖三角网；
S53：构造包含Um+e+V全部坐标点区域的三角形剖分，递归得到三角网的三角剖分，对三角网执行三角形内插值计算，对三角网内每个坐标点集Uk对应的时间序列抽水量预测结果和地下水水位预测结果进行进一步修正，对生成的初步预测结果进行优化更新。
可选的，S6具体包括如下步骤：
S61：Un＝{(x1,y1),(x2,y2),…,(xn-1,yn-1),(xn,yn)}区域中的网格点各中心(xi,yi)构造高斯核函数如下：
其中，σ为高斯核生成器中高斯分布的标准差，表征预测结果数据的离散程度；
S62：对修正后的抽水量预测结果和地下水水位预测结果进行高斯平滑，采用sobel函数对预测结果横向和纵向梯度进行计算，采用canny函数设置双阈值进行边缘检测，利用弱梯度值进行双阈值检测，对整个区域的预测结果进一步修正，得到最终预测结果；
S63：基于实测数据，生成的区域地下水水位分布模拟结果，对生成的预测结果采用高斯平滑函数进行平滑，进一步修正后得到最终预测结果，对采用SVR的回归算法的地下水水情预测结果和基于马尔可夫-蒙特卡洛法算的地下水水情预测结果进行模拟区域分布；
S64：基于不同的算法预测出的结果进行对比，统计分析预测的精度和计算效率，并且加入监测期间气温与降雨或/和降雪、地表沉降信息，分析不同条件对最终预测结果的影响，根据精度和计算效率分析，保存具备鲁棒性的区域地下水水情预测模型。
有益效果
本发明的实施例中所提供的一种基于地下水监测***的区域地下水情预测方法，可以实现整个监测区域内整体地下水水情的预测，不受观测井位置和观测井分布情况的影响，无论观测资料是否充足，只要在若干离散分布在监测区域围内的有限个观测井和抽水井(群)取得观测数据就可以进行预测，且可以通过现有的观测资料针对多个不同的预测方法同时进行精度和计算效率的对比分析，细化到预测区域各部分区域的地下水水位分布，实时地反映预测区域各部分地下水水情。此外，本发明可以实现通过构建抽水量预测函数建立抽水量-地下水水位泄降关系函数预测地下水水情，也可以减少抽水量和地下水水位之间的耦合，基于马尔可夫-蒙特卡洛法同时建立抽水量-地下水水位泄降关系预测函数。根据实际案例预测结果的精度分析，保存具备鲁棒性的区域地下水水情预测模型。解决了传统地下水水情的估计建模耗时长且需要大量水文地质参数，受限于复杂的地质条件、水情预测参数的不确定性与数学建模负担等因素，尚未形成较精确的区域地下水水情量化预测工具，无法细化到区域各部分地下水水情的精准预测，不能很好地服务于区域供水管理，尤其是根据单一化的预警参数模型对整个区域进行停水、供水、限制供水量、供水时间做出部署，会为居民的生活造成不便、甚至会影响企业生产或工程进度安排的问题。
附图说明
图1为本发明实施例的基于地下水监测***的区域地下水情预测的流程图；
图2为本发明实施例的地下水水位观测数据的区域模拟分布示意图；
具体实施方式
根据本发明的实施例，一种基于地下水监测***的区域地下水情预测方法，请参照图1，包括如下步骤：
S1：按照时间序列收集通过各观测井(群)水位观测设备采集地下水水位观测信息，观测井附近的抽水量观测信息，除此之外还需要收集气温与降雨(降雪)信息，并根据各观测井的坐标对提取的观测数据进行收集和整理。
利用已有的观测数据，通过支持向量回归(Support Vector Regression，SVR)模型构建抽水量预测函数和抽水量-地下水水位预测函数，从而对未来的地下水水位进行预测。
S1具体包括如下步骤：
S11：对应不同各观测井坐标按照时间序列建立一个训练样本集表示抽水量的历史观测数据，对应不同各观测井坐标按照时间序列建立一个训练样本集/>表示地下水水位的历史观测数据。
S12：在样本空间中训练模型函数，定一个宽度ε给定决策面位置，μ为决策面的法向量，决定了决策面的方向，引入松弛变量ξi和ηi，在如下约束条件yi∈{-1,+1}松弛变量ξi≥0，松弛变量ηi≥0，i＝0,1,2…t，i表示第i个时间点，t表示观测时间点总数，损失函数如下：
其中，C为>0的常数，ε给定决策面位置，μ表示决策面的法向量，ξi和ηi为松弛变量，满足优化目标，得到降水量预测模型函数f1(xa，ti)，其中ti表示时间，xa表示抽水量，同样在如下约束条件得到地下水水位预测模型函数f2(xa，xb)，xa表示抽水量，xb表示地下水水位，即可对未来时刻的抽水量和地下水水位进行实时预测。
S13：对抽水量观测数据地下水水位的历史观测数据/>进行线性回归，构建多种的抽水量-地下水水位线性回归模型(如最小二乘、多项式回归、对数线性回归)，计算出残差进行对比分析。将线性回归模型函数的预测结果同S12当中利用SVR预测的地下水水位结果进行对比，分析预测的精度。
S2：基于马尔可夫-蒙特卡洛法估计抽水井(群)与其对应的观测井的抽水-水位泄降关系，构建区域地下水情预测模型。
S2具体包括如下步骤：
S21：基于贝叶斯理论(Bayesian method)建立水情参数连续计算的概率模型，水情观测历史数据可以同时估计包括抽水量和地下水水位在内的多个水情参数的概率分布。令t0时刻及t0时刻前各个观测井及附近收集的水情观测数据为Xt，后验概率P(xa，xb，t|Xt)可以表示为如下：
其中，xa表示抽水量，xb表示地下水水位，t0预测起始时间，t为预测时间。
S22：随着大量带有权重的粒子来近似后验分布，经过迭代可以计算出估计结果：
取样本简化模型φt(xa，xb，t|Xt)符合q分布，q(·)通常被简化定义为先验分布P(φt|φt-1)的转移矩阵，i表示符合q分布的第i个粒子，记作t表示时刻，由i＝1，…，N开始循环迭代。后验概率分布可以近似归纳为如下：
其中，t时刻第i个粒子的权系数ⅴt i，其总数可以归一化为1，N为粒子总数，δ表示等比例概率。
随着每一步更新的观测数据Xt参与到水情参数的估计当中，所有粒子权重重新标准化为
其中，L(Xt|φt)＝L(Xt|xa，xb，t)作为新的权系数，通过大量带有权重的粒子，估计抽水井(群)与其对应观测井的抽水-水位泄降关系，构建区域地下水情预测模型。
S3：根据观测井的坐标和分布情况划定覆盖全部观测井(群)的预测区域，预测区域覆盖整个监测网区域，对该区域进行网格划分，为后续的泰森多边形结合欧几里德距离变换的初步预测、该区域的插值计算和模拟平滑提供初始条件，划分成n×n个网格，设置网格的横纵间距l和h，并且统一以单位矩阵函数进行初始化为。
S3具体包括如下步骤：
S31：根据观测井的坐标和分布情况划定覆盖全部观测井(群)的预测区域，预测区域覆盖整个监测网区域，尽可能全部利用由监测网提供的实时监测数据，对预测区域进行网格划分和观测井坐标数组的建立；
S32：根据预测区域范围划分的网格，创建网格点区域中心的二维数组以横纵坐标(xn，yn)为中心的n×n区域Un＝{(x1,y1),(x2,y2),…,(xn-1,yn-1),(xn,yn)}，将数组中所有点初始化为单位矩阵函数，其中l和h分别为x向与y方向相邻网格点之间的间距，xn＝l×n，yn＝h×n；
S33：构建一个对应观测井位置坐标的二维数组Um＝{(x1,y1),(x2,y2),…,(xm-1,ym-1),(xm,ym)}，(xm，ym)为Um中第m个观测井对应的水情预测结果；
经过本步骤，完成了预测区域网格的划分和观测井水情预测情况的数组建立，为后续初步预测、修正计算和平滑计算提供了条件，一次保证即便在监测网分布稀疏，缺乏充足的观测资料的情况下，依然可以对预测区域内任意位置的降水量和地下水水位进行和预测。
S4：对预测区域模拟初步预测结果：以相邻观测井作中垂线构造Voronoi多边形，N个观测井，即VN图，以观测井坐标为中心划分成n个泰森多边形，增加一个观测井，生成N+1个泰森多边形，即VN+1图，结合欧几里德距离变换对应坐标集合Ue并入到观测井坐标集合Um当中，形成新的二维数组Um+e再合并泰森多边形顶点坐标集合UV，构成二维数组表示集合观测井坐标、欧几里德距离变换的坐标和泰森多边形顶点坐标的坐标点集合Um+e+V，从而对整个预测区域进行模拟，得到整个区域的初步预测结果；
对整个预测区域进行模拟，得到整个区域的地下水水情分布的初步预测结果。
S4具体包括如下步骤：
S41：构造监测网区域的泰森多边形：根据N个观测井坐标对预测区域构造泰森多边形，以相邻观测井作中垂线生成泰森多边形，即VN图。如果新增加一个观测井，则生成新的VN图；
利用观测井水位预测结果，根据观测区域构造的VN图对覆盖整个监测网的预测区域进行模拟，对S3中生成的n×n个网格点赋值，令泰森多边形当中观测井i的泰森多边形Vi区域网格点赋值为观测井i的地下水位预测结果，同时将观测井附近抽水量预测结果赋值给Vi区域网格点，
S42：以Df(xi，yi)表示坐标(xi，yi)的欧几里德距离变换计算，以f(xj，yj)表示坐标(xj，yj)的消耗函数，(xi,yi)∈Un，(xj,yj)∈Un，Df(xi,yi)＝min(((xi,yi)-(xj,yj))2+f(xj,yj))遍历n×n区域Un＝{(x1,y1),(x2,y2),…,(xn-1,yn-1),(xn,yn)}，算得欧几里德距离变换对应坐标集合Ue并入到观测井坐标集合Um当中，形成新的二维数组Um+e，再合并泰森多边形顶点坐标集合UV，构成二维数组表示集合观测井坐标、欧几里德距离变换的坐标和泰森多边形顶点坐标的坐标点集合Um+e+V。
步骤S1和步骤S2当中地下水水情预测结果，根据观测区域构造的泰森多边形当中不同的泰森多边形进行赋值，以观测井坐标为中心泰森多边形赋值为该观测井的观测水位或预测水位，结合欧几里德距离变换对应坐标集合Ue形成新的二维数组Um+e再合并泰森多边形顶点坐标集合UV，构成二维数组表示集合观测井坐标、欧几里德距离变换对应坐标和泰森多边形顶点坐标的坐标点集合Um+e+V，从而对整个预测区域进行模拟，得到整个区域的初步预测结果；
S5：通过观测井坐标、欧几里德距离变换对应坐标和泰森多边形顶点坐标构建包含所有坐标点Um+e+V的带权三角剖分，对相邻三角形进行调整，通过最大凸壳点集合的带权三角剖对初步预测结果进行实时修正，由观测井构成两个具有共同边的三角形合成一个多边形，将所有坐标点包围起来，根据每个观测井预测结果对最大凸多边形点集合相应的水情预测结果进行逐点修正，对步骤S4中生成的初步预测结果进行实时修正；
S5具体包括如下步骤：
S51：将步骤S4中生成的二维数组，对应观测井位置坐标、欧几里德距离变换对应坐标和泰森多边形顶点坐标的二维数组Um+e+V，每个坐标点对应Um+e+V中的一个步骤S1和步骤S2当中观测井相应的水情预测结果，Um+e+V当中每个坐标点都对应一个时间序列的地下水水情预测结果，由Um+e+V中的点作为端点构成三角剖分网络的封闭线段；
S52：Um+e+V中的点作为端点构成三角剖分网络凸多边形，在网格点Un中包络的点集为Uk＝{(x1,y1),(x2,y2),…,(xk-1,yk-1),(xk,yk)}，其任意两点形的权重为ωij，将凸多边形分割成若干个互不相交的三角形，生成带权三角剖三角网。
S53：构造包含Um+e+V全部坐标点区域的三角形剖分，递归得到三角网对应的权值函数值，对三角网执行三角形内插值计算，对三角网内每个坐标点集Uk对应的时间序列抽水量预测结果和地下水水位预测结果进行进一步修正，对从而对步骤S4中生成的初步预测结果进行优化更新。
S6：采用高斯平滑函数对步骤S5的预测结果进行平滑，得到最终预测结果；对步骤S1采用不同的回归算法的地下水水情预测结果和步骤S2中基于马尔可夫-蒙特卡洛法算的的地下水水情预测结果进行统计分析，分析预测结果的精确度。
S61：Un＝{(x1,y1),(x2,y2),…,(xn-1,yn-1),(xn,yn)}区域中的网格点各中心(xi,yi)构造高斯核函数σ为高斯核生成器中高斯分布的标准差，表征预测结果数据的离散程度；
S62：对步骤S5当中修正后的抽水量预测结果和地下水水位预测结果进行高斯平滑，采用sobel函数对预测结果横向和纵向梯度进行计算，采用canny函数设置双阈值进行边缘检测，利用弱梯度值进行双阈值检测，对整个区域的预测结果进一步修正，得到最终预测结果。
S63：根据实测数据采用步骤S4～步骤S6生成的区域地下水水位分布模拟结果，对步骤S5生成的预测结果采用高斯平滑函数进行平滑，进一步修正后得到最终预测结果，对步骤S1采用支持向量回归等各种不同的回归算法的地下水水情预测结果和步骤S2中基于马尔可夫-蒙特卡洛法算的的地下水水情预测结果采用步骤S4～S6模拟区域分布。
S64：步骤S1和步骤S2当中可以载入不同的地下水水情预测算法，针对不同的算法预测出的结果进行对比，统计分析预测的精度和计算效率等，并且加入监测期间气温与降雨和/或降雪、地表沉降等信息，分析这些细节条件对最终预测结果的影响，根据精度和计算效率分析，保存具备鲁棒性的区域地下水水情预测模型。
如图2所示，以某城市施工区域、某日地下水位观测数据经过步骤S4～S6实现该区域的地下水位分布的数值模拟过程，得到该区域150m×40m范围内地下水水位分布，圆点表示观测井位置，颜色由浅到深表示观测井水位高度由低到高(-5～25m)，经步骤S4、S5得到经三角网内插值后的地下水水位分布，最终得到经步骤S6平滑后的区域地下水水位分布。
本发明区别于以往传统地下水水情数值预测模型或机器学习方法采用离散分布的单观测井进行预测，根据经验公式估计全区域的地下水水情，某些区域缺乏观测资料的情况下，这种估计对于整个区域的水情分布情况会存在很大的误差，使得地下水水情预测不够准确，尤其是根据单一化的预警参数模型对整个区域进行停水、供水、限制供水量、供水时间做出部署，会为居民的生活造成不便、甚至会影响企业生产或工程进度安排。本发明可以实现整个监测区域内整体地下水水情的预测，不受观测井位置和观测井分布情况的影响，无论观测资料是否充足，只要在若干离散分布在监测区域围内的有限个观测井和抽水井(群)取得观测数据就可以进行预测，且可以通过现有的观测资料针对多个不同的预测方法同时进行精度和计算效率的对比分析，细化到预测区域各部分区域的地下水水位分布，实时地反映预测区域各部分地下水水情。此外，本发明可以实现通过构建抽水量预测函数建立抽水量-地下水水位泄降关系函数预测地下水水情，也可以减少抽水量和地下水水位之间的耦合，基于马尔可夫-蒙特卡洛法同时建立抽水量-地下水水位泄降关系预测函数。根据实际案例预测结果的精度分析，保存具备鲁棒性的区域地下水水情预测模型。
本领域的技术人员容易理解的是，在不冲突的前提下，上述各有利方式可以自由地组合、叠加。
Claims (7)
1.一种基于地下水监测***的区域地下水情预测方法，其特征在于，所述方法包括如下步骤：
S1：通过地下水观测***中的观测数据，结合实时监测数据，基于支持向量回归，建立各个观测井抽水量预测函数和抽水量-地下水水位预测函数，对地下水水位进行实时预测；
S2：采用马尔可夫-蒙特卡洛法计算抽水井与其对应的观测井的抽水-水位泄降关系，利用观测网实时监测的观测数据，同时计算多个地下水水情参数，基于贝叶斯理论建立的概率分布，通过对样本进行迭代筛选，获取粒子来近似后验概率分布，同时计算各个观测井水位预测结果和观测井抽水量预测结果，避免各个预测参数之间存在耦合；
S3：根据观测井的坐标和分布情况，划定覆盖全部观测井的预测区域，预测区域覆盖整个监测网区域，对整个预测区域进行网格划分，生成n×n个网格点；
S4：构造监测网区域的泰森多边形，结合欧几里德距离变换对应坐标集合Ue并入到观测井坐标集合Um当中，形成二维数组Um+e再合并泰森多边形顶点坐标集合UV，构成二维数组表示集合观测井坐标、欧几里德距离变换坐标和泰森多边形顶点坐标的坐标点集合Um+e+V；
利用观测井水位预测结果，根据观测区域构造的泰森多边形对覆盖整个监测网的预测区域进行模拟，对生成的n×n个网格点赋值，设定泰森多边形中观测井i的泰森多边形Vi区域网格点赋值，为观测井i的地下水位预测结果，同时将观测井抽水量预测结果赋值给Vi区域网格点，生成该区域抽水量的初步预测结果，同时生成该区域地下水水位的初步预测结果；
S5：通过观测井坐标、欧几里德距离变换交点坐标以及泰森多边形顶点坐标构建包含所有散点的凸多边形点集Um+e+V，将凸多边形分割成互不相交的三角形，对相邻三角形进行调整获得带权三角剖分，通过三角形内插值计算，根据坐标点集中各点对应的初步预测结果进行修正，对生成的地下水水位-抽水量区域分布初步预测结果进行实时修正；
S6：采用高斯平滑函数对修正后的地下水水位-抽水量分布区域水情预测结果进行平滑，得到最终预测结果，采用若梯度双阈值检测进一步修正预测结果，对采用支持向量回归算法的地下水水情预测结果和采用马尔可夫-蒙特卡洛法算法的地下水水情预测结果进行统计分析，根据对比残差统计学分布比较分析预测结果的精确度，创建一个保存结果的数组，优化预测模型保存作为稳定的地下水水位预测***模型。
2.根据权利要求1所述的基于地下水监测***的区域地下水情预测方法，其特征在于，S1具体包括如下步骤：
S11：对应不同各观测井坐标按照时间序列建立一个训练样本集表示抽水量的历史观测数据，对应不同各观测井坐标按照时间序列建立一个训练样本集/>表示地下水水位的历史观测数据；
S12：在样本空间中训练模型函数，在如下约束条件yi∈{-1,+1}，松弛变量ξi≥0，松弛变量ηi≥0，i＝0,1,2…t，i表示第i个时间点；t表示观测时间点总数；损失函数为：
其中，C为>0的常数；ε给定决策面位置；μ表示决策面的法向量；ξi和ηi为松弛变量；满足优化目标，得到降水量预测模型函数f1(xa，ti)，其中，ti表示时间，xa表示抽水量；
在如下约束条件得到地下水水位预测模型函数f2(xa，xb)，其中，xa表示抽水量；xb表示地下水水位；即对未来时刻的抽水量和地下水水位进行实时预测；
S13：对抽水量观测数据和地下水水位的历史观测数据/>进行线性回归，构建多种的抽水量-地下水水位线性回归模型，计算得到残差，进行对比分析，将线性回归模型函数的预测结果与支持向量回归预测的地下水水位结果进行对比，分析预测的精度，优化预测模型。
3.根据权利要求1所述的基于地下水监测***的区域地下水情预测方法，其特征在于，步骤S2中，基于马尔可夫-蒙特卡洛法计算抽水井与其对应的观测井的抽水-水位泄降关系，避免各个预测参数之间存在耦合，具体包括如下步骤：
S21：基于贝叶斯理论建立水情参数连续计算的概率模型，水情观测历史数据可以同时计算包括抽水量和地下水水位的多个水情参数的概率分布；
设定t0时刻及t0时刻前各个观测井及收集的水情观测数据为Xt，后验概率P(xa,xb,t|Xt)可以表示为：
其中，xa表示抽水量；xb表示地下水水位；t0预测起始时间；t为预测时间；
S22：取样本简化模型φt(xa，xb，t|Xt)符合q分布，q(·)为先验分布P(φt|φt-1)的转移矩阵，i表示符合q分布的第i个粒子，记作t表示时刻，由i＝1，…，N循环迭代；
后验概率分布可以近似归纳为：
其中，t时刻第i个粒子的第一权系数总数可以归一化为1，/>N为粒子总数；δ表示等比例概率；
基于更新后的观测数据Xt，参与到水情参数的计算中，所有粒子权重重新标准化为：
其中，L(Xt|φt)＝L(Xt|xa，xb，t)作为第二权系数，通过若干带有权重的粒子，计算抽水井与其对应的观测井的抽水-水位泄降关系，构建区域地下水情预测模型。
4.根据权利要求1所述的基于地下水监测***的区域地下水情预测方法，其特征在于，S3具体包括如下步骤：
S31：根据观测井的坐标和分布情况划定覆盖全部观测井的预测区域，预测区域覆盖整个监测网区域，基于由监测网提供的实时监测数据，对预测区域进行网格划分和观测井坐标数组的建立；
S32：根据预测区域范围划分的网格，创建网格点区域中心的二维数组以横纵坐标(xn，yn)为中心的n×n区域Un＝{(x1,y1),(x2,y2),…,(xn-1,yn-1),(xn,yn)}，将数组中所有点初始化为单位矩阵函数，其中l和h分别为x向与y方向相邻网格点之间的间距，xn＝l×n，yn＝h×n；
S33：构建一个对应观测井位置坐标的二维数组Um＝{(x1,y1),(x2,y2),…,(xm-1,ym-1),(xm,ym)}，(xm，ym)为Um中第m个观测井对应的水情预测结果。
5.根据权利要求1所述的基于地下水监测***的区域地下水情预测方法，其特征在于，S4具体包括如下步骤：
S41：构造监测网区域的泰森多边形：根据N个观测井坐标对预测区域构造泰森多边形，以相邻观测井作中垂线生成泰森多边形，即VN图；
利用观测井水位预测结果，根据观测区域构造的VN图对覆盖整个监测网的预测区域进行模拟，对生成的n×n个网格点赋值，设定泰森多边形当中观测井i的泰森多边形Vi区域网格点赋值，为观测井i的地下水位预测结果，同时将观测井抽水量预测结果赋值给Vi区域网格点；
S42：以Df(xi，yi)表示坐标(xi，yi)的欧几里德距离变换计算，以f(xj，yj)表示坐标(xj，yj)的消耗函数，(xi,yi)∈Un，(xj,yj)∈Un，Df(xi,yi)＝min(((xi,yi)-(xj,yj))2+f(xj,yj))遍历n×n区域Un＝{(x1,y1),(x2,y2),…,(xn-1,yn-1),(xn,yn)}，计算得到欧几里德距离变换对应坐标集合Ue并入到观测井坐标集合Um当中，形成二维数组Um+e，再合并泰森多边形顶点坐标集合UV，构成二维数组表示集合观测井坐标、欧几里德距离变换坐标和泰森多边形顶点坐标的坐标点集合Um+e+V。
6.根据权利要求1所述的基于地下水监测***的区域地下水情预测方法，其特征在于，S5具体包括如下步骤：
S51：将生成的二维数组，对应观测井位置坐标、欧几里德距离变换坐标以及泰森多边形顶点坐标的二维数组Um+e+V，每个坐标点对应Um+e+V中的观测井相应的水情预测结果，Um+e+V当中每个坐标点都对应一个时间序列的地下水水情预测结果，由Um+e+V中的点作为端点构成三角剖分网络的封闭线段；
S52：Um+e+V中的点作为端点构成三角剖分网络凸多边形，在网格点Un中包络的点集为Uk＝{(x1,y1),(x2,y2),…,(xk-1,yk-1),(xk,yk)}，其任意两点形权重为ωij，将凸多边形分割成若干个互不相交的三角形，生成带权三角剖三角网；
S53：构造包含Um+e+V全部坐标点区域的三角形剖分，递归得到三角网的三角剖分，对三角网执行三角形内插值计算，对三角网内每个坐标点集Uk对应的时间序列抽水量预测结果和地下水水位预测结果进行进一步修正，对生成的初步预测结果进行优化更新。
7.根据权利要求1所述的基于地下水监测***的区域地下水情预测方法，其特征在于，S6具体包括如下步骤：
S61：Un＝{(x1,y1),(x2,y2),…,(xn-1,yn-1),(xn,yn)}区域中的网格点各中心(xi,yi)构造高斯核函数如下：
其中，σ为高斯核生成器中高斯分布的标准差，表征预测结果数据的离散程度；
S62：对修正后的抽水量预测结果和地下水水位预测结果进行高斯平滑，采用sobel函数对预测结果横向和纵向梯度进行计算，采用canny函数设置双阈值进行边缘检测，利用弱梯度值进行双阈值检测，对整个区域的预测结果进一步修正，得到最终预测结果；
S63：基于实测数据，生成的区域地下水水位分布模拟结果，对生成的预测结果采用高斯平滑函数进行平滑，进一步修正后得到最终预测结果，对采用SVR的回归算法的地下水水情预测结果和基于马尔可夫-蒙特卡洛法算的地下水水情预测结果进行模拟区域分布；
S64：基于不同的算法预测出的结果进行对比，统计分析预测的精度和计算效率，并且加入监测期间气温与降雨或/和降雪、地表沉降信息，分析不同条件对最终预测结果的影响，根据精度和计算效率分析，保存具备鲁棒性的区域地下水水情预测模型。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210500900.7A CN114970302B (zh) | 2022-05-10 | 2022-05-10 | 一种基于地下水监测***的区域地下水情预测方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210500900.7A CN114970302B (zh) | 2022-05-10 | 2022-05-10 | 一种基于地下水监测***的区域地下水情预测方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114970302A CN114970302A (zh) | 2022-08-30 |
| CN114970302B true CN114970302B (zh) | 2024-03-15 |
Family
ID=82981874
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210500900.7A Active CN114970302B (zh) | 2022-05-10 | 2022-05-10 | 一种基于地下水监测***的区域地下水情预测方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114970302B (zh) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115687883B (zh) * | 2022-11-18 | 2024-04-30 | 水利部信息中心 | 一种基于公里网格的浅层地下水蓄变量计算方法 |
| CN117390471B (zh) * | 2023-12-12 | 2024-03-26 | 山东省地质矿产勘查开发局第八地质大队(山东省第八地质矿产勘查院) | 一种水文地质数据管理方法及*** |
| CN117969789B (zh) * | 2024-03-29 | 2024-07-09 | 山东省地矿工程勘察院(山东省地质矿产勘查开发局八〇一水文地质工程地质大队) | 一种地下水资源的水情预测方法及监测*** |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2019076078A1 (zh) * | 2017-10-16 | 2019-04-25 | 中国环境科学研究院 | 地下水污染监测网多目标优化方法 |
| CN110580328A (zh) * | 2019-09-11 | 2019-12-17 | 江苏省地质工程勘察院 | 一种地下水位监测值缺失的修复方法 |
| CN114254561A (zh) * | 2021-12-17 | 2022-03-29 | 奥格科技股份有限公司 | 一种内涝预测方法、***及存储介质 |
-
2022
- 2022-05-10 CN CN202210500900.7A patent/CN114970302B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2019076078A1 (zh) * | 2017-10-16 | 2019-04-25 | 中国环境科学研究院 | 地下水污染监测网多目标优化方法 |
| CN110580328A (zh) * | 2019-09-11 | 2019-12-17 | 江苏省地质工程勘察院 | 一种地下水位监测值缺失的修复方法 |
| CN114254561A (zh) * | 2021-12-17 | 2022-03-29 | 奥格科技股份有限公司 | 一种内涝预测方法、***及存储介质 |
Non-Patent Citations (1)
| Title |
|---|
| 刘长生 ; .基于支持向量机理论的地下水水源地动态预测模型.勘察科学技术.2010,(第01期),第22-25 页. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114970302A (zh) | 2022-08-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN114970302B (zh) | 一种基于地下水监测***的区域地下水情预测方法 | |
| Hancock et al. | Spatial interpolation of large climate data sets using bivariate thin plate smoothing splines | |
| CN110263111B (zh) | 基于先验知识的土地利用/覆被信息时空监测方法 | |
| CN106780089B (zh) | 基于神经网络元胞自动机模型的永久性基本农田划定方法 | |
| CN114254561A (zh) | 一种内涝预测方法、***及存储介质 | |
| CN113610301B (zh) | 一种表征地面沉降时序演变的预测算法 | |
| Chen et al. | Local curvature entropy-based 3D terrain representation using a comprehensive Quadtree | |
| CN107133398B (zh) | 一种基于复杂网络的河流径流量预测方法 | |
| CN113051638B (zh) | 一种建筑高度优化配置方法及装置 | |
| CN113128055B (zh) | 一种基于产流系数的分布式水文模型空间率定方法 | |
| JP7153330B2 (ja) | 土砂災害予測装置、コンピュータプログラム、土砂災害予測方法及び地図情報 | |
| CN114723149A (zh) | 土壤墒情预测方法、装置、电子设备及存储介质 | |
| CN112381334B (zh) | 基于多因素模糊时间序列的高陡边坡形变趋势预测的方法 | |
| CN116738669A (zh) | 城市地表径流模拟方法、装置、设备以及存储介质 | |
| CN109063388A (zh) | 风环境模拟的微气候建筑选址设计方法 | |
| CN113987969A (zh) | 多高程尺度水流网的陆地水动态仿真模型 | |
| CN117291444B (zh) | 一种数字乡村业务管理方法及*** | |
| CN115182398B (zh) | 一种地震预警区域的地下水位及地表沉降预测方法 | |
| CN117236668A (zh) | 基于供给与消耗的区域水资源配置分析决策方法及*** | |
| CN116663395A (zh) | 基于参数优化的支持向量机回归的降雨量等值面生成方法 | |
| WO2019128018A1 (zh) | 确定风资源的方法和装置 | |
| Mikołajewski et al. | Development of cluster analysis methodology for identification of model rainfall hyetographs and its application at an urban precipitation field scale | |
| CN115438590A (zh) | 一种基于bp神经网络的降水预报订正方法及*** | |
| Ali et al. | Statistical downscaling of precipitation and temperature using gene expression programming | |
| CN117332544B (zh) | 矢量与栅格水文计算单元协同的城市雨洪模型建模方法 |
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 |