基于NSMC方法的某核工程...隙地下水数值模拟不确定分析_赵敬波.pdf

第 6 9 卷第 2 期2 0 2 3 年 3 月地质论评GEOLOGICAL EVIEWVol 69No 2March，2 0 2 3注:本文为核设施退役及放射性废物治理专项资助项目(编号:科工二司 2019 1496 号)的成果。收稿日期:2022-05-09;改回日期:2022-09-30;网络首发:2022-10-20;责任编辑:刘志强。Doi:1016509/jgeoreview202210081作者简介:赵敬波，男，1988 年生，博士，高级工程师，主要从事水文地质与高放废物地质处置相关研究;Email:zhaojingbobriug outlookcom。通讯作者:周志超，男，1983 年生，正高级工程师，主要高放废物地质处置相关研究;Email:zhouzhichao2006 163com。基于 NSMC 方法的某核工程场地裂隙地下水数值模拟不确定分析赵敬波1)，周志超1)，潘跃龙2)，叶浩3)，吴群2)，郭永海1)，李杰彪1)，付馨雨1)1)核工业北京地质研究院，国家原子能机构高放废物地质处置创新中心，北京，100029;2)深圳中广核工程设计有限公司，广东深圳，518000;3)长江勘测规划设计研究有限责任公司，武汉，430010内容提要:裂隙介质渗透结构表现为高度的非均质性与各项异性。为了科学有效地预测某核工程场地裂隙地下水的流动规律，揭示裂隙岩体地下水的渗流特性，笔者等采用 Pilot Point 调参方法与 null space Monte Carlo 方法(NSMC)，开展了裂隙岩体渗透结构的不确定性分析研究，构建了符合实际水文地质条件的多个渗流数值模型集合。结果表明:该方法获得的各个实现地下水位模拟结果能够与实际观测数据较好吻合，可反映工程场地裂隙地下水动力特征与流动趋势;各个实现的参数化渗透结构在空间上存在一定的差异性，但整体变化趋势是保持一致的，渗透参数的不确定性表现为在实测数据分布区域相对较低，钻孔空白区域相对较高;该方法可以弥补单一、确定性模拟结果在表征裂隙介质渗透结构方面的局限性，有效地降低模型参数的不确定性与随机性。此方法对进一步提升裂隙岩体渗流模拟精度与预测能力，深化裂隙地下水迁移规律的认识具有重要的意义。关键词:裂隙岩体;渗透结构;null space Monte Carlo 方法;不确定分析;Pilot Point 方法场地岩体的渗透性对核设施的安全性能评价极为重要。场地岩体由于裂隙的存在，其渗流特性表现极为显著的非均质、各向异性特征(张人权等，2011)。考虑到野外实测数据的有限性，实际裂隙介质渗透结构的复杂性，模型结构的高度概化，现阶段广泛采用的确定性、唯一的模拟结果存在诸多的不确定性(吴吉春和陆乐，2011;曾献奎等，2012)。因此，研究地下水模拟的不确定性，提升地下水数值模拟精度，科学合理地预测裂隙地下水流动规律，这对裂隙岩体地下工程最终安全性能评价具有重要的应用价值与现实意义。地下水模拟的不确定性按照来源可分为参数的不确定性、模型结构的不确定性和观测数据的不确定性 3 类。其中，以参数不确定性研究相对较为广泛(曾献奎，2012;宋凯等，2018)。目前关于地下水参数不确定性常见的分析方法主要包括通用似然不确定性估计(GLUE)、马尔科夫链蒙特卡洛方法(MCMC)、经典贝叶斯方法、矩方程法等(Beven andBinley，1992;Woodbury et al，2000;夏强，2011)。这类方法原理简单、实用性广、可靠性强，但计算工作量大，对于复杂的高维参数问题与高度非线性数值问题(例如:地下热水运移、CO2地质储存、污染物迁移与降解等)适应能力较差(Pappenberger andBeven，2006;Yoon et al，2013;莫绍星，2019)。与此同时，一些学者提出了地质统计相关方法，包括变异函数方法、转移概率方法等(Deutsch and Journel，1995;Carle and Fogg，1996)。这类方法原理是通过实测数据数理统计规律，获得成千上万次正演随机模型，统计所有模拟结果的集合，分析地下水模拟不确定性(ubin et al，1990;杨金忠，2000)。目前该方法在含水层非均质性方面得到了广泛应用，例如:美国亚利桑那 Avra 第四系含水层、芬兰与加拿大放射性废物处置及其他相关水力水电工程的裂隙岩体等(Clifton et al，1982;Park et al，2004;Blessent et al，2011;Chen Yifeng et al，2018)。对于裂隙介质而言，这类方法可实现渗透结构的非均质性、各向异性数值表征，一定程度上能够合理地模拟与预测裂隙地下水流动规律。然而该方法对模型不确定性分析工作所需的计算量同样相对较大，且对实测数据的空间分布及数量要求相对较高。为了模拟含水层非均质各向特征，提高模型不确定分析的数值计算效率，部分学者提出了 PilotPoint 调 参 方 法 与 null space Monte Carlo 方 法(NSMC)相结合的方式进行地下水流模拟不确定分析研究(Marsily et al，1984;Christensen et al，2008;Tonkin et al，2009)。其原理主要包括两部分内容，首先是在变异函数的基础上，采用 PilotPoint 调参方法与空间插值方法进行模型参数识别与校准，生成符合实际水文地质条件的高度参数化渗透结构模型，用于描述含水介质的非均质各向异性特征，目前已成功应用到不同尺度的地下水流数值模拟研究中(Carniato et al，2015;Muhammad etal，2018)。其次，在上述模拟结果的基础上，基于参数敏感性分析与模型反演的方式，生成一系列的随机实现，调整 Pilot Point 参数值，获得符合实际观测数据地下水模型集合，进行模型不确定性的综合分析与判定(Tonkin et al，2009;Yoon et al，2013;Herckenrath et al，2011)。考虑到该方法在含水层非均质性刻画、参数识别、数值计算效率等方面具有良好的适用性，后续学者研发了集模型参数识别、敏感性分析及不确定分析等功能于一体的数值并行程序，例如:PEST、ITOUGH、UCODE 等，其中以 PEST程序应用最为广泛(Poeter et al，1999;Finsterle etal，2012;Doherty，2015)。同时基于高性能数值计算平台或者云计算技术，采用该程序可实现上千数量模型参数的数值反演研究，使得该方法在含水层非均质各向异性刻画、污染物扩散模拟、复杂地表与地下水模型耦合模拟等方面得到一定地推广应用(Hayley et al，2014;Alberti et al，2018;Moeck etal，2019;Moges et al，2020)。然而此方法目前研究重点更多是侧重于沉积相含水层，未考虑模型参数渗透结构在垂向上变异性特征，且关于裂隙地下水模拟不确定分析的工程案例研究甚少。此外，还存在水力层析、自校正法、集合卡尔曼滤波法等其他数值方法，多数以理论探索为主，实际应用案例研究相对较少(Wen Xianhuan et al，1999;Yeh et al，2000;周海燕，2012)。笔者等以我国滨海地区某核设施花岗岩工程场地为研究对象，通过引入 NSMC 方法，开展裂隙地下水模拟不确定性分析研究，探讨此方法在裂隙介质渗流模拟方面的适用性;基于 PEST 并行程序与Pilot Point 调参方法，获取符合实际水文地质条件多个参数化渗透结构模型，综合评定研究区裂隙地下水流动特征，最终为本工程场地性能评价提供参考依据。同时为水利水电工程、放射性废物处置等其他工程提供技术支撑。1研究区概况11自然地理研究区处于我国南部沿海地区，年平均气温约为 223，多年平均降水量约为 2348 mm，降雨呈季节性波动。研究区所在山体呈北东东展布，标高分布在 17002655 m，南侧为海域，北侧为沟谷，常年存在地表径流的溪水，径流量随季节性波动较为明显。工程场地为典型基岩海蚀地貌(图 1)。图 1 研究区地形图Fig1 Topographic map of the study area12地质与水文地质条件研究区主要岩性以花岗岩为主，沟谷低洼处沉积了第四系松散堆积物和海相沉积物。根据花岗岩矿物成分含量及颗粒大小不同，可划分为粗粒二长花岗岩、斑状花岗岩、中细粒花岗岩及花岗斑岩。第四系沉积物可分为冲洪积砂土与粉质黏土、海相砂土，如图 2 所示。根据目前钻孔揭露的岩性信息可知，基岩垂向上可分为全风化层、强风化层、中风化层及微风化层。全风化厚度一般小于 118 m，强风化层厚度一般分布在 12298 m，中风化层厚度则为 13593 m，目前尚未有钻孔揭穿微风化层。第四系松散沉积物分为人工堆积层与冲洪积层的粉质黏土和砂土、海相沉积的砂土，最大厚度约为 5060 m。场地范围内无较大规模的断层穿过。根据地下水赋存条件，研究区地下水分为第四系松散孔隙水与基岩风化裂隙水。工程场地地下水位受地形影响较为明显，其中地形位置较高处钻孔CCZK01 揭露的地下水位为 2120 m，而地形较低处钻孔 KJZK01 揭露的水位则为 146 m，不同钻孔揭626地质论评2023 年露的地下水位相差较大。地下水的补给来源为大气降水，主要的排泄区域为场地山体北侧沟谷与南侧海域，少数地下水以泉的形式排泄，且研究区无人工开采的情景，地下水系统整体处于天然稳定状态。总体而言，工程场地具有相对较为独立完整的地下水系统。图 2 研究区地质简图与剖面图Fig2 Simplified geological map and profile of the study area2NSMC 方法原理考虑到目标函数的多解性、参数间的补偿作用及随机性，地下水流模拟结果往往存在“异参同效”的现象，模型参数的不确定性分析成为地下水数值模拟研究的重要组成部分。基于此，本研究整体思路是在前期工作的基础上，采用程序 PEST，通过Pilot Point 调参方法，获得单个确定的、合理的地下水参数化渗透结构模型。据此，根据渗透参数的数理统计规律，采用 NSMC 方法生成一系列随机模型，每个模型都采用 Pilot Point 方式进行识别与验证，选取符合实际水文地质条件多个参数化渗透结构模型集合，以此综合评估工程场地地下水动力场特征。首选定义观测数据与参数的满足线性关系采用式(1)表示，目标函数 为计算水位与观测水位之间差值的平方和采用式(2)表示结果如下:h=Xp+(1)=(Xp h)TQm(Xp h)(2)式中 p 为模型参数长度为 n 的列向量;h 为观测数据长度为 m 的列向量;X 为雅克比矩阵(mn)，为水位对参数的导数;矩阵 Qm为权重系数的对角线矩阵;为观测数据 h 的测量误差值。726第 2 期赵敬波等:基于 NSMC 方法的某核工程场地裂隙地下水数值模拟不确定分析p=(XTQmX)1XTQmh(3)NSMC 不确定分析方法要求式(3)采用截断奇异值分解(SVDassist 方法)进行求解，以此获得唯一的、确定性结果。其求解思路是为了确保数值计算的稳定性和高效性，划分为对观测水位较为敏感的参数集合和不敏感参数集合两部分进行矩阵的降维求解(Tonkin et al，2005)，针对矩阵 XTQmX 进行奇异值分解可表示为:XTQmX=V1V2E100E2V1V2T(4)式中 E1为对角矩阵，对角线的值为矩阵 XTQmX 的较大非零奇异值;E2同样也为对角矩阵，其包括较小的奇异值和零值，其 V1和 V2为与矩阵 E1和 E2维数相同的幺正矩阵。为了确保数值计算的稳定性和高效性，采用 SVDassist 方法处理时，考虑到矩阵E2对角线的奇异值接近为 0，这些值都近似设置为零，通过目标函数的最小化，联合式(3)和式(4)即可求解优化待估参数向量 p:p=V1E11VT1XTQmh=Gh(5)式(5)是用于约