一种非均质含水层中渗透系数随机场与孔隙度随机场协同构建方法

专利名称：一种非均质含水层中渗透系数随机场与孔隙度随机场协同构建方法

专利类型：发明专利

专利申请号：CN202310037190.3

专利申请（专利权）人：中国矿业大学
权利人地址：江苏省徐州市大学路1号

专利发明（设计）人：强静,张双圣,刘汉湖

专利摘要：本发明涉及环保和水资源技术领域，具体为一种非均质含水层中渗透系数随机场与孔隙度随机场协同构建方法，步骤为：步骤1：构建研究区的二维直角坐标体系，收集研究区含水层内勘探孔的渗透系数和孔隙度数据，分别计算渗透系数与孔隙度的均值与方差，并确定场地渗透系数和孔隙度的相关距离；步骤2：对研究区进行网格剖分，计算产生渗透系数对数正态分布随机场；步骤3：通过关联随机变量生成孔隙度正态分布随机场。本发明依据勘探数据，在构建渗透系数对数正态分布随机场的基础上构建孔隙度正态分布随机场，实现渗透系数与孔隙度在空间上的变异性及其内在关联性。

主权利要求：
1.一种非均质含水层中渗透系数随机场与孔隙度随机场协同构建方法，其特征在于，包括：步骤1：构建研究区的二维直角坐标体系，收集研究区含水层内勘探孔的渗透系数和孔隙度数据，分别计算渗透系数与孔隙度的均值与方差，并确定场地渗透系数和孔隙度的相关距离；
步骤2：对研究区进行网格剖分，计算产生渗透系数对数正态分布随机场；
步骤3：通过关联随机变量生成孔隙度正态分布随机场；
步骤2中的应根据勘探数据的变异性进行网格剖分，确定单元格的长度，并保证每个单元格内参数值的差异性不大，渗透系数对数正态分布随机场公式为：式中,K表示场地渗透系数,为一随机变量；E(K)表示场地渗透系数的均值；cv表示变异系数, D(K)表示场地渗透系数的方差；Lε1表示服从m维正态分布N(0,C)的随机变量 ,m为场地剖分的网格数 ,C＝(Cij)m×m表示其协方差矩阵 ,其元素表示网格i与网格j之间的空间相关性,Δxij与Δyij表示网格i与网格j在X方向和Y方向上的坐标之差,Lx和Ly分别表示场地渗透系数在X方向和Y方向上的相关距离；
步骤3中的关联随机变量的推导过程如下：
对公式(1)两侧取对数得：
根据对数正态分布性质可知,lnK服从正态分布,且lnK的均值 方差D2
(lnK)＝ln(cv+1),从而上式可变形如下：
由于孔隙度服从正态分布,参考(3)式的正态分布表达形式,一种孔隙度随机场计算公式如下：其中,n表示场地孔隙度,为一随机变量；E(n)表示场地孔隙度的均值；D(n)表示场地孔隙度的方差；Lε2仍然表示服从m维正态分布N(0,C)的随机变量；
由于同一研究场地渗透系数与孔隙度密切相关,可知渗透系数随机场与孔隙度随机场并不是独立的,即Lε1与Lε2不是相互独立的，由公式(2)可得服从正态分布N(0,C)，此时引入随机变量η和常数ρ(ρ∈[‑
1,1]),假设η也服从正态分布N(0,C),且Lε1与η相互独立,则 也服从正态分布N(0,C)，可令(4)式中 从而得到与渗透系数随机场相关联的孔隙度随机场计算公式：
公式(5)所给孔隙度随机场既与渗透系数随机场有关,又包含了孔隙度自身的随机性；
用ni表示网格i的孔隙度,Ki表示网格i的渗透系数(i＝1,2,L,m),根据公式(3)以及(5)计算ni与lnKi的协方差 得到 可知ρ表示孔隙度n与渗透系数对数lnK的相关系数,且ρ∈[‑1,1]；0＜ρ≤1表示正相关，‑1≤ρ＜0表示负相关，ρ＝0表示不相关，由于渗透系数与孔隙度通过空间随机变量Lε1产生关联,故把Lε1称之为关联随机变量；
规定n取值范围为[E(n)‑δ,E(n)+δ]，其中0＜δ＜E(n)，在孔隙度随机抽样时,根据正态分布的性质，限制标准差 规定,当n＜E(n)‑δ时,孔隙度赋值为E(n)‑δ；当n＞E(n)+δ时,孔隙度赋值为E(n)+δ。
2.根据权利要求1所述的一种非均质含水层中渗透系数随机场与孔隙度随机场协同构建方法，其特征在于：步骤1中的勘探孔应根据研究区场地的面积、形状、用途、污染状况的特点，采用对角线布点法、梅花形布点法、棋盘式布点法或“S”形布点法进行勘探孔的布设，相关距离可采用递推空间法、相关函数法、变异函数法计算获得。 说明书 : 一种非均质含水层中渗透系数随机场与孔隙度随机场协同构
建方法技术领域[0001] 本发明涉及环保和水资源技术领域，具体为一种非均质含水层中渗透系数随机场与孔隙度随机场协同构建方法。背景技术[0002] 随着现代应用数学以及计算机技术的发展和广泛应用，数值模拟逐渐成为研究地下水运动规律、定量评价地下水资源以及模拟一些水文地质过程发生发展的主要手段。在较大研究区域内，水文地质参数往往具有空间变异性，表现出非均质性，因此往往需要构建水文地质参数空间分布场。另一方面构建某一参数的单一空间分布场往往存在一定的风险性，不能充分体现参数不确定性对数值模拟结果的影响，为此需要构建水文地质参数的空间分布随机场。[0003] 针对非均质地层，传统的水文地质概念模型构建方法是根据岩性的相似性，人为地将研究区划分为若干个区域，然后对每个分区赋予特定的水文地质参数。但是通过该方法构建的概念模型的精确程度跟研究者对研究区的认知水平及水文地质勘探情况等因素密切相关，受人为影响较大，可能存在“模拟失效”的情况。[0004] 目前国际上对于渗透参数与孔隙度的空间变异性做了大量的研究工作，普遍认为含水层渗透系数与孔隙度是影响地下水水流与溶质运移的主要因素，而且渗透系数服从对数正态分布，孔隙度服从正态分布。对于特定岩性的土壤，渗透系数与孔隙度可能存在着正相关、负相关及不相关三种关系。因此如何在保证两个参数空间分布的前提下构建渗透系数与孔隙度的协同随机场，对研究岩性非均质性对地下水水流及溶质运移的影响具有重要意义。发明内容[0005] 本发明的目的在于提供一种非均质含水层中渗透系数随机场与孔隙度随机场协同构建方法，通过引入关联随机变量，将渗透系数随机场与孔隙度随机场进行空间关联，并实现正向关联、负向关联及不相关的人为调控。[0006] 为实现上述目的，本发明提供如下技术方案：一种非均质含水层中渗透系数随机场与孔隙度随机场协同构建方法，包括[0007] 步骤1：构建研究区的二维直角坐标体系，收集研究区含水层内勘探孔的渗透系数和孔隙度数据，分别计算渗透系数与孔隙度的均值与方差，并确定场地渗透系数和孔隙度的相关距离；[0008] 步骤2：对研究区进行网格剖分，计算产生渗透系数对数正态分布随机场；[0009] 步骤3：通过关联随机变量生成孔隙度正态分布随机场。[0010] 优选的，步骤1中的勘探孔应根据研究区场地的面积、形状、用途、污染状况的特点，采用对角线布点法、梅花形布点法、棋盘式布点法或“S”形布点法进行勘探孔的布设，相关距离可采用递推空间法、相关函数法、变异函数法计算获得。[0011] 优选的，步骤2中应根据勘探数据的变异性进行网格剖分，确定单元格的长度，并保证每个单元格内参数值的差异性不大，渗透系数对数正态分布随机场公式为：[0012][0013] 式中,K表示场地渗透系数,为一随机变量；E(K)表示场地渗透系数的均值；cv表示变异系数, D(K)表示场地渗透系数的方差；Lε1表示服从m维正态分布N(0,C)的随机变量,m为场地剖分的网格数,C＝(Cij)m×m表示其协方差矩阵,其元素表示网格i与网格j之间的空间相关性,Δxij与Δyij表示网格i与网格j在X方向和Y方向上的坐标之差,Lx和Ly分别表示场地渗透系数在X方向和Y方向上的相关距离。[0014] 优选的，步骤3中的关联随机变量的推导过程如下：[0015] 对公式(1)两侧取对数得：[0016][0017] 根据对数正态分布性质可知,lnK服从正态分布,且lnK的均值2方差D(lnK)＝ln(cv+1),从而上式可变形如下：[0018][0019] 由于孔隙度服从正态分布,参考(3)式的正态分布表达形式,本文提出一种孔隙度随机场计算公式如下：[0020][0021] 其中,n表示场地孔隙度,为一随机变量；E(n)表示场地孔隙度的均值；D(n)表示场地孔隙度的方差；Lε2仍然表示服从m维正态分布N(0,C)的随机变量；[0022] 由于同一研究场地渗透系数与孔隙度密切相关,可知渗透系数随机场与孔隙度随机场并 不是独 立的 ,即L ε 1与L ε 2 不是相 互独立的 .由 公式 (2) 可 得服从正态分布N(0,C).此时引入随机变量η和常数ρ(ρ∈[‑1,1]),假设η也服从正态分布N(0,C),且Lε1与η相互独立,则 也服从正态分布N(0,C)，可令(4)式中 从而得到与渗透系数随机场相关联的孔隙度随机场计算公式：[0023][0024] 公式(5)所给孔隙度随机场既与渗透系数随机场有关,又包含了孔隙度自身的随机性；[0025] 用ni表示网格i的孔隙度,Ki表示网格i的渗透系数(i＝1,2,L,m),根据公式(3)以及(5)计算ni与lnKi的协方差 得到 可知ρ表示孔隙度n与渗透系数对数lnK的相关系数,且ρ∈[‑1,1]；0＜ρ≤1表示正相关，‑1≤ρ＜0表示负相关，ρ＝0表示不相关。由于渗透系数与孔隙度通过空间随机变量Lε1产生关联,故把Lε1称之为关联随机变量；[0026] 规定n取值范围为[E(n)‑δ,E(n)+δ]，其中0＜δ＜E(n)，在孔隙度随机抽样时,为了保证绝大多数样本点落在取值范围内,根据正态分布的性质，可限制标准差 虽然限制 但仍然会有少量的点落在[E(n)‑δ,E(n)+δ]外面,因此规定,当n＜E(n)‑δ时,孔隙度赋值为E(n)‑δ；当n＞E(n)+δ时,孔隙度赋值为E(n)+δ。[0027] 与现有技术相比，本发明的有益效果是：[0028] 本发明提供一种非均质含水层中渗透系数随机场与孔隙度随机场协同构建方法，通过引入关联随机变量，提出的关联随机变量能够将渗透系数与孔隙度在空间上进行有效关联，并根据研究需要，分别构建正相关、负相关及不相关的空间关系。附图说明[0029] 图1为本发明的算例模型示意图。[0030] 图2为本发明的呈正相关关系的渗透系数随机场和孔隙度随机场示意图。[0031] 图3为本发明的呈负相关关系的渗透系数随机场和孔隙度随机场示意图。[0032] 图4为本发明的呈不相关关系的渗透系数随机场和孔隙度随机场示意图。[0033] 图5为本发明的渗透系数与孔隙度散点正相关示意图。[0034] 图6为本发明的渗透系数与孔隙度散点负相关示意图。[0035] 图7为本发明的渗透系数与孔隙度散点不相关示意图。具体实施方式[0036] 下文结合附图和具体实施例对本发明的技术方案做进一步说明。[0037] 本发明提出的一种非均质含水层中渗透系数随机场与孔隙度随机场协同构建方法，包括[0038] 步骤1：构建研究区的二维直角坐标体系，收集研究区含水层内勘探孔的渗透系数和孔隙度数据，分别计算渗透系数与孔隙度的均值与方差，并确定场地渗透系数和孔隙度的相关距离；[0039] 步骤2：对研究区进行网格剖分，计算产生渗透系数对数正态分布随机场；[0040] 步骤3：通过关联随机变量生成孔隙度正态分布随机场。[0041] 本发明中步骤1中的勘探孔应根据研究区场地的面积、形状、用途、污染状况的特点，采用对角线布点法、梅花形布点法、棋盘式布点法或“S”形布点法进行勘探孔的布设，相关距离可采用递推空间法、相关函数法、变异函数法计算获得。[0042] 步骤2中应根据勘探数据的变异性进行网格剖分，确定单元格的长度，并保证每个单元格内参数值的差异性不大。渗透系数对数正态分布随机场公式为：[0043][0044] 式中,K表示场地渗透系数,为一随机变量；E(K)表示场地渗透系数的均值；cv表示变异系数, D(K)表示场地渗透系数的方差；Lε1表示服从m维正态分布N(0,C)的随机变量,m为场地剖分的网格数,C＝(Cij)m×m表示其协方差矩阵,其元素表示网格i与网格j之间的空间相关性,Δxij与Δyij表示网格i与网格j在X方向和Y方向上的坐标之差,Lx和Ly分别表示场地渗透系数在X方向和Y方向上的相关距离。[0045] 步骤3中的关联随机变量的推导过程如下：[0046] 对公式(1)两侧取对数得：[0047][0048] 根据对数正态分布性质可知,lnK服从正态分布,且lnK的均值 方2差D(lnK)＝ln(cv+1),从而上式可变形如下：[0049][0050] 由于孔隙度服从正态分布,参考(3)式的正态分布表达形式,一种孔隙度随机场计算公式如下：[0051][0052] 其中,n表示场地孔隙度,为一随机变量；E(n)表示场地孔隙度的均值；D(n)表示场地孔隙度的方差；Lε2仍然表示服从m维正态分布N(0,C)的随机变量。[0053] 由于同一研究场地渗透系数与孔隙度密切相关,可知渗透系数随机场与孔隙度随机场并 不是独 立的 ,即L ε 1与L ε 2 不是相 互独立的 ，由 公式 (2) 可 得服从正态分布N(0,C)，此时引入随机变量η和常数ρ(ρ∈[‑1,1]),假设η也服从正态分布N(0,C),且Lε1与η相互独立,则 也服从正态分布N(0,C).可令(4)式中 从而得到与渗透系数随机场相关联的孔隙度随机场计算公式：[0054][0055] 公式(5)所给孔隙度随机场既与渗透系数随机场有关,又包含了孔隙度自身的随机性。[0056] 用ni表示网格i的孔隙度,Ki表示网格i的渗透系数(i＝1,2,L,m),根据公式(3)以及(5)计算ni与lnKi的协方差 得到 可知ρ表示孔隙度n与渗透系数对数lnK的相关系数,且ρ∈[‑1,1]；0＜ρ≤1表示正相关，‑1≤ρ＜0表示负相关，ρ＝0表示不相关。由于渗透系数与孔隙度通过空间随机变量Lε1产生关联,故把Lε1称之为关联随机变量.[0057] 以孔隙度为例构建正态分布随机场.假设场地孔隙度n为一随机变量,服从正态分布,均值为E(n),方差为D(n),则理论上,孔隙度n的取值范围是(‑∞,+∞),这显然不符合实际情况。[0058] 如果n没有给定范围，为规避孔隙度n随机取值导致n≤0或者n≥1,规定n≤0时用一个很小的正数代替(比如0.001)，n≥1时用一个接近且小于1的正数代替(比如0.999)。[0059] 如果规定n取值范围为[E(n)‑δ,E(n)+δ]，其中0＜δ＜E(n).在孔隙度随机抽样时,为了保证绝大多数样本点落在取值范围内,根据正态分布的性质，可限制标准差虽然限制 但仍然会有少量的点落在[E(n)‑δ,E(n)+δ]外面,因此规定,当n＜E(n)‑δ时,孔隙度赋值为E(n)‑δ；当n＞E(n)+δ时,孔隙度赋值为E(n)+δ。[0060] 具体的[0061] 假定研究区域为矩形区域,长510m,宽300m,目标含水层为厚度20m的砂质承压含水层,西部边界与东部边界为给定水头边界,其中东部水头为25m,西部水头为30m,北部边界与南部边界为隔水边界.天然状态下地下水流为自西向东流动的二维均质各向异性的承压稳定流。[0062] 以研究区西南角为坐标原点建立坐标系。[0063] 表1研究区文地质参数的概率分布及取值特征[0064][0065] 注:NA表示非必需参数[0066] 由于研究区分成A,B,C,D四个区,渗透系数均值分别为E(KA)＝1m/d,E(KB)＝23m/d,E(KC)＝15m/d,E(KD)＝10m/d；孔隙度均值分别为E(nA)＝0.27,E(nB)＝0.32,E(nC)＝0.22,E(nD)＝0.18.设定 LxA＝LyA＝LxB＝LyB＝LxC＝LyC＝LxD＝LyD＝20m，4个区孔隙度方差设定为 场地剖分网格数n＝17000,单元格长度为3m.分别取ρ＝0.8、ρ＝‑0.8及ρ＝0表示渗透系数与孔隙度正相关、负相关及不相关。根据公式(1)和(5),采用MonteCarlo方法产生400组符合给定分布特征的渗透系数随机场和孔隙度随机场，分别列举其中1组渗透系数随机场和孔隙度随机场。[0067] 分别对构建的渗透系数随机场和孔隙度随机场进行统计,研究区内的A、B、C、D四个区域的渗透系数与孔隙度之间均呈现正相关、负相关及不相关的关系，表明提出的基于关联随机变量的渗透系数对数正态分布随机场与孔隙度正态分布随机场协同构建方法是可行的，也是可靠的。[0068] 上述具体实施例仅仅是本发明的几种优选的实施例，基于本发明的技术方案和上述实施例的相关启示，本领域技术人员可以对上述具体实施例做出多种替代性的改进和组合。

专利地区：江苏

专利申请日期：2023-01-10

专利公开日期：2024-11-29

专利公告号：CN116029136B