技术领域
[0001] 本
发明涉及河口
淡水资源利用领域，具体涉及一种河口沿程淡水保证率的反演方法。
背景技术
[0002] 随着社会经济的发展，水资源对于各行各业而言作为一种越来越重要的资源。河流入海口地区常常受盐水入侵影响，
水体的
盐度大小影响着河流与湿地
生态系统，也影响了人类的水资源开发利用。很多大城市和人口聚居地均位于河口地区，淡水资源十分重要，由于盐水入侵的影响，河流入海口的取水规模选择和取水点选址对饮水安全与社会经济发展具有重要影响。
[0003] 河口盐度处于径流、潮位等水文因素和地形因素的共同作用下，并且常常受到人类取水影响，河口盐度
时空分布规律复杂，在全球海平面上升和
气候变化的大背景下，加上人类活动增多，河口地区径流、潮位等水文因素和地形因素的变化对盐水入侵的影响不可忽视。
[0004] 但相比于径流和潮位实测序列，对河口全河段盐度进行现场实测十分困难，河口长序列盐度实测资料往往十分欠缺，不便进行河口盐水入侵
风险的评估，更无法进行取水口的合理确定。因此需要一种河口沿程淡水保证率反演方法，为河口取水口的确定提供技术支持，保障河口地区的饮水安全。
发明内容
[0005] 针对
现有技术上存在的不足，本发明的目的在于提出一种基于河口径流、潮位序列和地形资料的沿程淡水保证率反演方法，通过数值模拟得到各河段不发生盐水入侵的临界径流-潮位控制线，引入Copula函数研究河口地区的径流-潮位关联特性，并考虑人类取水导致的径流减小影响，基于径流-潮位Copula联合分布函数研究了河口地区的盐水入侵及各河段淡水保证率，得出了不同取水规模下的沿程淡水保证率，并得出了各设计取水保证需求下的临界取水断面。
[0006] 为了实现上述目的，本发明是通过如下的技术方案来实现：
[0007] 一种河口沿程淡水保证率的反演方法，包括以下步骤：
[0008] A.根据河口历史淡水径流和潮位资料，对河口逐日水文和潮汐特征进行解析，基于常用且典型的分布函数(例如：线性、均匀、正态、指数、PIII，Gompertz分布函数等)对径流和潮位累计
频率曲线进行拟合，得出径流、潮位的边缘分布特性，并进一步考虑人为取水的影响，得到人为取水影响后的径流分布函数，取水后的径流分布函数即为天然径流方程发生平移，平移的单元数与取水规模相对应；
[0009] B.根据步骤A所得的河口径流、潮位的边缘分布特性，考虑两者之间的相互影响和关联特性，建立径流-潮位二维联合分布及其联合风险概率模型；
[0010] C.基于河口实测地形，构建三维水动
力和盐度数学模型，所述三维水动力和盐度数学模型上游边界为淡水径流，下游边界为外海潮位；
[0011] D.基于步骤C所得的三维水动力和盐度数学模型，根据河口径流和潮位分布范围，以径流、潮位的逐项组合来设置计算工况，进行各工况下河口盐水入侵三维数值模拟，以径流、潮位的逐项组合来进行盐水入侵模拟，可得径流、潮位逐项组合下的最大盐水入侵距离(以盐度为0.45‰作为盐水入侵
阈值)，反推得不同盐水入侵距离下的临界径流-潮位控制线，即河口各河段沿程发生盐水入侵的临界径流-潮位控制线，并得各河段不发生盐水入侵的径流和潮位区间，即淡水保证区；
[0012] E.基于步骤B得出的径流-潮位二维联合分布和步骤D得出的河口各河段沿程临界径流-潮位控制线，计算出盐水入侵的河口沿程淡水保证率；
[0013] F.根据步骤E得出的河口沿程淡水保证率，对比河口取水的设计供水保证率，得出不同取水规模下的临界取水断面，即取水口的临界布置
位置。
[0014] 上述一种河口沿程淡水保证率的反演方法，步骤B中所述的径流-潮位二维联合分布采用Copula函数，并包括以下步骤：
[0015] B1)令H为具有边际分布F和G的联合分布函数，那么存在一个Copula函数C，满足：
[0016] H(x1,x2)＝C(F(x1),G(x2)) (1)
[0017] 其中，x1,x2为边缘随机变量，F为变量x1的边缘分布函数，G为变量x2的边缘分布函数。
[0018] 常用的二元Copula函数总体上可以分为三类，二元正态Copula函数、二元t-Copula函数和阿基米德Copula函数；其中含有一个参数的阿基米德Copula函数应用最广，包括Gumbel Copula函数、Clayton Copula函数和Frank Copula函数。因此，步骤B1中的径流-潮位二维联合分布采用Copula函数，选用Gumbel Copula函数、Clayton Copula函数、Frank Copula函数、二元正态Copula函数以及二元t-Copula函数中的任意一种或者任意几种组合使用，几种二元Copula函数如下所示：
[0019] (2-1)Gumbel Copula函数
[0020] C(x1,x2)＝exp{-[(-lnx1)θ+(-lnx2)θ]1/θ}，θ∈[1,∞) (2)
[0021] 式中，θ为相关参数，x1,x2为边缘随机变量，下同；由式(2)可知，Gumbel Copula函数仅能够适用于变量存在正相关的联合分布；
[0022] (2-2)Clayton Copula函数
[0023] C(x1,x2)＝(x1-θ+x2-θ-1)-1/θ，θ∈(0,∞) (3)
[0024] 与Gumbel Copula函数一样，Clayton Copula函数仅能够适用于变量存在正相关的联合分布；
[0025] (2-3)Frank Copula函数
[0026]
[0027] Frank Copula函数对变量的相关性没有要求，正负相关的随机变量都能用它来描述；
[0028] (2-4)正态Copula函数
[0029]
[0030] 其中，Φ-1(.)是标准一元正态分布函数Φ(.)的逆函数，ρ为Φ-1(u)和Φ-1(v)的线性相关系数，r,s为积分变量，下同；正态Copula函数具有对称性，无法捕捉到非对称的相关关系；
[0031] (2-5)t-Copula函数
[0032]
[0033] 其中，Tv-1(.)是
自由度为ν的一元t分布函数Tv(.)的逆函数，二元t-Copula函数也具有对称性。
[0034] B2)对Copula函数进行的参数确定、拟合检验与评价；
[0035] 对Copula函数的参数估计方法采用矩估计法或极大似然估计法；
[0036] 理论上，对于Copula函数的拟合检验来说，传统的单变量检验方法都适合；对Copula函数的拟合检验方法采用Kolmogorov-Smirnov检验(柯尔莫哥洛夫-斯米尔诺夫检，简称：K-S检验)或Q-Q图；K-S检验是一种常用的非参数检验方法，二维检验统计量D的定义如下式(7)所示：
[0037]
[0038] 式中，n为样本个数，k为评价样本中的数据个数，Ck为联合观测样本xk＝(x1k,x2k)的Copula值，mk为联合观测值样本中满足条件x≤xk的联合观测值个数。
[0039] Q-Q图”即“分位数-分位数图”，它可直观的表达变量的经验频率与理论频率的拟合情况。
[0040] 评价Copula函数拟合程度的方法主要有三种：离差平方和准则(OLS)法、AIC信息准则法及Genest-Rivest法，本发明采用离差平方和准则(OLS)法，其公式为：
[0041]
[0042] 式(8)中，k为评价样本中的数据个数，pi和pei分别为联合分布的理论频率和经验频率。
[0043] 上述一种河口沿程淡水保证率的反演方法，步骤C中所建立的三维水动力和盐度数学模型采用MIKE3、EFDC、FVCOM中的一种，其控制方程为：
[0044] 连续性方程
[0045]
[0046] x和y方向的水平动量方程
[0047]
[0048]
[0050]
[0051] 式中t为时间(s)；x、y、z为笛卡尔坐标(m)；η为水面高程(m)；d为静水深(m)；h＝η+d为总水深(m)；u、v、w为x、y、z方向的流速分量(m/s)； 为
科里奥利力(N)；Ω为旋转
角速度，为所在地区纬度；g为重力
加速度(m2/s)；ρ为水
密度(kg/m3)；sxx、sxy、syx、syx为
辐射应力张量分量；νt为垂直湍粘滞系数(m2/s)；pa为
大气压(pa)；ρo为参考密度(m3/s)；S为源项(s-1)；s为盐度(‰)；us、vs为源项进入临近水体的流速分量(m/s)，Fu、Fv、Fs分别为x方向水平动量、y方向水平动量以及盐度的扩散项；ss为源项盐度(‰)。
[0052] 上述一种河口沿程淡水保证率的反演方法，步骤E中所述的临界径流-潮位控制线采用若干种工况组合试算的方法，具体的工况组合为径流、潮位的逐项组合。
[0053] 上述一种河口沿程淡水保证率的反演方法，步骤F所述的河口沿程淡水保证率采用
迭代求解的方式计算得到。
[0054] 本发明的有益效果为：
[0055] (1)本发明为一种耦合径流潮位联合分布和数值模型，并考虑人类取水需求的河口沿程淡水保证率反演方法，是针对河口径流和潮位实测序列较完整，而长序列盐度实测资料往往欠缺的现状不足，提出的一种基于历史径流、潮位资料、河口地形的盐水入侵风险评估方法，可用于任何具有长期历史数据的河口，具有良好的移植性和适用性；
[0056] (2)本发明的沿程淡水保证率反演方法有利于精准量化河口各河段的沿程淡水保证率，能够实现河口沿程淡水保证率的精准量化，为河口盐水入侵风险的确定提供技术支持，适于在河口地区的淡水资源利用中推广，十分方便进行河口取水口的合理选址，具有良好的应用性；
[0057] (3)本发明可通过计算机语言编译，即三维水动力和盐度数学模型采用MIKE3、EFDC、FVCOM中的任意一种工程
软件来构建，能够平面直观展示河口的盐水入侵风险，也可输出为文本表格。
附图说明
[0058] 下面结合附图和具体实施方式来详细说明本发明：
[0060] 图2是本发明实施例一中所建立的三维水动力和盐度数学模型；
[0061] 图3是本发明实施例一中径流、潮位逐项组合下的盐水入侵距离；
[0062] 图4是本发明实施例一中河口各河段临界径流-潮位控制线；
[0063] 图5是本发明实施例一中不同取水规模下各河段淡水保证率。
具体实施方式
[0064] 为使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施方式，进一步阐述本发明。
[0065] 如图1至图5所示，本发明的河口沿程淡水保证率的反演方法，包括以下步骤：
[0066] A.根据河口历史淡水径流和潮位资料，对河口逐日水文和潮汐特征进行解析，基于常用且典型的分布函数(例如：线性、均匀、正态、指数、PIII，Gompertz分布函数等)对径流和潮位累计频率曲线进行拟合，得出径流、潮位的边缘分布特性，并进一步考虑人为取水的影响，得到人为取水影响后的径流分布函数，取水后的径流分布函数即为天然径流方程发生平移，平移的单元数与取水规模相对应；
[0067] B.根据步骤A所得的河口径流、潮位的边缘分布特性，考虑两者之间的相互影响和关联特性，建立径流-潮位二维联合分布及其联合风险概率模型；
[0068] C.基于河口实测地形，构建三维水动力和盐度数学模型，三维水动力和盐度数学模型上游边界为淡水径流，下游边界为外海潮位；
[0069] D.基于步骤C所得的三维水动力和盐度数学模型，根据河口径流和潮位分布范围，以径流、潮位的逐项组合来设置计算工况，进行各工况下河口盐水入侵三维数值模拟，以径流、潮位的逐项组合来进行盐水入侵模拟，可得径流、潮位逐项组合下的最大盐水入侵距离(以盐度为0.45‰作为盐水入侵阈值)，反推得不同盐水入侵距离下的临界径流-潮位控制线，即河口各河段沿程发生盐水入侵的临界径流-潮位控制线，并得各河段不发生盐水入侵的径流和潮位区间，即淡水保证区；
[0070] E.基于步骤B得出的径流-潮位二维联合分布和步骤D得出的河口各河段沿程临界径流-潮位控制线，计算出盐水入侵的河口沿程淡水保证率；
[0071] F.根据步骤E得出的河口沿程淡水保证率，对比河口取水的设计供水保证率，得出不同取水规模下的临界取水断面，即取水口的临界布置位置。
[0072] 进一步地，步骤B中的径流-潮位二维联合分布采用Copula函数，并包括以下步骤：
[0073] B1)令H为具有边际分布F和G的联合分布函数，那么存在一个Copula函数C，满足：
[0074] H(x1,x2)＝C(F(x1),G(x2)) (1)
[0075] 其中，x1,x2为边缘随机变量，F为变量x1的边缘分布函数，G为变量x2的边缘分布函数。
[0076] 常用的二元Copula函数总体上可以分为三类，二元正态Copula函数、二元t-Copula函数和阿基米德Copula函数；其中含有一个参数的阿基米德Copula函数应用最广，包括Gumbel Copula函数、Clayton Copula函数和Frank Copula函数。因此，步骤B1中的径流-潮位二维联合分布采用Copula函数，选用Gumbel Copula函数、Clayton Copula函数、Frank Copula函数、二元正态Copula函数以及二元t-Copula函数中的任意一种或者任意几种组合使用，几种二元Copula函数如下所示：
[0077] (2-1)Gumbel Copula函数
[0078] C(x1,x2)＝exp{-[(-lnx1)θ+(-lnx2)θ]1/θ}，θ∈[1,∞) (2)
[0079] 式中，θ为相关参数，x1,x2为边缘随机变量；下同；由式(2)可知，Gumbel Copula函数仅能够适用于变量存在正相关的联合分布；
[0080] (2-2)Clayton Copula函数
[0081] C(x1,x2)＝(x1-θ+x2-θ-1)-1/θ，θ∈(0,∞) (3)
[0082] 与Gumbel Copula函数一样，Clayton Copula函数仅能够适用于变量存在正相关的联合分布；
[0083] (2-3)Frank Copula函数
[0084]
[0085] Frank Copula函数对变量的相关性没有要求，正负相关的随机变量都能用它来描述；
[0086] (2-4)正态Copula函数
[0087]
[0088] 其中，Φ-1(.)是标准一元正态分布函数Φ(.)的逆函数，ρ为Φ-1(u)和Φ-1(v)的线性相关系数，r,s为积分变量，下同；正态Copula函数具有对称性，无法捕捉到非对称的相关关系；
[0089] (2-5)t-Copula函数
[0090]
[0091] 其中，Tv-1(.)是自由度为ν的一元t分布函数Tv(.)的逆函数，二元t-Copula函数也具有对称性。
[0092] B2)对Copula函数进行的参数确定、拟合检验与评价；
[0093] 对Copula函数的参数估计方法采用矩估计法或极大似然估计法；
[0094] 理论上，对于Copula函数的拟合检验来说，传统的单变量检验方法都适合；对Copula函数的拟合检验方法采用Kolmogorov-Smirnov检验(柯尔莫哥洛夫-斯米尔诺夫检，简称：K-S检验)或Q-Q图；K-S检验是一种常用的非参数检验方法，二维检验统计量D的定义如下式(7)所示：
[0095]
[0096] 式中，n为样本个数，k为评价样本中的数据个数，，Ck为联合观测样本xk＝(x1k,x2k)的Copula值，mk为联合观测值样本中满足条件x≤xk的联合观测值个数。
[0097] Q-Q图”即“分位数-分位数图”，它可直观的表达变量的经验频率与理论频率的拟合情况。
[0098] 评价Copula函数拟合程度的方法主要有三种：离差平方和准则(OLS)法、AIC信息准则法及Genest-Rivest法，本发明采用离差平方和准则(OLS)法，其公式为：
[0099]
[0100] 式(8)中，k为评价样本中的数据个数，pi和pei分别为联合分布的理论频率和经验频率。
[0101] 再进一步地，步骤C中所建立的三维水动力和盐度数学模型采用MIKE3、EFDC、FVCOM中的一种，其控制方程为：
[0102] 连续性方程
[0103]
[0104] x和y方向的水平动量方程
[0105]
[0106]
[0107] 盐度对流扩散方程
[0108]
[0109] 式中t为时间(s)；x、y、z为笛卡尔坐标(m)；η为水面高程(m)；d为静水深(m)；h＝η+d为总水深(m)；u、v、w为x、y、z方向的流速分量(m/s)； 为科里奥利力(N)；Ω为旋转角速度，为所在地区纬度；g为
重力加速度(m2/s)；ρ为水密度(kg/m3)；sxx、sxy、syx、syx为辐射应力张量分量；νt为垂直湍粘滞系数(m2/s)；pa为大气压(pa)；ρo为参考密度(m3/s)；S为-1源项(s )；s为盐度(‰)；us、vs为源项进入临近水体的流速分量(m/s)，Fu、Fv、Fs分别为x方向水平动量、y方向水平动量以及盐度的扩散项；ss为源项盐度(‰)。
[0110] 再进一步地，步骤F的河口沿程淡水保证率采用迭代求解的方式计算得到。
[0111] 实施例一
[0112] 本实例以某河口的盐水入侵研究为例，分别基于Copula函数和三维数学模型进行了径流-潮位关联特性研究和盐水入侵敏感性分析，并根据不同取水规模引起的河口径流变化，进行动态环境下的沿程淡水保证率等不确定性研究。按照图1的简化流程图，具体实施步骤如下：
[0113] (1)对河口历史日径流和潮位数据进行频率分析，表明某河口径流累计频率曲线取对数后近似服从Gompertz分布，逐日高潮位分布近似服从均值μ为1.24,标准差σ为0.27的正态分布，
[0114]
[0115] x＝lg(q+M)
[0116] 其中q为径流(m3/s)，F为径流边缘分布函数，M为取水规模，当M＝0时为现阶段天然径流下的径流累计频率曲线。
[0117] h～N(μ,σ2)，其中μ＝1.24，σ＝0.27
[0118]
[0119] 其中h为潮位(m3/s)，G为潮位边缘分布函数，N为正态分布，基于Copula函数进行径流、潮位联合分布分析。根据逐日平均径流和潮位资料，选取较常用的Gumbel，Frank，Clayton和t-Copula函数分析多维水文联合分布规律，分别基于秩估计进行参数估计，得出四种Copula函数的参数估计值和拟合优度结果，如表1所示。
[0120] C(q,h)＝exp{-[(-lnq)1.28+(-lnh)1.28]1/1.28}
[0121] 表1二维变量Copula函数参数估计表
[0122]项目 Gumbel Frank Clayton t-Copula
参数大小 θ＝1.28 θ＝2.07 θ＝0.57 ρ＝0.34
AIC -105154 -98182 -90576 -99475
R2 0.9487 0.9122 0.8637 0.9155
[0123] (2)根据某河口地形资料，基于MIKE3模型构建了三维水动力和盐度数学模型，建模范围包括上游大坝至入海口之间的26km河段，如图2所示。上游边界为淡水径流，下游开边界为入海口水位，陆地边界采用已建堤防工程或护岸工程。基于河口地形情况，摩阻高度取值为0.002～0.005m，其中上游粗糙度大于下游，滩地粗糙度大于主槽。盐度水平扩散系数为0.12m3/s，垂向扩散系数为0.0001m3/s。计算计算时间步长采用2s。先加入流量、水位边界条件进行水动力计算，待河道流场稳定后加入盐度边界进行计算。待模拟区域盐度场稳定后，上下游边界根据计算工况进行变动，由此得出盐度场的时空变化分布。
[0124] (3)根据1980～2012年某河口逐日径流、高潮位序列，逐日径流范围为6.5～8670m3/s，逐日高潮位范围为-0.72～3.05m；计算工况采用以径流、潮位的逐项组合，总计
320组计算工况，涵盖了某河口全部逐日径流和高潮位范围，其中径流计算边界采用5～
9000m3/s，分别为5,10,20,40,…,200,250,…,500,600,…,1000,2000,…9000m3/s，共32项径流边界；潮位计算边界采用-1m～3.5m，分别为-1,-0.5,0,0.5,1,1.5,2,2.5,3,3.5m，共
10项潮位边界。
[0125] 根据径流、潮位的逐项组合的计算工况设置成果，以步骤2所得的三维数值模型进行模拟，得出盐度场稳定后的盐度场时空分布特征，可得不同径流、潮位逐项组合下的盐水入侵距离，如图3所示。由此可得河口不同断面的临界径流-潮位控制线，如图4所示为分别距离河口5km,10km,15km,20km和25km处的临界径流-潮位控制线。
[0126] (4)基于潮位、径流边缘分布函数和已择优的Copula联合分布函数，计算淡水保证区的概率值，即得出不同取水规模下(5m3/s,10m3/s,15m3/s,20m3/s,25m3/s,30m3/s)各断面处的淡水保证率，见图5所示。
[0127] 根据不同提水规模下各河段淡水保证率，可得某河口不同取水保证需求(80％,85％,90％,95％,99％)下的临界取水断面，如表2所示。
[0128] 表2不同取水保证需求下的临界取水断面(km)
[0129]
[0130] 以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的
权利要求书及其等效物界定。
