CN116680991B - 一种粘土-砂介质中dnapl污染源区识别方法及系统 - Google Patents

一种粘土-砂介质中dnapl污染源区识别方法及系统 Download PDF

Info

Publication number
CN116680991B
CN116680991B CN202310920494.4A CN202310920494A CN116680991B CN 116680991 B CN116680991 B CN 116680991B CN 202310920494 A CN202310920494 A CN 202310920494A CN 116680991 B CN116680991 B CN 116680991B
Authority
CN
China
Prior art keywords
data
cec
dnapl
field
cvae
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
Application number
CN202310920494.4A
Other languages
English (en)
Other versions
CN116680991A (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.)
Nanjing University
Original Assignee
Nanjing University
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 Nanjing University filed Critical Nanjing University
Priority to CN202310920494.4A priority Critical patent/CN116680991B/zh
Publication of CN116680991A publication Critical patent/CN116680991A/zh
Application granted granted Critical
Publication of CN116680991B publication Critical patent/CN116680991B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/38—Processing data, e.g. for analysis, for interpretation, for correction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/27—Regression, e.g. linear or logistic regression
    • 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/25—Design optimisation, verification or simulation using particle-based methods
    • 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/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/045—Combinations of networks
    • G06N3/0455—Auto-encoder networks; Encoder-decoder networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/0464—Convolutional networks [CNN, ConvNet]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/047—Probabilistic or stochastic networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/048—Activation functions
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/08—Learning methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • 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
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30—Assessment of water resources

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Artificial Intelligence (AREA)
  • Software Systems (AREA)
  • Molecular Biology (AREA)
  • Mathematical Physics (AREA)
  • Computing Systems (AREA)
  • General Health & Medical Sciences (AREA)
  • Computational Linguistics (AREA)
  • Biophysics (AREA)
  • Biomedical Technology (AREA)
  • Health & Medical Sciences (AREA)
  • Computer Hardware Design (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Geometry (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • Geophysics (AREA)
  • Remote Sensing (AREA)
  • Medical Informatics (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Probability & Statistics with Applications (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种粘土‑砂介质中DNAPL污染源区识别方法及系统，属于污染水文地质勘探技术领域。所述方法包括：第一步，基于DCIP数据的地质统计反演初步刻画含水层渗透系数K、DNAPL的饱和度S N ，第二步，利用地质统计和多相流模拟生成一系列的K、S N 场训练卷积变分自编码器CAVE网络，然后根据CVAE进行多源数据融合反演，获得K、S N 的最优估计场。本发明将DCIP数据与钻孔观测数据相结合，通过基于深度学习的CVAE联合反演策略，实现复杂粘土‑砂介质环境下DNAPL污染源区结构及其相关不确定性的精细推估。

Description

一种粘土-砂介质中DNAPL污染源区识别方法及系统
技术领域
本发明涉及污染水文地质勘探技术领域，具体涉及一种粘土-砂介质中DNAPL污染源区识别方法及系统。
背景技术
重非水相液体（Dense Non-Aqueous Phase Liquids，DNAPL）导致的地下水污染在许多工业场地中普遍存在。当DNAPL被释放到地下时，会被困在孔隙中形成不连续的离散DNAPL或在低渗透区域上方形成池状DNAPL，从而形成长期的地下水污染源。由于DNAPL空间分布受到含水层非均质性的影响，其污染源区结构非常不规则且复杂。因此，为了准确评估DNAPL的溶质通量并高效低耗地修复DNAPL污染，须首先精细刻画出DNAPL饱和度（SN）的空间分布与含水层渗透系数（K）的非均质结构。
许多方法已经被用于DNAPL污染场地中来刻画DNAPL污染源区结构。由于采样点稀疏，传统方法（例如钻孔采样、下游监测井和直接贯入技术）只能提供有限空间分辨率的DNAPL含量和岩性非均质性信息。为克服传统方法的缺陷，近年来地球物理方法被越来越多地应用到DNAPL污染场地调查中。作为非侵入式方法，地球物理方法可较廉价地提供大量连续数据。作为一种广泛应用的地球物理方法，直流电阻率（DC resistivity）方法已成功应用于DNAPL污染源区探测，溶解的NAPL污染羽刻画，以及相应水力特性的推估。
然而，对于复杂的非均质环境（例如粘土-砂介质），直流电阻率法对DNAPL源区的刻画精度是有限的。这是因为直流电阻率测量受孔隙水特性（如流体电导率）和界面电化学特性（如与岩性相关的阳离子交换容量，简称CEC）的共同影响。因此，在粘土-砂环境中，由于固体-孔隙水界面存在表面电导率，直流电阻率对于区分DNAPL和周围土壤特别是粘土介质方面效果较差。例如，直流电阻率无法区分粘土中的高饱和DNAPL和砂土中的中等饱和DNAPL。因为DNAPL倾向于在粘土地层之上蓄积，其在粘土-砂介质中的赋存很常见。在砂土和粘土混合的情况下，DNAPL形成的污染物空间分布结构在很大程度上取决于土壤的非均质性。进而又影响DNAPL污染源区的寿命。
与直流电阻率相比，时域激发极化（induced polarization，IP）方法可以通过测量含水层的电容特性（即极化率）来提供有关岩性和颗粒表面特性的额外信息，可以用来区分孔隙水及界面特性。此外，大量涌现的岩性物理模型被用来解释一系列土壤和岩石被油或溶解金属部分填充的极化机制。因此，通过结合直流电阻率和IP方法（下称DCIP）来进行DNAPL污染场地调查的案例逐渐增多。测得的极化率同样能够提供控制DNAPL污染源区结构的岩性分布信息，尤其是在粘土环境中。
尽管相比只用直流电阻率法具有一定的优势，DCIP方法可能仍然只能提供低分辨率的源区结构图像。主要原因有三个：（1）测量误差、干扰和岩性物理参数的不确定性；（2）较深区域的数据敏感度较低；（3）由于在传统地球物理反演中采用了简化的各向同性平滑，刻画的源区结构图像缺少高度不规则的复杂空间特征。
为了解决前两个问题，可以将低分辨率但空间广泛的DCIP数据与高精度但空间稀疏的钻孔数据整合在联合反演框架内；这样既可以有效地刻画出大尺度特征，又可以识别出影响所收集观测数据的更精细特征，从而最大限度地利用不同数据集中的互补信息。
第三个问题是关于传统的反演方法可能导致对DNAPL污染源区结构的推估过于平滑。为更好地体现高度不规则SN场的分布，可以应用卷积变分自编码器（ConvolutionalVariational AutoEncoder，简称CVAE）。CVAE是一种卷积神经网络，通过从训练样本中学习，可以有效地捕捉到图像中高度复杂的空间分布模式。Kang等人在“ HydrogeophysicalCharacterization of Nonstationary DNAPL Source Zones by Integrating aConvolutional Variational Autoencoder and Ensemble Smoother.Water Resources Research, 57, e2020WR028538. doi: 10.1029/2020WR028538"中展示了CVAE在参数化非均质土壤条件下的不规则DNAPL分布特性方面的效果。
虽然基于CVAE的联合反演框架已被证明能够通过整合水文地质和地球物理数据合理有效地刻画出DNAPL污染源区结构，但现有研究仅融入了单独的直流电阻率数据进行源区结构的推估，其岩性物理假设忽略了固相（例如粘土介质）的表面电导率贡献。这个假设仅对无粘土砂质介质有效。对于复杂的粘土-砂环境，单独的直流电阻率法无法区分DNAPL污染物与周围土壤对电导率的影响，限制了识别DNAPL源区结构的能力。由于极化率测量可以更好地描述岩性，以及区分DNAPL和粘土，DCIP方法有望改善SN和K场的成像精度。这也将改善利用这些地球物理数据集的联合反演性能，例如在复杂的粘土-砂介质环境中基于CVAE的源区结构反演策略。
发明内容
发明目的：为了克服现有技术在粘土-砂介质环境下无法有效区分DNAPL污染物与周围土壤的电信号的缺陷，本发明提供一种粘土-砂介质中DNAPL污染源区识别方法及系统，将DCIP数据与钻孔观测数据相结合，通过基于深度学习的CVAE联合反演策略，实现复杂粘土-砂介质环境下DNAPL污染源区结构及其相关不确定性的精细推估。
技术方案：为了实现上述发明目的，本发明提出的粘土-砂介质中DNAPL污染源区识别方法，包括以下步骤：
第一步，基于DCIP数据的地质统计反演初步刻画含水层渗透系数K、DNAPL的饱和度S N ，包括以下步骤：
（1）利用高密度电阻率调查获得的电势数据V推估含水层的体电导率σ，然后利用推估的体电导率σ图像和激发极化调查获得的视极化率数据M a 来推估标准化后的极化率M n ；
（2）利用岩性物理模型将体电导率σ和标准化后的极化率M n 转换为体积含水率θ和阳离子交换能力CEC；
（3）根据CEC与钻孔观测的含水层孔隙度φ拟合建立场地的CEC-φ关系，并根据CEC与钻孔观测的含水层渗透系数K拟合建立CEC-K关系；
（4）通过建立的CEC-φ与CEC-K关系，将CEC转化为φ和K场，进而获得S N 场；
第二步，利用地质统计和多相流模拟生成一系列的K、S N 场训练卷积变分自编码器CVAE网络，所述CVAE网络包括编码器部分和解码器部分，编码器部分用于将S N 和K的图像转化为潜变量β，解码器部分用于将潜变量β解码为S N 和K图像；然后进行多源数据融合的联合反演：将CVAE的潜变量β作为未知参数，运行水文地球物理正演模拟模型，对第一步利用DCIP直接反演得到的K、S N 场图像和钻孔观测数据进行数据同化，其中钻孔观测数据包括钻孔K、S N 和下游水相浓度数据，最后通过CVAE，将数据同化最终更新的β解码转化为K、S N 的最优估计场。
本发明还提供一种粘土-砂介质中DNAPL污染源区识别系统，包括：
粗略推估模块，基于DCIP数据的地质统计反演初步刻画含水层渗透系数K、DNAPL的饱和度S N ，所述粗略推估模块包括：
第一推估单元，利用高密度电阻率调查获得的电势数据V推估含水层的体电导率σ；
第二推估单元，利用推估的体电导率σ图像和激发极化调查获得的视极化率数据M a 来推估标准化后的极化率M n ；
第一转换单元，利用岩性物理模型将体电导率σ和标准化后的极化率M n 转换为体积含水率θ和阳离子交换能力CEC；
第一拟合单元，根据CEC与钻孔观测的含水层孔隙度φ拟合建立适合场地的CEC-φ关系；
第二拟合单元，根据CEC与钻孔观测的含水层渗透系数K拟合建立CEC-K关系；以及
第二转换单元，通过建立的CEC-φ与CEC-K关系，将CEC转化为φ和K场，进而获得S N 场；
精细推估模块，利用深度学习对粗略推估模块获得的K、S N 进行精细化推估，所述精细推估模块包括：
CVAE训练单元，利用地质统计和多相流模拟生成一系列的K、S N 场训练卷积变分自编码器CVAE网络，所述CVAE网络包括编码器部分和解码器部分，编码器部分用于将S N 和K的图像转化为潜变量β，解码器部分用于将潜变量β解码为S N 和K图像；以及
联合反演单元，将CVAE的潜变量β作为未知参数，运行水文地球物理正演模拟模型，对第一推估模块利用DCIP直接反演得到的K、S N 图像和钻孔观测数据进行数据同化，其中钻孔观测数据包括钻孔K、S N 和下游水相浓度数据，最后通过CVAE，将数据同化最终更新的β解码转化为K、S N 的最优估计场。
本发明还提供一种计算机设备，包括：一个或多个处理器；存储器；以及一个或多个程序，其中所述一个或多个程序被存储在所述存储器中，并且被配置为由所述一个或多个处理器执行，所述程序被处理器执行时实现如上所述的粘土-砂介质中DNAPL污染源区识别方法的步骤。
本发明还提供一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现如上所述的粘土-砂介质中DNAPL污染源区识别方法的步骤。
有益效果：现有技术在仅融入直流电阻率数据的情况下，忽略了固相（例如粘土介质）的表面电导率贡献，因此仅适用于不存在粘土的砂质含水层条件。对于复杂的粘土-砂介质含水层环境，该技术尚无法有效区分DNAPL污染物与周围土壤的电信号。本发明引入DCIP方法，在两步联合反演框架下，将DCIP数据与钻孔观测数据相结合，通过基于深度学习的CVAE联合反演策略，实现复杂粘土-砂介质环境下DNAPL污染源区结构及其相关不确定性的精细推估。相比直流电阻率方法，DCIP方法可以通过测量极化率来提供有关岩性和颗粒表面特性的额外信息，可以区分孔隙水及界面特性，进而区分DNAPL污染物与周围土壤。
附图说明
图1为本发明提出的两步联合反演框架流程图；
图2为本发明的CVAE-ESMDA联合反演框架流程图；
图3为本发明实施例中理想粘土-砂介质含水层的参考参数场，其中（a）-（h）展示了不同参数相对应的三维表面图像；
图4为本发明实施例中DCIP获取观测数据示意，其中（a）为DCIP的钻孔示意图和电极分布，（b）为DCIP调查的敏感度图；
图5为本发明实施例中参考和从算例Geosta-DCIP（同时使用电阻率和极化率数据）和Geosta-DC（仅使用直流电阻率数据）中推估的电导率（σ）、标准化极化率（M n ）、含水率（θ）、阳离子交换能力（CEC）场；
图6为本发明实施例中渗透系数K、孔隙度与CEC的拟合示意，其中（a）为观测lnK与DCIP得到的CEC拟合；（b）为观测的孔隙度与DCIP得到的CEC拟合；
图7为本发明实施例中CVAE-Bore、CVAE-DCIP和CVAE-Joint算例的K和S N 推估场。
具体实施方式
下面结合附图对本发明的技术方案作进一步说明。
参照图1，本发明提出了一个两步联合反演框架：第一步，将直流电阻率法和IP结合起来，对S N 和K的空间分布特征进行粗略估计，以及第二步，这些DCIP推估的S N 和K被用作约束数据，随后与钻孔观测数据相结合（即观测K、S N 和下游DNAPL水相浓度数据），在基于深度学习的CVAE联合反演策略下，提供DNAPL污染源区结构及其相关不确定性的精细推估。
第一步Step A，基于DCIP数据的地质统计反演初步刻画K、S N 。
（1）对体电导率σ和标准化的极化率M n 的地质统计反演。
首先，利用水文地球物理场地调查获得的电势数据V推估含水层的体电导率σ，然后利用推估的体电导率σ图像和场地调查获得的视极化率数据M a 来推估标准化后的极化率M n （其中M n = M·σ，M为极化率）。二者的反演分别采用PCGA（Principal ComponentGeostatistical Approach）和线性反演算法。σ和M的控制方程分别如下：
(1)
(2)
根据公式(1)和公式(2)，如已知M a 与V的观测数据，即可由地质统计学反演算法求解出σ和M的空间分布。式中，δ为狄拉克函数，r代表点电流源（其电流值为I，单位安培），表示微分算子。(2)式中的Ψ算子代表(1)式所体现的V-σ的控制方程。
（2）利用岩性物理模型将σ和M n 转换为θ和CEC。
采用以下岩性物理模型将σ和M n 转换为体积含水率θ和阳离子交换能力CEC：
(3)
(4)
其中σ w 为孔隙水电导率（S/m），本发明中视为已知；R为引入的无量纲常数，R =λ/B≈ 0.10；m为胶结指数；ρ g 表示固相的质量密度（对于砂和粘土矿物，2650 kg/m3）；λ表示极化作用下抗衡离子的视迁移率（m2V-1s-1）；B表示表面传导的抗衡离子的视迁移率（m2V-1s-1）。
（3）通过拟合钻孔数据建立CEC-φ与CEC-K关系。
由于缺乏含水层孔隙度φ的非均质性信息，即使上一步推出了体积含水率θ，也无法得知水的饱和度S w 和DNAPL的饱和度S N （S w =θ/φ，S N = 1 -S w ，适用于饱水带）。然而，由于CEC和φ都与粘土含量cl有关，因此可以通过拟合DCIP推估的CEC与钻孔观测的φ来建立特地场地的CEC-φ经验关系。此外，CEC对于推估渗透系数K也十分有用。在本发明方法中，与CEC-φ相同，同样采用拟合的方式获得CEC-K关系。渗透率k和渗透系数K可以根据公式k=K ζ/ρ进行换算，ρ和ζ分别表示孔隙水的密度和动力粘滞系数。
本发明通过最小二乘线法性拟合出CEC与钻孔观测的φ、K之间的CEC-φ、CEC-K关系。最小二乘法是常用的线性回归解法，它通过最小化误差的平方和来寻找数据的最佳函数匹配。最小二乘法的目的是找到因变量y与自变量x之间的函数关系y=f(x)，这里自变量x为CEC，因变量y为钻孔观测的φ、K。
（4）将CEC和θ转化为K和S N 场。
通过建立的CEC-φ与CEC-K关系，可以将CEC转化为φ和K场，进而可以获得S N 场（S w =θ/φ，S N = 1 -S w ，适用于饱水带）。
第二步Step B，基于深度学习的联合反演。
通过上一步的反演，初步获得了K、S N 场的分布图像。然而由于岩性物理模型的不确定性，以及地质统计先验的过度简化假设，初步反演获得的图像分辨率较低。为了进一步提高刻画精度，本发明采用基于CVAE-ESMDA（Convolutional Variational AutoEncoder -Ensemble Smoother with Multiple Data Assimilation）的联合反演框架来融合以下多源信息：（i）通过上一步DCIP反演获得的低分辨率但连续的K、S N 信息；（ii）精度高但稀疏的钻孔观测数据，包括K，S N 和下游水相浓度数据；（iii）CVAE提供的复杂DNAPL入渗模式特征。
图2展示了CVAE-ESMDA的联合反演框架。本发明中使用的神经网络为卷积变分自编码器（CVAE）。由于希望神经网络（CVAE）学习到S N 与K的空间分布特征，网络训练的输入图像与输出图像均为S N 和K的三维空间分布图。本发明中所使用的CVAE包含两部分：编码器部分（encoder）和解码器部分（decoder）。编码器部分由六个卷积层（convolutional layer）组成，其中第一个卷积层含有16个滤波器，每个滤波器可处理的S N 和K三维图片像素的大小为（3x3x3），步长值为2，激活函数为ReLU。第二个卷积层有32个滤波器，每个滤波器可处理的S N 和K三维图片像素的大小为（3x3x3），步长值为2，激活函数为ReLU。第三个卷积层有64个滤波器，每个滤波器可处理的S N 和K图片像素的大小为（3x3x3），步长值为2，激活函数为ReLU。第四个卷积层有128个滤波器，每个滤波器可处理的S N 和K图片像素的大小为（3x3x3），步长值为1，激活函数为ReLU。第五个卷积层有4个滤波器，每个滤波器可处理的S N 和K图片像素的大小为（3x3x3），步长值为1，激活函数为线性函数。第六个卷积层有4个滤波器，每个滤波器可处理的S N 和K图片像素的大小为（3x3x3），步长值为1，激活函数为线性函数。
经过编码器后，S N 和K的图像将被转化为潜变量β（latent variable），该潜变量服从标准正态分布，且该潜变量的维度（个数）远远低于原始S N 和K图像的像素点数，因而达到参数降低维度的目的。
而解码器主体部分是编码器部分的镜像。解码器由六个反卷积层构成：第一个反卷积层有128个滤波器，每个滤波器可处理的S N 和K图片像素的大小为（3x3x3），步长值为1，激活函数为ReLU。第二个反卷积层有64个滤波器，每个滤波器可处理的S N 和K图片像素的大小为（3x3x3），步长值为2，激活函数为ReLU。第三个反卷积层有32个滤波器，每个滤波器可处理的S N 和K图片像素的大小为（3x3x3），步长值为2，激活函数为ReLU。第四个反卷积层有16个滤波器，每个滤波器可处理的S N 和K图片像素的大小为（3x3x3），步长值为2，激活函数为ReLU。第五个反卷积层有1个滤波器，可处理的S N 和K图片像素的大小为（3x3x3），步长值为1，激活函数为线性函数。第六个反卷积层有10个滤波器，每个滤波器可处理的S N 和K图片像素的大小为（3x3x3），步长值为1，激活函数为softmax。
在第二步中，首先，通过仿真模拟符合迁移规律的图像来训练CVAE网络，本发明中利用地质统计和多相流模拟生成一系列的K、S N 场，训练CVAE网络。上述网络的训练过程可方便地在python环境下的Keras包中实现。神经网络中的权重参数均经由训练自动获得。经此训练过程，CVAE网络可对S N 和K场进行参数化，即学习到了针对某个具体场地的S N 和K场的空间分布特征。然后，将CVAE的中间潜变量β作为未知参数，运行水文地球物理正演模拟模型，对DCIP直接反演得到的K、S N 图像和钻孔观测数据（包括钻孔K、S N 和下游水相浓度数据）进行数据同化。最后，通过CVAE，将数据同化算法ESMDA最终更新的β解码，转化为K、S N 的最优估计场。
详细的CVAE-ESMDA反演步骤如下，以下描述中K/S N 指的是K和S N 。
（a）从均值为0，标准差为1的标准正态分布N(0, I)中采样获得中间潜变量β的初始估计样本集合；
（b）利用CVAE解码器，将中间潜变量β的估计样本集合经过解码器解码成相应的K、 S N 场的估计样本集合，以此解码的S N 场将拥有DNAPL入渗模式的特征；
（c）对于（b）中获取的K、S N 场估计样本集合中的每一个样本，获取特定钻孔位置处的K、S N 和下游水相DNAPL浓度c数据，以及K、S N 场图像。其中下游水相DNAPL浓度c数据由溶质运移模型（MT3DMS）模拟得到；
（d）对于第i次迭代，通过对比观测值（钻孔观测的K、S N 和浓度c数据，以及DCIP推估的K、S N 场图像）与上述（c）中获取的K、S N 和浓度c数据，运用数据同化算法ESMDA对潜变量β进行校正；
（e）重复步骤（b）-（d），直至满足ESMDA的最大迭代次数，得到中间潜变量β的最优估计；
（f）利用CVAE解码器将中间潜变量β的最优估计经过解码器解码为K、S N 场。
以上描述了本发明方法的一般实施过程，下面通过具体实施例来进一步演示方法并说明提出的反演框架的性能。实施例中以数值案例来演示，先假定一组具有物理现实意义的待求参数场（K、θ、S N 、φ等），这些参数之间是有相关性的，因此需要先给出符合物理规律的待求参数场，再进行方法的验证。在实施例中，在具有复杂DNAPL源区结构的三维非均质含水层中进行验算。参照图3，研究区20m×20m×5m被剖分为40×40×10=16,000个离散网格。每个网格的长度为0.5m。地下水位的深度为0.5m。地下水位以上包气带的水饱和度设置为0.4。边界由恒定水力水头条件控制，因此在x轴正方向上平均水力梯度为0.01。如图3所示，x、y、z分别表示三维坐标系下待测区域某点的三维位置。含水层底部设置为无通量边界。图3中（a）-（h）展示了对应不同参数的三维表面图像，其中K、S N 、θ、CEC、cl、σ、M、M n 分别表示渗透系数、DNAPL饱和度、体积含水率、阳离子交换能力、粘土含量、体电导率、极化率以及标准化后的极化率。
本实施例中应用了一种基于谱的方法来生成理想的渗透系数K场，相关参数列于表1中。所生成参考K场（图3中的(a)）的平均值为5.17×10-7m/s，方差为2.0（对于lnK），代表一种具有粘土-砂混合介质的强非均质含水层。基于参考K场，利用SIP（StochasticInvasion Percolation）随机逾渗模拟方法来生成S N 的分布场。得到的参考DNAPL饱和度场如图3中的(b)所示。需要说明的是，该实施例中假定污染源区仅含有单一成分污染物，即油相三氯乙烯（TCE）。
表1数值算例的参考模型参数
在本实施例中，利用无粘土砂的渗透率下限k sd （单位m2）来计算体积粘土含量（cl，无量纲）。k sd 可以通过以下经验方程得到：
(5)
其中d sd 表示砂的颗粒直径（单位m）；φ sd 表示砂的孔隙度；m是胶结系数。渗透率大于k sd 认为是不含粘土的砂介质，否则介质被认为含有一定量的粘土。
基于粘土-砂混合介质的渗透性（k，单位m 2），利用以下方程式来计算体积粘土含量cl：
(6)
其中，φ cl 表示纯粘土的孔隙度。因此，根据前文生成的参考K场，计算得到粘土含量cl场。
总孔隙度f可以从粘土含量cl计算得到：
(7)
进而可以计算得到体积含水率θ=φ·S w =φ·（1-S N ）。根据第一步中所提的岩性物理模型，可以从参考的S N 和cl场计算生成体电导率σ和极化率M场。粘土-砂混合介质的等效阳离子交换能力CEC根据以下公式计算：
(8)
其中，CEC sd 代表砂的阳离子交换能力(C/kg)，CEC cl 代表粘土的阳离子交换能力(C/kg)。相关参数列于表1，生成的参考参数场见图3。
图4为本发明实施例中DCIP获取观测数据示意，其中（a）为DCIP的钻孔示意图和电极分布，（b）为DCIP调查的敏感度图。其中在污染源区附近设置钻孔，以及下游截面设置用于收集水相浓度数据的钻孔。为了获得理想算例的观测数据，在污染源区附近设定了五个钻孔，如图4的(a)中的区域中心处五个圆柱体所示。从这5个钻孔的不同深度采集了50个K和S N 的直接测量值（z= -0.5, -1.0, …, -5.0 m），z表示深度。为了建立经验CEC-φ关系，从这些钻孔中另外收集了35个φ测量值（z= -0.5, -1.0, …, -3.5 m）。从与x轴方向正交的下游截面收集了50个溶解的DNAPL浓度水样，钻孔位置如图4的(a)中的区域边缘处五个深色圆柱体所示。
为了收集粘土-砂混合物中DNAPL污染源区的DCIP响应，在地面上设置了180个电极，如图4的(a)中的黑色矩形点所示。沿每条直线，电极间距为1m，而直线之间的偏移距离（或间距）为5m。DCIP测量采用pole-dipole配置记录，得到6,444个V和6,444个Ma的测量数据。DCIP调查的灵敏度图像如图4中的(b)所示，是根据基于累积灵敏度的方法计算的。
根据图3中生成的各参数参考场，以及DNAPL溶解和溶质运移模型（可以用comsol等模拟软件模拟），可以收集多源观测数据集（即K、S N 、φ、V、M a 和c）。将标准偏差分别为0.2ln（m/s）、0.005、0.001 g/L、0.15 V和0.003的高斯噪声添加到参考lnK、φ、c、V和M a 数据中，以表示真实的观测数据情况。
由于DNAPL饱和度的测量很困难且容易出错，因此本发明实施例中不使用精确的DNAPL饱和度值，而是用三个类别来表示：没有DNAPL区域（S N =0）、低饱和度区域（0<S N ≤0.2）和高饱和度区域（S N >0.2）。
用于参数化的CVAE网络结构如图2中的(a)所示。训练的输入数据为三维K和S N 图像，尺寸为40×40×10×2=32,000。对于CVAE的训练，将K和S N 重构问题处理为分类任务，来表示高度不规则的K和S N 场。具体来说，连续S N 值分为10个离散类（即SN = 0, 0<SN<0.06,0.06≤ SN<0.12, 0.12≤ SN<0.18, 0.18≤ SN<0.30, 0.30≤ SN<0.40, 0.40≤ SN<0.50, 0.50≤ SN<0.60, 0.60≤ SN<0.70, SN≥0.70）。通过由六个完全卷积层组成的编码器，将所得的离散S N 和K图像参数化为潜在变量β。训练集包含由SIP算法生成的40,000个K和S N 样本。另外生成4000个测试样本，测试CVAE的重构性能。反演的参考S N 和K图像（图3中的(a)-(b)）未包含在训练集中，而是来自测试集。请注意，在实际情景中，K场的纵向相关长度（λ x ,λ z ）、DNAPL释放点的位置（x 0, y 0, z 0）、释放速率（v 0）和总释放质量（M 0 ）是未知的。因此，在生成训练样本时，所有这七个变量都被视为随机变量，即m、m、/>m、/>m、/>m、/>m3/s、tons。
通过使用学习率为1×10-3、batch size为32的Adam优化程序对CVAE网络进行了训练。epochs数定为200。将潜在变量维度定为500。这些超参数（例如，epochs数、学习率和batch size）是通过试错试验来设置的。使用NVIDIA Tesla V 100 GPU进行的训练耗时4.6小时。
表2归纳了本发明实施例中设计的五个理想算例。首先算例Geosta-DC和Geosta-DCIP用来验证从直流电阻率拓展到DCIP的益处。为了证明融合多源数据集（即DCIP估计的K、S N 和钻孔数据）的优势，设计了CVAE-Bore、CVAE-DCIP和CVAE-Joint对比算例，在CVAE-ESMDA反演框架中融入了不同的数据集。
表2 数值实验五种算例的参数设置
算例 G观测数据 未知参数
Geosta-DC 6,444V 16,000σ
Geosta-DCIP 6,444V/M a + 35 φ 16,000K/S N
CVAE-DCIP 从Geosta-DCIP获取的K/S N 数据 500 z
CVAE-Bore 50K/S N /c 500 z
CVAE-Joint ①从Geosta-DCIP获取的K/S N 数据②50K/S N /c 500 z
对于Geosta-DC和Geosta-DCIP算例，在应用PCGA进行σ反演时，主成分个数设置为100。对于应用CVAE-ESMDA的算例，最大迭代次数设置为N I =4，集合大小设置为N e =500。以归一化均方根误差（NRMSE）作为反演精度评价准则：
(9)
其中，N u 为研究区中的未知参数数量；x和分别为未知参数的参考值（即真实值）和最优估计值；x max 和x min 分别为参数的最大和最小值。
图5展示了表2所示Geosta-DC和Geosta-DCIP算例中推估的lnσ和M n 场。请注意，直流电阻率法只能提供σ的估计（即图5中的第一列），而DCIP可以提供σ、M n 、θ和CEC（即第一至第四列）的估计。在算例Geosta-DC中，单独的直流电阻率法只能提供σ的推估结果（图5中的(g)）。推估的lnσ可以重现参考电导率场的主要空间模式（例如地表附近的低电导率区域和x= 7~13 m,y= 16~20 m的电阻性砂介质区域）。然而其不能揭示DNAPL区域的显著低导电率异常情况，部分可能原因是PCGA中使用的高斯先验协方差模型偏离了lnσ的真实参考空间分布。此外，因为σ受孔隙水和岩性的共同制约，仅使用σ区分不了绝缘DNAPL和背景电阻性砂介质区域。
图5显示了DCIP估计的M n 可以体现真实M n 场的大尺度特征，但高估了x= 0~10 m,y=0~8 m区域的M n 值。基于DCIP推估的σ和M n 场，可以利用岩性物理模型计算出θ和CEC的分布。DCIP得出的θ可以识别出DNAPL源区的位置（图5第三列中的黑色点缀矩形）。此外，推估的CEC场成功地划出了粘土区域（对应于高CEC值）和砂质区域（低CEC值）。
针对算例Geosta-DCIP，通过将DCIP推导的CEC和钻孔观测的φ、K拟合，得到了CEC-φ与CEC-K的特定场地经验关系。图6的(a)和(b)分别展示了CEC-K、CEC-φ关系的拟合结果，最终确定性系数（R 2）高于0.7。需要说明的是，虽然拟合的CEC-φ、CEC-K关系是基于假想数据，但基于实验室的实验研究也表明，CEC与多孔介质的渗透率、孔隙度和粘土含量高度相关。
利用建立的CEC-K、CEC-φ关系，DCIP得到的CEC、θ可以转换为K和S W （S W =θ/ φ）。相应的K、S N 推估结果如图5中（k）-（l）所示。DCIP推估的K场反映了参考K分布的大尺度空间特征（如x= 5~20 m,y= 12~20 m,z= -2~-4 m的高渗区域），但是高估了低渗透区域的程度（如x= 0~13 m,y= 0~10 m,z= -2~-4 m）。推估的S N （即1-S W ）场成功地揭示了高S N 区的位置和地表附近的非饱和层。尽管如此，推估的S N 场相对平稳，总的来说高估了DNAPL源区的扩散程度。主要原因可归纳为：（1）DCIP数据中包含的信息不完整（不足以实现高分辨率DNAPL成像）；（2）地质统计反演方法中考虑的是简化先验协方差模型；（3）CEC-K、CEC-φ关系的误差。
图7展示了CVAE-Bore、CVAE-DCIP和CVAE-Joint（即本发明）算例的K和S N 推估场。三维S N 图像中的灰色圆柱体表示钻孔的位置。表3总结了每一个算例的NRMSEs。在算例CVAE-Bore中，仅使用钻孔数据（即K、S N 和c）。推估的K揭示了非均质参考K场的主要分布模式。然而，由于钻孔采样的空间分辨率有限，仅靠钻孔数据无法重现参考污染源区结构的形状和大小。例如，除研究区中心的钻孔（x= 10 m,y= 10 m）外，在其它钻孔中均未检测到DNAPL污染物。所以，推估的污染源区结构仅局限于这五个钻孔内的区域。
表3 各算例K/S N 推估场的NRMSE以及计算时间
算例 NRMSE S NRMSE K 计算时间 (min)
Geosta-DC / / 241
Geosta-DCIP 0.25 0.53 243
CVAE-DCIP 0.19 0.41 1
CVAE-Bore 0.09 0.36 315
CVAE-Joint 0.08 0.31 326
在CVAE-DCIP算例中，K的最优估计可以刻画出参考K分布的大尺度特征，推估的S N 可以捕获y轴上高DNAPL饱和度（S N >0.4）的范围。此外，CVAE-DCIP推估的K/S N 可获得比Geosta-DCIP算例结果更高的精度（见表3）。原因主要有：（1）CVAE可以更好地表示高度不规则的K/S N 场；（2）在CVAE-DCIP中，仅使用来自DCIP高灵敏度区域的更可靠的K/S N 数据（见图4中的(b)）来避免灵敏度降低造成的误差。尽管如此，CVAE-DCIP的结果仍然高估了低渗透区域的范围，以及污染源区的总体扩散程度。也就是说，由于DCIP数据分辨率较低，无法仅靠反演该类型数据准确捕捉源区的形状和范围。
在CVAE-Joint算例中，融合了所有的三种数据信息，即DCIP估计的K/S N 、钻孔观测数据和CVAE提供的K/S N 的准确先验。结果表明，与使用单个数据集相比，CVAE-Joint更加精细地刻画出了参考K/S N 场的分布模式。在所有三种信息来源中，DCIP提供了低分辨率但在空间上连续的K/S N 信息。作为补充，少量的高精度钻孔数据为DCIP直接反演的过度扩散源区结构提供了硬约束。最后，CVAE提供的K/S N 场先验特征可以迫使污染源区结构的推估结果保有DNAPL入渗的空间分布模式。因此，利用三种互补的数据信息，可以反演得到相近于参考的污染源区模式。K/S N 推估场的NRMSE（见表3）也表明了CVAE-Joint的优势。
基于与方法实施例相同的技术构思，本发明还提供一种粘土-砂介质中DNAPL污染源区识别系统，包括：
粗略推估模块，基于DCIP数据的地质统计反演初步刻画含水层渗透系数K、DNAPL的饱和度S N ，所述粗略推估模块包括：
第一推估单元，利用高密度电阻率调查获得的电势数据V推估含水层的体电导率σ；
第二推估单元，利用推估的体电导率σ图像和激发极化调查获得的视极化率数据M a 来推估标准化后的极化率M n ；
第一转换单元，利用岩性物理模型将体电导率σ和标准化后的极化率M n 转换为体积含水率θ和阳离子交换能力CEC；
第一拟合单元，根据CEC与钻孔观测的含水层孔隙度φ拟合建立适合场地的CEC-φ关系；
第二拟合单元，根据CEC与钻孔观测的含水层渗透系数K拟合建立CEC-K关系；以及
第二转换单元，通过建立的CEC-φ与CEC-K关系，将CEC转化为φ和K场，进而获得S N 场；
精细推估模块，利用深度学习对粗略推估模块获得的K、S N 进行精细化推估，所述精细推估模块包括：
CVAE训练单元，利用地质统计和多相流模拟生成一系列的K、S N 场训练卷积变分自编码器CVAE网络，所述CVAE网络包括编码器部分和解码器部分，编码器部分用于将S N 和K的图像转化为潜变量β，解码器部分用于将潜变量β解码为S N 和K图像；以及
联合反演单元，将CVAE的潜变量β作为未知参数，运行水文地球物理正演模拟模型，对第一推估模块利用DCIP直接反演得到的K、S N 图像和钻孔观测数据进行数据同化，其中钻孔观测数据包括钻孔K、S N 和下游水相浓度数据，最后通过CVAE，将数据同化最终更新的β解码转化为K、S N 的最优估计场。
其中，第一推估单元利用如下控制方程由电势数据V推估含水层的体电导率σ，其中/>为微分算子，δ为狄拉克函数，r代表点电流源，其电流值为I，单位安培。
第二推估单元根据推估极化率M，根据M n =M·σ获得标准化后的极化率M n ，其中Ψ算子代表V-σ的控制方程。
第一转换单元利用如下岩性物理模型：和来转换θ和CEC，其中σ w 为孔隙水电导率；R为无量纲常数，R =λ/B≈ 0.10；m为胶结指数；ρ g 表示固相的质量密度；λ表示极化作用下抗衡离子的视迁移率；B表示表面传导的抗衡离子的视迁移率。
第一拟合单元和第二拟合单元分别根据钻孔观测的φ、K用最小二乘线性拟合出CEC-φ关系、CEC-K关系。
第二转换单元根据S w =θ/φ获得水的饱和度S w ，根据S N = 1 -S w 获得DNAPL的饱和度S N 。
联合反演单元进行多源数据融合的联合反演具体步骤包括：
（a）从均值为0，标准差为1的标准正态分布N(0, I)分布中采样获得潜变量β的初始估计样本集合；
（b）利用CVAE解码器，将潜变量β的估计样本集合解码成相应的K、S N 场的估计样本集合，以此解码的S N 场将拥有DNAPL入渗模式的特征；
（c）对于K、S N 场估计样本集合中的每一个样本，获取指定钻孔位置处的K、S N 和下游水相DNAPL浓度c数据，以及K、S N 场图像，其中下游水相DNAPL浓度c数据由溶质运移模型模拟得到；
（d）对于第i次迭代，将钻孔观测的K、S N 和浓度c数据，以及DCIP推估的K、S N 场图像作为观测值，通过对比观测值与上述（c）中获取的K、S N 和浓度c数据，运用数据同化算法对潜变量β进行校正；
（e）重复步骤（b）-（d），直至满足算法的最大迭代次数，得到潜变量β的最优估计；
（f）利用CVAE解码器将潜变量β的最优估计解码成K、S N 场。
其中，钻孔观测数据的获取包括：在污染源区附近设定多个钻孔，从多个钻孔的不同深度采集K和S N 的直接测量值，并收集f测量值；从与x轴方向正交的下游截面收集溶解的DNAPL浓度水样。
本发明还提供一种计算机设备，包括：一个或多个处理器；存储器；以及一个或多个程序，其中所述一个或多个程序被存储在所述存储器中，并且被配置为由所述一个或多个处理器执行，所述程序被处理器执行时实现如上所述的粘土-砂介质中DNAPL污染源区识别方法的步骤。
本发明还提供一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现如上所述的粘土-砂介质中DNAPL污染源区识别方法的步骤。
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、计算机设备或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质（包括但不限于磁盘存储器、CD-ROM、光学存储器等）上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法的流程图来描述的。应理解可由计算机程序指令实现流程图中的每一流程以及流程图中的流程的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程中指定的功能的系统。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令系统的制造品，该指令系统实现在流程图一个流程或多个流程中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程中指定的功能的步骤。

Claims (10)

1.一种粘土-砂介质中DNAPL污染源区识别方法，其特征在于，包括以下步骤：
第一步，基于DCIP数据的地质统计反演初步刻画含水层渗透系数K、DNAPL的饱和度S N ，包括以下步骤：
（1）利用高密度电阻率调查获得的电势数据V推估含水层的体电导率σ，然后利用推估的体电导率σ图像和激发极化调查获得的视极化率数据M a 来推估标准化后的极化率M n ；
（2）利用岩性物理模型将体电导率σ和标准化后的极化率M n 转换为体积含水率θ和介质阳离子交换能力CEC；
（3）根据CEC与钻孔观测的含水层孔隙度φ拟合建立场地的CEC-φ关系，并根据CEC与钻孔观测的含水层渗透系数K拟合建立CEC-K关系；
（4）通过建立的CEC-φ与CEC-K关系，将CEC转化为φ和K场，进而获得S N 场；
第二步，利用地质统计和多相流模拟生成一系列的K、S N 场训练卷积变分自编码器CVAE网络，所述CVAE网络包括编码器部分和解码器部分，编码器部分用于将S N 和K的图像转化为潜变量β，解码器部分用于将潜变量β解码为S N 和K图像；然后进行多源数据融合的联合反演：将CVAE的潜变量β作为未知参数，运行水文地球物理正演模拟模型，对第一步利用DCIP直接反演得到的K、S N 场图像和钻孔观测数据进行数据同化，其中钻孔观测数据包括钻孔K、S N 和下游水相浓度数据，最后通过CVAE，将数据同化最终更新的β解码转化为K、S N 的最优估计场。
2.根据权利要求1所述的方法，其特征在于，所述步骤（1）中，由电势数据V推估含水层的体电导率σ的控制方程如下：，/>为微分算子，δ为狄拉克函数，r代表点电流源，其电流值为I，单位安培；
由体电导率σ图像和视极化率数据M a 推估标准化后的极化率M n 包括：根据推估极化率M，根据M n = M·σ获得标准化后的极化率M n ，其中Ψ算子代表V-σ的控制方程。
3.根据权利要求1所述的方法，其特征在于，所述步骤（2）中，岩性物理模型如下：
，
，
其中σ w 为孔隙水电导率；R为无量纲常数，R =λ/B≈ 0.10；m为胶结指数；ρ g 表示固相的质量密度；λ表示极化作用下抗衡离子的视迁移率；B表示表面传导的抗衡离子的视迁移率。
4.根据权利要求1所述的方法，其特征在于，所述步骤（3）中，根据钻孔观测的φ、K，用最小二乘线性拟合出CEC-φ关系、CEC-K关系。
5.根据权利要求1所述的方法，其特征在于，所述步骤（4）中，根据S w =θ/φ获得水的饱和度S w ，根据S N = 1 - S w 获得DNAPL的饱和度S N 。
6.根据权利要求1所述的方法，其特征在于，多源数据融合的联合反演包括：
（a）从均值为0，标准差为1的标准正态分布N(0, I)分布中采样获得潜变量β的初始估计样本集合；
（b）利用CVAE解码器，将潜变量β的估计样本集合解码成相应的K、S N 场的估计样本集合，以此解码的S N 场将拥有DNAPL入渗模式的特征；
（c）对于K、S N 场估计样本集合中的每一个样本，获取指定钻孔位置处的K、S N 和下游水相DNAPL浓度c数据，以及K、S N 场图像，其中下游水相DNAPL浓度c数据由溶质运移模型模拟得到；
（d）对于第i次迭代，将钻孔观测的K、S N 和浓度c数据，以及DCIP推估的K、S N 场图像作为观测值，通过对比观测值与上述（c）中获取的K、S N 和浓度c数据，运用数据同化算法对潜变量β进行校正；
（e）重复步骤（b）-（d），直至满足算法的最大迭代次数，得到潜变量β的最优估计；
（f）利用CVAE解码器将潜变量β的最优估计解码成K、S N 场。
7.根据权利要求1所述的方法，其特征在于，钻孔观测数据的获取包括：
在污染源区附近设定多个钻孔，从多个钻孔的不同深度采集K和S N 的直接测量值，并收集φ测量值；从与x轴方向正交的下游截面收集溶解的DNAPL浓度水样。
8.一种粘土-砂介质中DNAPL污染源区识别系统，其特征在于，包括：
粗略推估模块，基于DCIP数据的地质统计反演初步刻画含水层渗透系数K、DNAPL的饱和度S N ，所述粗略推估模块包括：
第一推估单元，利用高密度电阻率调查获得的电势数据V推估含水层的体电导率σ；
第二推估单元，利用推估的体电导率σ图像和激发极化调查获得的视极化率数据M a 来推估标准化后的极化率M n ；
第一转换单元，利用岩性物理模型将体电导率σ和标准化后的极化率M n 转换为体积含水率θ和阳离子交换能力CEC；
第一拟合单元，根据CEC与钻孔观测的含水层孔隙度φ拟合建立适合场地的CEC-φ关系；
第二拟合单元，根据CEC与钻孔观测的含水层渗透系数K拟合建立CEC-K关系；以及
第二转换单元，通过建立的CEC-φ与CEC-K关系，将CEC转化为φ和K场，进而获得S N 场；
精细推估模块，利用深度学习对粗略推估模块获得的K、S N 进行精细化推估，所述精细推估模块包括：
CVAE训练单元，利用地质统计和多相流模拟生成一系列的K、S N 场训练卷积变分自编码器CVAE网络，所述CVAE网络包括编码器部分和解码器部分，编码器部分用于将S N 和K的图像转化为潜变量β，解码器部分用于将潜变量β解码为S N 和K图像；以及
联合反演单元，将CVAE的潜变量β作为未知参数，运行水文地球物理正演模拟模型，对第一推估模块利用DCIP直接反演得到的K、S N 图像和钻孔观测数据进行数据同化，其中钻孔观测数据包括钻孔K、S N 和下游水相浓度数据，最后通过CVAE，将数据同化最终更新的β解码转化为K、S N 的最优估计场。
9.一种计算机设备，其特征在于，包括：
一个或多个处理器；
存储器；以及
一个或多个程序，其中所述一个或多个程序被存储在所述存储器中，并且被配置为由所述一个或多个处理器执行，所述程序被处理器执行时实现如权利要求1-7中任一项所述的粘土-砂介质中DNAPL污染源区识别方法的步骤。
10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现如权利要求1-7中任一项所述的粘土-砂介质中DNAPL污染源区识别方法的步骤。
CN202310920494.4A 2023-07-26 2023-07-26 一种粘土-砂介质中dnapl污染源区识别方法及系统 Active CN116680991B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310920494.4A CN116680991B (zh) 2023-07-26 2023-07-26 一种粘土-砂介质中dnapl污染源区识别方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310920494.4A CN116680991B (zh) 2023-07-26 2023-07-26 一种粘土-砂介质中dnapl污染源区识别方法及系统

Publications (2)

Publication Number Publication Date
CN116680991A CN116680991A (zh) 2023-09-01
CN116680991B true CN116680991B (zh) 2023-11-17

Family

ID=87784040

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310920494.4A Active CN116680991B (zh) 2023-07-26 2023-07-26 一种粘土-砂介质中dnapl污染源区识别方法及系统

Country Status (1)

Country Link
CN (1) CN116680991B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117610319B (zh) * 2024-01-23 2024-04-30 南京大学 融合浓度和电导率的渗透系数与弥散度识别方法、系统

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109541149A (zh) * 2018-12-29 2019-03-29 广东工业大学 一种生物型大气重金属污染检测判别装置和应用
CN112149353A (zh) * 2020-09-24 2020-12-29 南京大学 基于卷积神经网络识别dnapl污染物在地下含水层分布的方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109541149A (zh) * 2018-12-29 2019-03-29 广东工业大学 一种生物型大气重金属污染检测判别装置和应用
CN112149353A (zh) * 2020-09-24 2020-12-29 南京大学 基于卷积神经网络识别dnapl污染物在地下含水层分布的方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
土工测试与勘察技术研究进展;蔡正银;周宏磊;蔡国军;郭万里;朱洵;;土木工程学报(第05期);全文 *
泥质砂岩储层有效介质电阻率模型进展与展望;韩进国;梁宇;黄布宙;杜利明;王璐;陈琦;华北;伊善强;燕军利;;矿产勘查(第01期);全文 *
砂岩型铀矿地浸浸出液溶质运移过程的综合地球物理监测;何柯;李建华;赵远程;魏文博;叶高峰;王刚;;物探与化探(第03期);全文 *

Also Published As

Publication number Publication date
CN116680991A (zh) 2023-09-01

Similar Documents

Publication Publication Date Title
CN112149353B (zh) 基于卷积神经网络识别dnapl污染物在地下含水层分布的方法
Tian-chyi et al. Uniqueness, scale, and resolution issues in groundwater model parameter identification
Eppstein et al. Simultaneous estimation of transmissivity values and zonation
Muniruzzaman et al. Experimental investigation of the impact of compound‐specific dispersion and electrostatic interactions on transient transport and solute breakthrough
Hopmans et al. How useful are small-scale soil hydraulic property measurements for large-scale vadose zone modeling?
CN116680991B (zh) 一种粘土-砂介质中dnapl污染源区识别方法及系统
Kang et al. Coupled hydrogeophysical inversion to identify non-Gaussian hydraulic conductivity field by jointly assimilating geochemical and time-lapse geophysical data
Smal et al. An automatic segmentation algorithm for retrieving sub-resolution porosity from X-ray tomography images
Kang et al. Coupled hydrogeophysical inversion of DNAPL source zone architecture and permeability field in a 3D heterogeneous sandbox by assimilation time-lapse cross-borehole electrical resistivity data via ensemble Kalman filtering
Mattern et al. Estimating travel time of recharge water through a deep vadose zone using a transfer function model
Esmaeilpour et al. Scale-dependent permeability and formation factor in porous media: Applications of percolation theory
Aghasi et al. A geometric approach to joint inversion with applications to contaminant source zone characterization
Kang et al. Characterization of DNAPL source zones in clay-sand media via joint inversion of DC resistivity, induced polarization and borehole data
Han et al. Characterization of the non-Gaussian hydraulic conductivity field via deep learning-based inversion of hydraulic-head and self-potential data
Kang et al. Integration of deep learning‐based inversion and upscaled mass‐transfer model for DNAPL mass‐discharge estimation and uncertainty assessment
CN115292890A (zh) 基于多源辅助数据开发的场地土壤污染物浓度三维空间预测方法
Sreeparvathy et al. Application of ERT, saline tracer and numerical studies to delineate preferential paths in fractured granites
CN116312824B (zh) 针对地下水污染物非菲克弥散的含水层异质性识别方法
CN117610319B (zh) 融合浓度和电导率的渗透系数与弥散度识别方法、系统
Yang et al. Electrical resistivity tomography through reinforced concrete floor
CN116819647B (zh) 一种基于交叉梯度结构约束的水文地球物理数据融合方法
Jiao et al. Functional parameterization for hydraulic conductivity inversion with uncertainty quantification
Qi Analysis of Vadose Zone Well Injection Performance
CN116165101A (zh) 一种低渗透孔隙-裂隙多尺度结构中非水相污染物的时空传输追踪方法
Capozzoli et al. Integrated geophysical and hydraulic methodologies for the study of contaminant transport process in the subsoil: A sand box experiment

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