CN104698507A - 一种采煤沉陷区水资源效应的定量方法 - Google Patents

Abstract

本发明提供一种采煤沉陷区水资源效应的定量方法。本发明采煤沉陷区水资源效应的定量方法，包括：获取采煤沉陷区中积水区的第一积水面积、第一积水水位；获取采煤沉陷区中积水区的第二积水面积、第二积水水位；根据第一积水面积与第二积水面积获得采煤沉陷区中积水区的平均积水面积，根据第一积水水位与第二积水水位获得采煤沉陷区中积水区的平均积水水位；根据平均积水面积、平均积水水位、地下水水位通过水量平衡关系定量采煤沉陷区水资源效应。本发明实施例不需要长时间大量开展各种观测实验才能推理出采煤沉陷区的水资源构成和水循环过程，节省了工作成本，降低了研究难度。

Description

一种采煤沉陷区水资源效应的定量方法

技术领域

本发明涉及水利技术，尤其涉及一种采煤沉陷区水资源效应的定量方法。

背景技术

在煤炭开采过程中，当采空区面积扩大到一定范围后，岩层移动发展到地表，地表从原有标高向下沉降，引起地表沉陷，从而在采空区上方地表形成一个比采空区面积大得多的沉陷区域。在一些地下水潜水位较高的地区，地表一旦沉陷则很容易形成常年积水区，引起地表沼泽化，影响农耕，对土地利用的破坏作用十分显著。

但采煤沉陷区中的水是可用的水资源，具有开发利用价值，这对解决地区水资源短缺问题意义十分重要。对采煤沉陷区的水资源进行开发利用，首先需要获取采煤沉陷区水资源效应的定量数据，包括其水资源构成和水循环过程。传统的方式是开展各种观测实验，包括对采煤沉陷区形成的积水区进行河道入流量、降水量、蒸发量、地下水位、积水水位等水文观测，并借助实验数据资料通过水文推理来研究采煤沉陷区的水循环通量的收支状况，从而获得对采煤沉陷区水资源构成和水循环过程的认识，并根据采煤沉陷区水源补给情况评价其是否具有工程可利用价值等。

然而采煤沉陷区的空间尺度一般都比较大，环境条件比较复杂多变，通过传统实验手段的方法进行水资源效应的定量其成本高，难度大。

发明内容

本发明提供一种采煤沉陷区水资源效应的定量方法，以克服现有技术中所存在的不足。

本发明实施例提供一种采煤沉陷区水资源效应的定量方法，包括：

获取所述采煤沉陷区中积水区的第一积水面积、第一积水水位；所述第一积水面积为时间步长初的积水面积，所述第一积水水位为时间步长初的积水水位；

获取所述采煤沉陷区中积水区的第二积水面积、第二积水水位；所述第二积水面积为时间步长末的积水面积，所述第二积水水位为时间步长末的积水水位；

根据所述第一积水面积与所述第二积水面积获得所述采煤沉陷区中积水区的平均积水面积，根据所述第一积水水位与所述第二积水水位获得所述采煤沉陷区中积水区的平均积水水位；

根据所述平均积水面积、所述平均积水水位、地下水水位通过水量平衡关系对所述采煤沉陷区水资源效应进行定量。

进一步的，所述采煤沉陷区还包括：未积水区；

所述水资源效应包括：所述积水区的水面降水量，所述采煤沉陷区的上游河道汇入流量，所述未积水区的产流流量，地下水向所述积水区的渗出补给流量，从所述未积水区渗出、并流入所述积水区的地下水流量，所述积水区的水面蒸发量，所述积水区的渗漏量，所述采煤沉陷区的人工取水量，所述采煤沉陷区的排泄流量。

进一步的，所述水量平衡关系为：

$V^n-V^{n-1}=\mathrm{\Δt}(P+Q_{\mathrm{si}}+\mathrm{rnf}+G_{\mathrm{gw}}^{\mathrm{in}}+G_{\mathrm{non}}^{\mathrm{in}})-\mathrm{\Δt}(E+G_{\mathrm{gw}}^{\mathrm{out}}+W+Q_{\mathrm{so}})$

其中，Vⁿ为第二积水水量，所述第二积水水量为所述第二积水面积与所述第二积水水位对应的积水量；V^n-1为第一积水水量，所述第一积水水量为所述第一积水面积与所述第一积水水位对应的积水量；Δt为时间步长；P为所述积水区的水面降水量；Q_si为所述采煤沉陷区的上游河道汇入流量；rnf为所述未积水区的产流流量；为地下水向所述积水区的渗出补给流量；为从所述未积水区渗出，并流入所述积水区的地下水流量；E为所述积水区的水面蒸发量；为所述积水区的渗漏量；W为所述采煤沉陷区的人工取水量；Q_so为所述采煤沉陷区的排泄流量。

进一步的，所述第二积水面积、所述第二积水水位为通过所述第一积水面积、所述第一积水水位和时间步长内所述采煤沉陷区的水量平衡关系迭代计算后的收敛值。

进一步的，所述获取所述采煤沉陷区中积水区的第一积水面积、第一积水水位；所述第一积水面积为时间步长初的积水面积，所述第一积水水位为时间步长初的积水水位之前，还包括：

对所述采煤沉陷区进行空间离散，形成多个采煤沉陷区差分单元；对所述采煤沉陷区进行时间离散；

所述获取所述采煤沉陷区中积水区的第一积水面积；所述第一积水面积为时间步长初的积水面积，具体为：

获取所述每个采煤沉陷区差分单元的第一积水面积；所述第一积水面积为时间步长初的所述每个采煤沉陷区差分单元的积水面积之和；

所述获取所述采煤沉陷区中积水区的第二积水面积；所述第二积水面积为时间步长末的积水面积，具体为：

获取所述每个采煤沉陷区差分单元的第二积水面积；所述第二积水面积为时间步长末的所述每个采煤沉陷区差分单元的积水面积之和。

进一步的，所述对所述采煤沉陷区进行空间离散，包括：

对所述采煤沉陷区的地表水系统空间离散和对所述采煤沉陷区的地下水系统空间离散。

进一步的，所述根据所述第一积水面积与所述第二积水面积获得所述采煤沉陷区中积水区的平均积水面积，根据所述第一积水水位与所述第二积水水位获得所述采煤沉陷区中积水区的平均积水水位，具体为：

根据获得所述积水区的平均积水水位，其中，为所述积水区的平均积水水位；为所述第一积水水位；为所述第二积水水位，β为隐式加权因子，取值为0-1之间；

根据获得所述积水区的平均积水面积，其中，为所述积水区的平均积水面积；为所述第一积水面积；为所述第二积水面积；β为隐式加权因子，取值为0-1之间。

进一步的，还包括：根据所述平均积水水位与所述采煤沉陷区的底部高程确定所述多个采煤沉陷区差分单元的积水状态；所述积水状态包括：完全积水状态、完全未积水状态、半积水状态；

所述半积水状态差分单元的积水面积比为：

$R_{\mathrm{pnd}}=({\stackrel{\‾}{A}}_w-\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{AFCell}}_k)/\underset{p=1}{\overset{P}{\mathrm{\Σ}}}\mathrm{AH}{\mathrm{Cell}}_p$

其中，为所述积水区的平均积水面积；K为处于完全积水状态差分单元的总个数；AFCell_k为第k个完全积水状态差分单元的面积；P为处于半积水状态差分单元的总个数；AHCell_p为第p个半积水状态差分单元的面积。

进一步的，所述地下水水位通过三维地下水动力学方程获得；

所述三维地下水动力学方程为：

$\frac{\∂}{\∂x}(K_{\mathrm{xx}}\·\frac{\∂h_a}{\∂x})+\frac{\∂}{\∂y}(K_{\mathrm{yy}}\·\frac{\∂h_a}{\∂y})+\frac{\∂}{\∂z}(K_{\mathrm{zz}}\·\frac{\∂h_a}{\∂z})-W=S_s\frac{\∂h_a}{\∂t};$

其中，h_a为所述地下水水位；K_xx，K_yy和K_zz为渗透系数在X，Y和Z方向上的分量；W为单位体积流量；S_s为孔隙介质的贮水率；t为时间。

进一步的，所述积水区的水面降水量，所述采煤沉陷区的上游河道汇入流量，所述未积水区的产流流量，地下水向所述积水区的渗出补给流量，从所述未积水区渗出、并流入所述积水区的地下水流量，所述积水区的水面蒸发量，所述积水区的渗漏量，所述采煤沉陷区的人工取水量，所述采煤沉陷区的排泄流量具体为：

所述积水区的水面降水量为：

$P={\stackrel{\‾}{A}}_w\·P_{\mathrm{day}}^w$

其中，P为所述积水区的水面降水量；为所述积水区的平均积水面积；为所述积水区单位面积内的水面降水强度；

所述采煤沉陷区的上游河道汇入流量为：

$Q_{\mathrm{si}}=\mathrm{\Σ}{Q}_{\mathrm{si}}^i$

其中，Q_si为所述采煤沉陷区的上游河道汇入流量；第i个上游河道的汇入流量；

所述未积水区的产流流量为：

$\mathrm{rnf}=\frac{{(P_{\mathrm{day}}^{\mathrm{nw}}-0.2S)}^2}{(P_{\mathrm{day}}^{\mathrm{nw}}+0.8S)}\×{\stackrel{\‾}{A}}_{\mathrm{nw}},{\stackrel{\‾}{A}}_{\mathrm{nw}}=A_r-{\stackrel{\‾}{A}}_w;$

其中，rnf为所述未积水区的产流流量；为所述未积水区的平均面积；为所述未积水区单位面积内的陆面降水强度；S为土壤的最大可能滞留量；A_T为所述采煤沉陷区的总面积；为所述积水区的平均积水面积；

地下水向所述积水区的渗出补给流量为：

$Q_{\mathrm{gw}}^{\mathrm{in}}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{q}_k^{\mathrm{in}}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{C}_{s,k}(h_{a,k}-{\stackrel{\‾}{h}}_l);$

其中，为地下水向所述积水区的渗出补给流量；K为所述积水区的地下水向采煤沉陷区渗出的差分单元总数；为第k个积水区的差分单元处地下水向所述采煤沉陷区的渗出补给量；C_s,k为第k个积水区的差分单元处采煤沉陷区与地下水含水层间的导水系数；h_a,k为第k个积水区的差分单元处的地下水水位；为所述积水区的平均积水水位；

从所述未积水区渗出，并流入所述积水区的地下水流量为：

$G_{\mathrm{non}}^{\mathrm{in}}=\underset{m}{\overset{M}{\mathrm{\Σ}}}{q}_m^{\mathrm{in}}=\underset{m}{\overset{M}{\mathrm{\Σ}}}{C}_{s,m}(h_{a,m}-{\mathrm{Btm}}_m);$

其中，为从所述未积水区渗出，并流入所述积水区的地下水流量；M为所述未积水区的地下水渗出的差分单元总数；为第m个未积水区的差分单元处地下水向采煤沉陷区的渗出补给量；C_s,m为第m个未积水区的差分单元处采煤沉陷区与地下水含水层间的导水系数；h_a,m为第m个积水区的差分单元处的地下水水位；Btm_m为第m个未积水区的差分单元处的底部高程；

所述积水区的水面蒸发量为：

$E={\stackrel{\‾}{A}}_w\·{\mathrm{ET}}_{\mathrm{day}}^w;$

其中，E为所述积水区的水面蒸发量；为所述积水区的平均积水面积；为所述积水区单位面积内的水面蒸发强度；

所述积水区的渗漏量为：

$Q_{\mathrm{gw}}^{\mathrm{out}}=Q_{\mathrm{gw}}^{\mathrm{out},1}+Q_{\mathrm{gw}}^{\mathrm{out},2};$

其中， $Q_{\mathrm{gw}}^{\mathrm{out},1}=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{q}_n^{\mathrm{out}}=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{C}_{s,n}({\stackrel{\‾}{h}}_l-h_{a,n});$

为地下水水位高于所述采煤沉陷区的底部高程时所述积水区向所述地下水含水层的渗漏量；N为所述差分单元的总数；为第n个积水区的差分单元处采煤沉陷区的渗漏量；C_s,n为第n个积水区的差分单元处采煤沉陷区与地下水含水层间的导水系数；h_a,n为第n个积水区的差分单元处的地下水水位；为所述积水区的平均积水水位；

其中， $Q_{\mathrm{gw}}^{\mathrm{out},2}=\underset{s=1}{\overset{S}{\mathrm{\Σ}}}{q}_s^{\mathrm{out}}=\underset{s=1}{\overset{S}{\mathrm{\Σ}}}{C}_{s,s}({\stackrel{\‾}{h}}_l-{\mathrm{Btm}}_s);$

为地下水水位低于所述采煤沉陷区的底部高程时所述积水区向所述地下水含水层的渗漏量；S为所述差分单元的总数；为第s个积水区的差分单元处采煤沉陷区的渗漏量；C_s,s为第s个积水区的差分单元处采煤沉陷区与地下水含水层间的导水系数；Btm_s为第s个差分单元处的底部高程；为所述积水区的平均积水水位；

所述采煤沉陷区的人工取水量为：

W＝w_irri+w_indu+w_life+w_eco；

其中，W为所述采煤沉陷区的人工取水量；w_irri为采煤沉陷区的农业灌溉取水量；w_indu为所述采煤沉陷区的工业取水量；w_life为所述采煤沉陷区的生活取水量；w_eco为所述采煤沉陷区的生态环境取水量；

所述采煤沉陷区的排泄流量为：

$Q_{\mathrm{so}}=\mathrm{\σ\ϵm}{B}_0\sqrt{2}g{H}_0^{3/2},\mathrm{\σ}=2.31\frac{h_s}{H_0}{\{1-\frac{h_s}{H_0}\}}^{0.4};$

其中，Q_so为所述采煤沉陷区的排泄流量；σ为堰流淹没系数；ε为堰流侧收缩系数；m为堰流流量系数；B₀为闸孔总净宽；H₀计入行近流速水头的堰上水深；h_s为下游淹没水深；g为重力加速度；

本发明提供的一种采煤沉陷区水资源效应的定量方法，通过围绕采煤沉陷区的水量平衡关系，对步长时间步长内采煤沉陷区各项水量补给项和排泄项的计算，再结合时间步长初采煤沉陷区的积水水量，便可通过水量平衡关系式计算得到采煤沉陷区在时间步长末的积水水量。将上个时间步长末计算得出的积水水量作为当前时间步长初的积水水量，再通过逐时间步长演算的方式，就可以获得在长时间尺度下采煤沉陷区的积水水量变化过程，进而对其水资源效应进行评价，实现了仅需采用一般性的气象、降水及用水统计数据和水文地质数据，便能够对采煤沉陷区复杂多变的地表水与地下水转化关系进行自动识别和计算，从而对采煤沉陷区的水资源的构成和水循环过程进行定量评价。不需要长时间大量开展各种观测实验才能推理出采煤沉陷区的水资源构成和水循环过程，节省了研究成本，降低了研究难度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明采煤沉陷区水资源效应的定量方法实施例一的流程图；

图2为本发明采煤沉陷区水资源效应的定量方法中采煤沉陷区的平面示意图；

图3为本发明采煤沉陷区水资源效应定量方法实施例二的流程图；

图4为本发明采煤沉陷区水资源效应的定量方法中通过子流域格式离散后的采煤沉陷区的平面示意图；

图5为本发明采煤沉陷区水资源效应的定量方法中通过有限差分格式离散后的采煤沉陷区平面示意图；

图6为本发明采煤沉陷区水资源效应的定量方法中通过有限差分格式离散后的采煤沉陷区剖面示意图；

图7为本发明采煤沉陷区水资源效应的定量方法中采煤沉陷区差分单元积水状态示意图；

图8为本发明采煤沉陷区水资源效应的定量方法中地下水向采煤沉陷区中积水区渗出的示意图；

图9为本发明采煤沉陷区水资源效应的定量方法中采煤沉陷区差分单元的地下水向积水区渗出的示意图；

图10为本发明采煤沉陷区水资源效应的定量方法中地下水向采煤沉陷区中未积水区渗出的示意图；

图11为本发明采煤沉陷区水资源效应的定量方法中采煤沉陷区差分单元的地下水向未积水区渗出的示意图；

图12为本发明采煤沉陷区水资源效应的定量方法中地下水水位高于采煤沉陷区的底部高程时积水区向地下水含水层渗漏的示意图；

图13为本发明采煤沉陷区水资源效应的定量方法中地下水水位低于采煤沉陷区的底部高程时积水区向地下水含水层渗漏的示意图；

图14为本发明采煤沉陷区水资源效应的定量方法中某矿区及采煤沉陷区的子流域空间离散图；

图15为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区中积水区的水面降水量模拟计算数据图；

图16为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区的上游河道汇入流量模拟计算数据图；

图17为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区中未积水区的产流流量模拟计算数据图；

图18为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区中地下水向积水区的渗出补给流量模拟计算数据图；

图19为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区中从未积水区渗出、并流入积水区的地下水流量模拟计算数据图；

图20为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区中积水区的水面蒸发量模拟计算数据图；

图21为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区中积水区的渗漏量模拟计算数据图；

图22为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区的人工取水量模拟计算数据图；

图23为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区的排泄流量模拟计算数据图；

图24为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区的未积水区降水入渗量模拟计算数据图；

图25为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区的未积水区潜水蒸发量模拟计算数据图；

图26为本发明采煤沉陷区水资源效应的定量方法在水文模拟完成后对该矿区各个采煤沉陷区水资源效应的评价数据成果图；

图27为本发明采煤沉陷区水资源效应的定量方法中某矿区的第一采煤沉陷区在2007年汛期积水水位实测结果与本发明的模拟计算结果之间的对比图；

图28为本发明采煤沉陷区水资源效应的定量方法实施例三的流程图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

图1为本发明采煤沉陷区水资源效应的定量方法实施例一的流程图，如图1所示，本实施例的方法可以包括：

S101、获取所述采煤沉陷区中积水区的第一积水面积、第一积水水位；所述第一积水面积为时间步长初的积水面积，所述第一积水水位为时间步长初的积水水位。

在本发明中，采用某个矿区的采煤沉陷区来举例说明。图2为本发明采煤沉陷区水资源效应的定量方法中采煤沉陷区的平面示意图，如图2所示，该采煤沉陷区的面积分为两个部分，一部分面积当前覆盖有水面，称为积水区21，如图2中标记21所示；一部分面积为尚未积水的陆面，称为未积水区，如图2中标记22所示。首先任何时刻采煤沉陷区都具有统一的积水水位，即积水区的积水水位，该积水水位不考虑因上游河道的水量汇入动量、风浪、下泄等过程干扰引起的水位在积水区空间分布上的局部不均。采煤沉陷区的总面积在模拟计算过程中是固定的，为积水区面积和未积水区面积之和，但积水区面积和未积水区面积之间随着采煤沉陷区的积水水位的高低变化是可以动态相互转变的，采煤沉陷区的积水水位下降时，原先积水的部分面积露出地表，成为未积水区面积；而当采煤沉陷区的积水水位上升时，原先未积水的部分面积被水面淹没，而成为积水区面积。采煤沉陷区作为区域汇水范围内的低洼地带，可以自由接纳来自于上游多条河道的汇入水量，上游河道的如图2中标记23所示。采煤沉陷区的下游一般通过闸门控制，因此通常有一个或数量有限的几个主要的下泄通道，下泄通道如图2中标记24所示。

其中，第一积水面积、第一积水水位为在时间步长初的采煤沉陷区中积水区的积水面积、积水水位，这些水文数据可以通过文献收集、网络查询等方式取得。具体举例来说明，例如评价水资源效应有一个评价期，比如10年。在本发明中，10年的评价期将被分为若干个时间段(时间步长)，比如按1天一个时间步长的话，大约有3650个时间步长。第一积水面积、第一积水水位需要在评价期刚开始时(第一个时间步长的时段初)给出初值，这个初值是通过现场调查、遥感数据、文献收集、网络查询等方式取得的。后续每个时间步长的第一积水面积、第一积水水位，将通过迭代计算的方式确定，无需再给出。因为第一个时间步长的第二积水面积、第二积水水位为第二个时间步长的第一积水面积、第二积水水位，如此递推。

S102、获取所述采煤沉陷区中积水区的第二积水面积、第二积水水位；所述第二积水面积为时间步长末的积水面积，所述第二积水水位为时间步长末的积水水位。

根据步骤101中已知的第一积水面积、第一积水水位，采用迭代算法，得出时间步长末的收敛值，即第二积水面积、第二积水水位，时间步长末的积水面积、积水水位。

S103、根据所述第一积水面积与所述第二积水面积获得所述采煤沉陷区中积水区的平均积水面积，根据所述第一积水水位与所述第二积水水位获得所述采煤沉陷区中积水区的平均积水水位。

根据步骤101中已知的第一积水面积及步骤102中计算得到的第二积水面积，获得采煤沉陷区中积水区在时间步长内的平均积水面积。根据步骤101中已知的第一积水水位及步骤102中计算得到的第二积水水位，获得采煤沉陷区中积水区在时间步长内的平均积水水位。

S104、根据所述平均积水面积、所述平均积水水位、地下水水位通过水量平衡关系定量所述采煤沉陷区水资源效应。

其中，采煤沉陷区及周边地区的地下水水位通过现有成熟的技术可以求得，再将步骤103中获得的采煤沉陷区中积水区在时间步长内的平均积水面积、平均积水水位及地下水水位围绕水量平衡关系，进而获取在各个时间步长内的采煤沉陷区水资源效应。

本步骤可以用于对采煤沉陷区的积水过程有关的各水量分项进行分析，并得出描述采煤沉陷区的积水过程的数学控制方程，即水量平衡关系，作为采煤沉陷区的水文模拟的依据。根据采煤沉陷区的水文模拟原型，采煤沉陷区上游有多条河道的水量汇入；采煤沉陷区一部分面积积水、一部分未积水，那么积水区面积上将接受降水和蒸发，同时积水区还可以跟当地地下水发生相互作用，当地地下水可以对积水区进行补给，积水区的水量也可能发生渗漏；未积水区在降雨时将形成产流，并汇入积水区；未积水区的地下水也可能渗出地表，并汇入到积水区；采煤沉陷区的水量可以被人工引流利用，或通过闸坝工程下泄。

其中，模拟原型是指对将要进行水文模拟的系统进行的数学和物理概化，包括对其主要物理特征的描述和相关模拟计算条件的设定。在本发明中，主要涉及对采煤沉陷区基本水文结构特征的描述，以及对其积水动态过程的考虑。

本实施例，通过围绕采煤沉陷区的水量平衡关系，对在时间步长内采煤沉陷区各项水量补给项和排泄项的计算，再结合时间步长初采煤沉陷区的积水水量，便可通过水量平衡关系式计算得到采煤沉陷区在时间步长末的积水水量。将上个时间步长末计算得出的积水水量作为当前时间步长初的积水水量，再通过逐时间步长演算的方式，就可以获得在长时间尺度下采煤沉陷区的积水水量变化过程，进而对其水资源效应进行评价，实现了仅需采用一般性的气象、降水及用水统计数据和水文地质数据，便能够对采煤沉陷区复杂多变的地表水与地下水转化关系进行自动识别和计算，从而对采煤沉陷区的水资源的构成和水循环过程进行定量评价。从而不需要长时间大量开展各种观测实验才能推理出采煤沉陷区的水资源构成和水循环过程，节省了研究成本，降低了研究难度。

下面采用几个具体的实施例，对图1所示方法实施例的技术方案进行详细说明。

图3为本发明采煤沉陷区水资源效应定量方法实施例二的流程图，如图3所示，本实施例的方法可以包括：

S201、对所述采煤沉陷区进行空间离散，形成多个采煤沉陷区差分单元；对所述采煤沉陷区进行时间离散，形成多个时间片段单元。

对采煤沉陷区进行空间离散，以准确评价采煤沉陷区周边河道汇入过程，以及地下水在采煤沉陷区积水过程中的作用和影响。一方面采煤沉陷区周边河道可能有多条，而且各条河道的汇入时间、汇入水量均不同，需要通过分布式计算确定各条河道的汇入过程和水量。另一方面地下水在采煤沉陷区的积水机制中起到了重要的作用，要弄清这个问题，包括地下水在采煤沉陷区积水的通量，地下水水位对采煤沉陷区积水的调节作用等，将涉及到采煤沉陷区积水与地下水之间转化的计算问题。在本步骤中，地表水系统和地下水系统的计算分别采用不同的空间离散策略，并将采煤沉陷区的离散思路考虑在内，构成了面向采煤沉陷区水文模拟的整体离散技术思路。

空间离散化的目的是进行分布式计算，即将采煤沉陷区整体空间进行剖分，分成更小的计算单元，这些单元在水力参数上可以不同，在空间上具有水分流动关联关系，通过对每个单元逐个进行模拟再通过整体关联计算，可以取得比一般集总式模拟更为合理和准确的评价计算效果。由于采煤沉陷区的水文过程同时涉及地表水及地下水的影响，都需要进行分布式计算，而地表水分布式计算和地下水分布式计算各有特点和要求，此外还需考虑采煤沉陷区在其中的离散问题，因此需针对本发明的目的构建采煤沉陷区的地表水系统空间离散和所述采煤沉陷区的地下水系统空间离散。

地表水系统空间离散：地表水系统包括河道、水库等要素，通过空间的自然汇流关系形成地表水汇流网络。目前在对于地表水的分布式水文模拟方面，有多种空间离散方式，包括地形单元、差分单元和子流域单元三种。地形单元空间离散的地表水模型包括TOPMODEL，GBHM等，差分单元空间离散的地表水模型包括SHE等。相对地形单元和差分单元，由于子流域单元离散方式比较灵活，特别是对大尺度的水文模拟，具有计算效率高、尺度效应小的特点，应用较其他两种方式广泛得多。因此在本发明中，对应于地表水系统的空间离散，采用子流域格式的离散方式。

地下水系统空间离散：地下水分布式模拟方面目前主要为有限差分、有限元、有效体积等几种方式。在本发明中，为达到提高计算效率的目的，采用的是有限差分格式的离散方法。

采煤沉陷区既涉及地表水模拟过程，也涉及地下水模拟过程，因此其空间离散既要满足地表水汇流模拟的要求，又要满足地下水数值模拟的要求。从地表水离散的角度，本发明将单个采煤沉陷区离散为一个或多个特殊的子流域，采煤沉陷区的面积范围为这些子流域的面积范围。图4为本发明采煤沉陷区水资源效应的定量方法中通过子流域格式离散后的采煤沉陷区的平面示意图，如图4所示，其中，标记41表示一般的子流域单元，标记42表示子流域单元的边界，标记43表示一般子流域内的主河道，标记44表示采煤沉陷区范围内的子流域单元。采煤沉陷区范围内的子流域单元与其他一般子流域单元相同之处在于其属于采煤沉陷区的汇流系统，周边其他一般的子流域的主河道可以向这些子流域汇流。不同之处在于，这些子流域内部没有主河道，但包括水面(积水区)和陆面(未积水区)，而其他一般子流域包括陆面和主河道，因此这些子流域并不像其他一般子流域那样模拟降水产流过程和子流域内部汇流过程，而是模拟水面降水、水面蒸发、未积水区降水产流、未积水区潜水蒸发过程。从地下水离散角度，采煤沉陷区将被离散成为一系列的有限差分单元。该系列有限差分单元属于采煤沉陷区地下水有限差分单元的一部分，但采煤沉陷区所包含的有限差分单元具有特殊性，这类单元有可能处于积水状态，因此这类单元不仅包含有地下水的流动，还有采煤沉陷区积水和地下水之间的水量转化。对于一个具体的采煤沉陷区，图5为本发明采煤沉陷区水资源效应的定量方法中通过有限差分格式离散后的采煤沉陷区平面示意图。如图5所示，其中标记51为一般地下水有限差分单元，标记52为采煤沉陷区所包含的地下水有限差分单元。图6为本发明采煤沉陷区水资源效应的定量方法中通过有限差分格式离散后的采煤沉陷区剖面示意图。如图6所示，在垂向空间上模拟区域可能划分为多层有限差分单元，标记61为一般地下水有限差分单元，标记62为离散后各有限差分单元的层号。采煤沉陷区的底部高程位置是关键数据，因为它代表了采煤沉陷区的空间结构特征，图6中标记63表示了采煤沉陷区的实际底部形状和高程位置。在空间离散化的处理方式下，采煤沉陷区的洼底高程数据被离散到各个采煤沉陷区有限差分单元上，每个采煤沉陷区有限差分单元都具有其单元面积范围内的平均高程值，如图6中标记64所示。只有采煤沉陷区底部高程所在层位的有限差分单元才是采煤沉陷区有限差分单元，如标记65所示，标记66为该采煤沉陷区有限差分单元的所在的层位，如果采煤沉陷区有限差分单元不位于第一层，则其上方的所有有限差分单元都将被定义为无效单元，如图6中标记67所示，无效单元不参与地下水水位数值计算。

对所述采煤沉陷区进行时间离散，形成多个时间片段，具体为：用于将一段时期内的采煤沉陷区的水文模拟工作分为小的时间步长(或称时间片段)，进行逐时间步长的模拟。对采煤沉陷区水文规律的研究期一般至少以单个水文年为基础，有时甚至要研究数年乃至数十年的变化趋势。然而由于水分在不同介质内的流动具有显著的尺度效应和特定的运动规律，一般不能以太长的时间尺度作为模拟时间步长，否则会造成模拟失真。因此需要对研究期进行时间片段划分，也就是进行水文模拟的时间离散化处理。

水分在不同介质中的运动速度和规律不一样，一般而言，水分的运动速度越快，一般需要以更小的时间步长(或称时间片段)去进行模拟计算，否则就会造成模拟失真。地表水(河道、湖泊等)中水分的运动速度较快，因此需要以日、小时、分钟等时间步长才能合理描述其运动过程。一般仅对小尺度的渠道水力学过程，才需要采用小时或分钟这样的时间步长进行模拟计算，而对于流域尺度的河道汇流过程，由于一般汇流时间相对较长，采用日时间步长的情况较为常见；水分在土壤中运动的日动态效应比较明显，采用日时间步长是比较合适的；水分在地下水含水层中的运动比较缓慢，一般数日、旬或月才能觉察其状态的变化，因此可以采用较长的时间步长。本发明对地表水、土壤水及地下水均有涉及，需要选取三者中较小的时间步长作为统一的时间步长，即采用日时间步长对模拟期进行离散化。

S202、获取所述每个采煤沉陷区差分单元的第一积水面积、第一积水水位；所述第一积水面积为时间步长初的所述每个采煤沉陷区差分单元的积水面积之和，所述第一积水水位为时间步长初的采煤沉陷区的积水水位。

根据水文数据获取在时间步长初的每个采煤沉陷区差分单元的积水面积，即第一积水面积，获取在时间步长初的每个采煤沉陷区差分单元的积水水位，即第一积水水位。其中，水文数据都是一些较为普通的数据，可通过文献收集、网络查询等方式取得。

S203、获取所述每个采煤沉陷区差分单元的第二积水面积、第二积水水位；所述第二积水面积为时间步长末的所述每个采煤沉陷区差分单元的积水面积之和，所述第二积水水位为时间步长末的采煤沉陷区的积水水位。

具体的，所述第二积水面积、所述第二积水水位为通过所述第一积水面积、所述第一积水水位和时间步长内所述采煤沉陷区的水量平衡关系迭代计算后的收敛值。

S204、根据所述第一积水面积与所述第二积水面积获得所述采煤沉陷区中积水区的平均积水面积，根据所述第一积水水位与所述第二积水水位获得所述采煤沉陷区中积水区的平均积水水位。

在采煤沉陷区的积水水量确定之后，通过积水水位-积水水量关系以及积水面积-积水水量关系确定采煤沉陷区的积水水位和积水面积。模拟计算过程中由于对采煤沉陷区的底部高程进行了离散化处理，采煤沉陷区的积水水位-积水水量关系以及积水面积-积水水量关系将以散点关系分布。由于模拟计算过程中采煤沉陷区的积水水量是连续的数据，因此需在相邻散点之间考虑插值技术，从而可以根据采煤沉陷区的积水水量确定其积水水位和积水面积。

实际情况下，采煤沉陷区洼地底部是连续性变化的，因此积水水量与积水水位之间、积水水量与积水面积之间的数量关系也是连续变化的。而在采煤沉陷区被空间离散后，采煤沉陷区是通过一系列的差分单元来近似代表的，其积水水量-积水水位、积水水量-积水面积的数量关系被离散成为多个不连续的数据点。在模拟过程中积水水量有可能是任意值，当积水水量不在这些数据点上时，对应的积水水位和积水面积需要进行确定。此处的技术要点是进行线性插值，方法为，若通过空间离散关系得到的两个相邻的积水水位-积水水量离散点分别为(h₁,V(h₁))和(h₂,V(h₂))，若采煤沉陷区当前积水水量为V_n，其大小正处于V(h₁)和V(h₂)之间，那么可通过以下线性插值公式确定V_n对应的积水水位h_n：

$h_n=h_1+\frac{V_n-V\left(h_1\right)}{V\left(h_2\right)-V\left(h_1\right)}(h_2-h_1)$

同理若(S₁,V(S₁))以及(S₂,V(S₂))为通过网格关系得到的两个已知的积水面积-积水水量离散点，若采煤沉陷区当前积水水量为V_n，其大小处于V(S₁)和V(S₂)之间，那么可以通过以下线性插值公式确定V_n对应的积水面积S_n：

$S_n=S_1+\frac{V_n-V\left(S_1\right)}{V\left(S_2\right)-V\left(S_1\right)}(S_2-S_1)$

具体的，根据获得所述积水区的平均积水水位，其中，为所述积水区的平均积水水位；为所述第一积水水位；为所述第二积水水位，β为隐式加权因子，取值为0-1之间；

根据获得所述积水区的平均积水面积，其中，为所述积水区的平均积水面积；为所述第一积水面积；为所述第二积水面积；β为隐式加权因子，取值为0-1之间。

若β＝0，上述平均积水水位的公式与平均积水面积的公式为纯显式格式，该格式下采煤沉陷区在时间步长内的平均积水水位和积水面积以时间步长初的值为准，而时间步长初的平均积水面积与平均积水水位是已知的，所以在模拟过程中无需迭代试算，因而计算速度较快，但其计算的稳定性与时间步长的大小有关系，时间步长的时间尺度越小，则稳定性越好，反之则较差；若β＝1，该式为纯隐式格式，该格式下采煤沉陷区在时间步长内的平均积水水位和积水面积以时间步长末的值为准，而时间步长末的值是未知的，所以在模拟过程中需反复进行迭代试算，直到计算收敛为止。该格式的计算稳定性与时间步长的时间尺度无关，是无条件稳定的，但由于需反复迭代试算，因此计算速度较慢；在0<β<1时，该公式体现出显式和隐式的综合特点，是一种混合格式，随着β增大，计算过程从偏向显式逐渐过渡到偏向隐式，意味着当天湖泊的水平衡计算从主要依据时间步长初的积水区的积水水位和积水面积过渡到主要依据当天时间步长末的积水区的积水水位和积水面积。混合格式同样也需要反复进行迭代计算直到计算收敛，计算速度介于纯显式格式和纯隐式格式之间，同时在β大于0.5时，可以确保模拟计算的稳定性和具有较高的计算效率。

S205、根据所述平均积水水位与所述采煤沉陷区的底部高程确定所述多个采煤沉陷区差分单元的积水状态；所述积水状态包括：完全积水状态、完全未积水状态、半积水状态。

在平均积水水位下，被空间离散的各采煤沉陷区差分单元被淹没的状态不一样，因此在计算过程中应该区别对待。根据采煤沉陷区差分单元底部高程的离散数据和采煤沉陷区的平均积水水位数据，可将采煤沉陷区差分单元的积水状态分为三种，一种是完全积水状态的采煤沉陷区差分单元，二是完全未积水状态的采煤沉陷区差分单元，三是半积水状态的采煤沉陷区差分单元。

具体的，单个采煤沉陷区被离散成为一系列的差分单元。由于空间离散的原因，原本连续性的底部高程也被离散成间断的数据，每个采煤沉陷区差分单元的底部高程数据等于该差分单元面积范围内的底部高程平均值。根据常规的数值算法，如果用采煤沉陷区的平均积水水位和某个采煤沉陷区差分单元的底部高程来判断某个差分单元是否积水，其结果只有两种：在采煤沉陷区平均积水水位高于该差分单元的底部高程时，该差分单元处于完全积水状态；当低于该差分单元的底部高程时，该差分单元处于完全未积水状态。但根据本发明的研究，仅考虑采煤沉陷区的两种状态在模拟情况较为复杂时通常会造成数值震荡，从而导致模拟计算不收敛。其原因在于，采煤沉陷区的积水水量在模拟过程中是连续的数据，据此反算的采煤沉陷区的积水水位和积水面积根据线性插值也是连续的数据，而若采煤沉陷区差分单元只有两种状态，则根据采煤沉陷区差分单元统计的积水面积只能是间断的数据。因此当采煤沉陷区的积水水位因积水水量变化的原因跨越某个积水水量-积水水位离散点时，即使积水水位的变幅很小，根据采煤沉陷区差分单元统计的积水面积也将出现显著的间断变化。采采煤沉陷区的积水面积是影响采煤沉陷区中积水区的水面蒸发、地下水渗漏/补给等的重要参数，积水面积的显著变化会导致下次迭代中采煤沉陷区中积水区的积水水量发生显著变化，积水水量变化又引起采煤沉陷区中积水区的平均积水水位发生显著变化，这样模拟过程将形成围绕某个积水水位离散点的数值震荡，难以收敛。在本发明中，为了避免数值震荡发生，对采煤沉陷区差分单元的积水状态进行进一步处理，在完全未积水状态、完全积水状态之外增加了半积水状态，共有三种状态。以采煤沉陷区的平均积水水位所处的两个相邻底部高程离散点构成的区间作为判断基础，假设这两个相邻底部高程离散点的高程分别为H_Low和H_Up，底部高程小于H_Low的差分单元全部为完全积水状态单元，底部高程大于或等于H_Up的差分单元全部为完全未积水状态单元，底部高程等于H_Low的差分单元为半积水状态单元。图7为本发明采煤沉陷区水资源效应的定量方法中采煤沉陷区差分单元积水状态示意图。如图7所示，图中标记71表示一般的地下水差分单元，标记72表示采煤沉陷区的积水水位，标记73表示采煤沉陷区的实际底部高程，标记74表示空间离散后的采煤沉陷区差分单元上的底部高程，标记75表示处于完全未积水状态的采煤沉陷区差分单元，标记76表示处于半积水状态的采煤沉陷区差分单元，标记77表示处于完全积水状态的采煤沉陷区差分单元。半积水状态的采煤沉陷区差分单元需进一步确定其积水面积比，即该差分单元的积水面积占该单元总面积的比例，从而其积水面积部分采用与完全积水状态的采煤沉陷区差分单元相同的规则进行计算，其未积水面积的部分则采用与完全未积水状态的采煤沉陷区差分单元相同的规则进行计算。在本发明中，通过如下思路确定半积水状态的采煤沉陷区差分单元的积水面积比：采煤沉陷区的积水面积通过其积水水量-面积关系式是可以确定的，那么将积水面积扣除完全积水状态的采煤沉陷区差分单元的面积之和，剩下的积水面积等于半积水状态的采煤沉陷区差分单元的积水面积之和。将半积水状态的采煤沉陷区差分单元的积水面积之和除以半积水状态的采煤沉陷区差分单元的总面积，即可得半积水状态的采煤沉陷区差分单元的积水面积比。

具体的，所述半积水状态的采煤沉陷区差分单元的积水面积比为：

$R_{\mathrm{pnd}}=({\stackrel{\&OverBar;}{A}}_w-\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{AFCell}}_k)/\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}\mathrm{AH}{\mathrm{Cell}}_p$

其中，为所述积水区的平均积水面积；K为处于完全积水状态的采煤沉陷区差分单元的总个数；AFCell_k为第k个完全积水状态的采煤沉陷区差分单元的面积；P为处于半积水状态的采煤沉陷区差分单元的总个数；AHCell_p为第p个半积水状态的采煤沉陷区差分单元的面积。

本步骤中将采煤沉陷区各差分单元的积水状态分为完全积水状态、完全未积水状态、半积水状态，并进行针对性的处理，从而有效避的免了在计算模拟过程中的数值震荡现象。

S206、获取所述采煤沉陷区的地下水水位。

地下水分布式模拟的基本思想是把连续的模拟区域用有限个离散点构成的网格单元来代替，这些离散点称作网格单元的节点。通过把模拟区域上的连续变量的地下水三维动力学方程用在节点上定义的离散变量函数来近似，从而进行地下水的分布式计算。在具体应用时，网格单元可有多种形式的划分，包括有限差分法、有限单元法和有限体积法等。

无论是以上哪种方式的空间离散，地下水分布式模拟的原理都基于三维地下水动力学方程。

具体的，所述地下水水位通过三维地下水动力学方程获得。

所述三维地下水动力学方程为：

$\frac{\&PartialD;}{\&PartialD;x}(K_{\mathrm{xx}}\&CenterDot;\frac{\&PartialD;h_a}{\&PartialD;x})+\frac{\&PartialD;}{\&PartialD;y}(K_{\mathrm{yy}}\&CenterDot;\frac{\&PartialD;h_a}{\&PartialD;y})+\frac{\&PartialD;}{\&PartialD;z}(K_{\mathrm{zz}}\&CenterDot;\frac{\&PartialD;h_a}{\&PartialD;z})-W=S_s\frac{\&PartialD;h_a}{\&PartialD;t};$

其中，h_a为所述地下水水位；K_xx，K_yy和K_zz为渗透系数在X，Y和Z方向上的分量；W为单位体积流量；S_s为孔隙介质的贮水率；t为时间。

使用网格单元划分法对被模拟区域进行空间划分后，将以上三维地下水动力学方程应用到每个网格单元上，得一个线性方程组。这个方程组可用矩阵的形式表示为：

[A]{h_a}＝{q}；

其中[A]为水头的系数矩阵；{h_a}为所求的水位矩阵；而{q}表示各个方程中所包含的所有常数项和已知项。{q}有时也称为右侧项。再对所有网格单元的方程进行统一的求解，从而模拟地下水在空间上的流动过程。

以上矩阵形式的方程组的规模由网格单元的数量决定，有时甚至包含有上百万个未知数。这样通常需要采用迭代的方法进行求解，如强隐式法(Strongly Implicit Procedure，简称：SIP)、分层逐次超松弛法(SliceSuccessive Overrelaxation，简称：SSOR)或预调共轭梯度法(PreconditionedConjugate-gradient，简称：PCG)。

地下水数值计算的目的是预测采煤沉陷区的地下水水头变化，这取决于煤沉陷区的初始水头分布、煤沉陷区的边界条件、煤沉陷区的各种水文地质参数的分布以及各种外界源和汇的分布与强度。计算时总是从初始水头开始，每一步求出每个时间步长结束时的水头值，并用该值作为下一时间步长的初始值，不断重复这样的过程，直至所要求的时间结束为止。

S207、根据所述平均积水面积、所述平均积水水位、所述地下水水位通过水量平衡关系定量所述采煤沉陷区水资源效应。

其中，所述水资源效应包括：所述积水区的水面降水量，所述采煤沉陷区的上游河道汇入流量，所述未积水区的产流流量，地下水向所述积水区的渗出补给流量，从所述未积水区渗出、并流入所述积水区的地下水流量，所述积水区的水面蒸发量，所述积水区的渗漏量，所述采煤沉陷区的人工取水量，所述采煤沉陷区的排泄流量。

具体的，所述水量平衡关系为：

$V^n-V^{n-1}=\mathrm{\&Delta;t}(P+Q_{\mathrm{si}}+\mathrm{rnf}+G_{\mathrm{gw}}^{\mathrm{in}}+G_{\mathrm{non}}^{\mathrm{in}})-\mathrm{\&Delta;t}(E+G_{\mathrm{gw}}^{\mathrm{out}}+W+Q_{\mathrm{so}});$

其中，Vⁿ为第二积水水量，所述第二积水水量为所述第二积水面积与所述第二积水水位的乘积，即时间步长初的积水水量；V^n-1为第一积水水量，所述第一积水水量为所述第一积水面积与所述第一积水水位的乘积，即时间步长末的积水水量；Δt为时间步长；P为所述积水区的水面降水量；Q_si为所述采煤沉陷区的上游河道汇入流量；rnf为所述未积水区的产流流量；为地下水向所述积水区的渗出补给流量；为从所述未积水区渗出，并流入所述积水区的地下水流量；E为所述积水区的水面蒸发量；为所述积水区的渗漏量；W为所述采煤沉陷区的人工取水量；Q_so为所述采煤沉陷区的排泄流量。

下面逐一对上述的各个分项进行说明：

所述积水区的水面降水量：

是指直接降落到采煤沉陷区中积水区表面的降水量，该项水量通过时间步长内积水区的平均水面面积与降水量的乘积确定。所述积水区的水面降水量的计算公式为：

$P={\stackrel{\&OverBar;}{A}}_w\&CenterDot;P_{\mathrm{day}}^w;$

其中，P为所述积水区的水面降水量；为所述积水区的平均积水面积；为所述积水区单位面积内的水面降水强度。

所述采煤沉陷区的上游河道汇入流量：

是指采煤沉陷区周边与采煤沉陷区有汇流关系的河道的汇入量。该水量可由分布式水文模型进行产汇流计算得出。分布式水文模型是对自然界的水文现象进行抽象和概化，建立水文过程的数学结构与逻辑结构并进行近似计算表达的一种工具，进行产汇流计算是分布式水文模型的基本功能，可以通过降水和气象数据对研究范围内的每条河流的流量进行模拟和预测。同时分布式水文模型还可以分析空间上的汇流路径，确定与采煤沉陷区有关的具体汇流河道。将这些汇流河道的流量进行统计，可以得到时间步长内进入采煤沉陷区的总河道入流量。所述采煤沉陷区的上游河道汇入流量的计算公式为：

$Q_{\mathrm{si}}=\mathrm{\&Sigma;}{Q}_{\mathrm{si}}^i;$

其中，Q_si为所述采煤沉陷区的上游河道在时间步长内的汇入流量；第i个上游河道在时间步长内的汇入流量。

所述未积水区的产流流量：

是指在降雨过程中，在采煤沉陷区中未积水区上的降水强度大于该区土壤的入渗速度，来不及下渗的降水产生地表径流，并流入采煤沉陷区中积水区的水量。未积水区的产流汇入量。在本发明中，所述未积水区的产流流量的计算公式为：

$\mathrm{rnf}=\frac{{(P_{\mathrm{day}}^{\mathrm{nw}}-0.2S)}^2}{(P_{\mathrm{day}}^{\mathrm{nw}}+0.8S)}\&times;{\stackrel{\&OverBar;}{A}}_{\mathrm{nw}},{\stackrel{\&OverBar;}{A}}_{\mathrm{nw}}=A_r-{\stackrel{\&OverBar;}{A}}_w;$

其中，rnf为所述未积水区在时间步长内的产流流量；为所述未积水区在时间步长内的平均面积；为所述未积水区在时间步长内的单位面积内的陆面降水强度；S为土壤的最大可能滞留量；A_T为所述采煤沉陷区的总面积；为所述积水区在时间步长内的平均积水面积。

地下水向所述积水区的渗出补给流量：

是指采煤沉陷区中积水区的地下水渗出发生在地下水水位高于积水区的积水水位时，图8为本发明采煤沉陷区水资源效应的定量方法中地下水向采煤沉陷区中积水区渗出的示意图。如图8所示，图中标记81表示地面，标记82表示地下水含水层，标记83表示采煤沉陷区底部沉积物，标记84表示采煤沉陷区中积水区的积水水位，标记85表示地下水水位，标记86表示地下水向采煤沉陷区积水区渗出。采煤沉陷区的水量下泄，或人工大量抽取，或降水导致周边地下水位短时期内升高等，均会形成该种情形，此时地下水含水层与采煤沉陷区间的水位高差将驱使地下水从积水区底部渗出补给到采煤沉陷区。针对以上物理过程，进行分布式模拟时，地下水含水层向采煤沉陷区中积水区的地下水渗出发生在含水层的地下水水头高于采煤沉陷区中积水区的积水水位的采煤沉陷区差分单元处，图9为本发明采煤沉陷区水资源效应的定量方法中采煤沉陷区差分单元的地下水向积水区渗出的示意图。如图9所示，其中标记91代表采煤沉陷区差分单元，标记92为该单元的地下水水位，标记93为采煤沉陷区中积水区的积水水位，标记94为该单元上的采煤沉陷区的底部沉积物，标记95为底部高程，由于该单元上的地下水水位高于采煤沉陷区中积水区的积水水位，因此水分从地下水含水层向采煤沉陷区渗出，其方向示意为标记96。由于采煤沉陷区具有多个差分单元，因此在计算地下水含水层向采煤沉陷区中积水区的总渗出量时，需要对出现该情况的采煤沉陷区差分单元处的地下水渗出补给流量进行累加，地下水向所述积水区的渗出补给流量的计算公式为：

$Q_{\mathrm{gw}}^{\mathrm{in}}=\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{q}_k^{\mathrm{in}}=\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{C}_{s,k}(h_{a,k}-{\stackrel{\&OverBar;}{h}}_l);$

其中，为地下水向所述积水区的渗出补给流量；K为所述积水区的地下水向采煤沉陷区渗出的差分单元总数；为第k个积水区的差分单元处地下水向所述采煤沉陷区的渗出补给量；C_s,k为第k个积水区的差分单元处采煤沉陷区与地下水含水层间的导水系数；h_a,k为第k个积水区的差分单元处的地下水水位；为所述积水区的平均积水水位。

从所述未积水区渗出，并流入所述积水区的地下水流量：

是指采煤沉陷区中未积水区的地下水渗出发生在地下水水位高于未积水区的底部高程时，由于地下水的顶托作用，将会形成一定面积的渗出面，在渗出面处地下水渗出地表并流入采煤沉陷区中的积水区，类似于泉涌过程。图10为本发明采煤沉陷区水资源效应的定量方法中地下水向采煤沉陷区中未积水区渗出的示意图。如图10所示，标记101为采煤沉陷区附近的地下水水位，标记102为采煤沉陷区中积水区的积水水位，标记103为采煤沉陷区的沉积物，标记104为采煤沉陷区渗出面处的地下水渗出。通过对采煤沉陷区进行差分网格空间离散，基于以上物理过程可建计算概念图。图11为本发明采煤沉陷区水资源效应的定量方法中采煤沉陷区差分单元的地下水向未积水区渗出的示意图，如图11所示，标记111表示未积水状态的采煤沉陷区差分单元，标记112为积水状态的采煤沉陷区差分单元，标记113为未积水状态的采煤沉陷区差分单元的地下水水位，标记114为该未积水状态的采煤沉陷区差分单元的底部高程，标记115表示该未积水状态的采煤沉陷区差分单元的沉积物，标记116表示由于该单元处的地下水水位高于其底部高程，地下水将从该未积水状态的采煤沉陷区差分单元渗出，标记117表示渗出的水量流入采煤沉陷区中积水区。因为通常该过程较短，这里不考虑渗出的水量沿采煤沉陷区底部陆面的流经时间，认为该过程瞬时完成。未积水区的下水渗出量根据采煤沉陷区差分单元处的地下水水位与底部高程间的差值计算，该方法在地下水数值模拟研究中常用来处理自由排水边界。同样类似于积水区地下水渗出计算，需要对出现这种情况的采煤沉陷区差分单元进行水量统计，从而得出从所述未积水区渗出，并流入所述积水区的地下水流量，该计算公式为：

$G_{\mathrm{non}}^{\mathrm{in}}=\underset{m}{\overset{M}{\mathrm{\&Sigma;}}}{q}_m^{\mathrm{in}}=\underset{m}{\overset{M}{\mathrm{\&Sigma;}}}{C}_{s,m}(h_{a,m}-{\mathrm{Btm}}_m);$

其中，为从所述未积水区渗出，并流入所述积水区的地下水流量；M为所述未积水区的地下水渗出的差分单元总数；为第m个未积水区的差分单元处地下水向采煤沉陷区的渗出补给量；C_s,m为第m个未积水区的差分单元处采煤沉陷区与地下水含水层间的导水系数；h_a,m为第m个积水区的差分单元处的地下水水位；Btm_m为第m个未积水区的差分单元处的底部高程。

所述积水区的水面蒸发量：

是指水分从采煤沉陷区中积水区通过自然蒸散作用消耗的水量，不包括未积水区的陆面蒸发量。积水区的水面蒸发量通过在时间步长内的积水区的平均积水面积与水面蒸发强度的乘积来确定的，具体的计算公式为：

$E={\stackrel{\&OverBar;}{A}}_w\&CenterDot;{\mathrm{ET}}_{\mathrm{day}}^w;$

其中，E为所述积水区在时间步长内的水面蒸发量；为所述积水区在时间步长内的平均积水面积；为所述积水区在时间步长内的单位面积内的水面蒸发强度。

所述积水区的渗漏量：

在采煤沉陷区中积水区的平均积水水位高于地下水水位时，采煤沉陷区将发生渗漏。根据采煤沉陷区中积水区的平均积水水位、底部高程以及与地下水水位之间的关系，采煤沉陷区发生渗漏有两种形式，一种为与地下水水位有关的渗漏，一种与地下水水位无关的渗漏。本发明中，积水区的渗漏量的计算公式为：

$Q_{\mathrm{gw}}^{\mathrm{out}}=Q_{\mathrm{gw}}^{\mathrm{out},1}+Q_{\mathrm{gw}}^{\mathrm{out},2};$

其中，为地下水水位高于所述采煤沉陷区的底部高程时所述积水区向所述地下水含水层的渗漏量；为地下水水位低于所述采煤沉陷区的底部时所述积水区向所述地下水含水层的渗漏量。

图12为本发明采煤沉陷区水资源效应的定量方法中地下水水位高于采煤沉陷区的底部高程时积水区向地下水含水层渗漏的示意图。如图12所示，标记121表示采煤沉陷区中积水区的采煤沉陷差分单元，标记122表示该差分单元的地下水水位，高于该差分单元的底部高程，标记123表示采煤沉陷区中积水区的平均积水水位，标记124为该差分单元的底部高程，标记125为该差分单元上的沉积物。由于该差分单元的地下水水位低于采煤沉陷区中积水区的平均积水水位，且该差分单元的地下水水位高于底部高程，因此采煤沉陷区中积水区的积水向地下水含水层渗漏，如标记126所示。在这种情况，可由采煤沉陷区中积水区的平均积水水位、采煤沉陷区差分单元处的地下水水位、采煤沉陷洼区差分单元处的导水系数计算采煤沉陷洼区的渗漏量。同样由于数值离散，具有以上条件的采煤沉陷区差分单元可能有多个，计算时需累加，其计算公式为：

$Q_{\mathrm{gw}}^{\mathrm{out},1}=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{q}_n^{\mathrm{out}}=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{C}_{s,n}({\stackrel{\&OverBar;}{h}}_l-h_{a,n});$

为地下水水位高于所述采煤沉陷区的底部时所述积水区向所述地下水含水层的渗漏量；N为所述差分单元的总数；为第n个积水区的差分单元处采煤沉陷区的渗漏量；C_s,n为第n个积水区的差分单元处采煤沉陷区与地下水含水层间的导水系数；h_a,n为第n个积水区的差分单元处的地下水水位；为所述积水区在时间步长内的平均积水水位。

图13为本发明采煤沉陷区水资源效应的定量方法中地下水水位低于采煤沉陷区的底部高程时积水区向地下水含水层渗漏的示意图。地下水水位低于采煤沉陷区的底部高程时积水区向地下水含水层渗漏，此时采煤沉陷区的渗漏强度与地下水水位无关，而与底部高程和采煤沉陷区中积水区的平均积水水位有关，属于稳定渗漏。如图13所示，其中标记131表示采煤沉陷区中积水区的采煤沉陷差分单元，标记132表示该差分单元的地下水水位，且低于该底部高程，标记133表示采煤沉陷区中积水区的平均积水水位，标记134表示该差分单元的底部高程，标记135为该差分单元处的沉积物，标记136表示采煤沉陷区的水量向地下水含水层渗漏过程。在这种情况，可由采煤沉陷区中积水区的平均积水水位、部高程、采煤沉陷区差分单元处的导水系数计算采煤沉陷区的渗漏量，并对所有具有该种情况的采煤沉陷区差分单元进行累加，其计算公式为：

$Q_{\mathrm{gw}}^{\mathrm{out},2}=\underset{s=1}{\overset{S}{\mathrm{\&Sigma;}}}{q}_s^{\mathrm{out}}=\underset{s=1}{\overset{S}{\mathrm{\&Sigma;}}}{C}_{s,s}({\stackrel{\&OverBar;}{h}}_l-{\mathrm{Btm}}_s);$

为地下水水位低于所述采煤沉陷区的底部时所述积水区向所述地下水含水层的渗漏量；S为所述差分单元的总数；为第s个积水区的差分单元处采煤沉陷区的渗漏量；C_s,s为第s个积水区的差分单元处采煤沉陷区与地下水含水层间的导水系数；Btm_s为第s个差分单元处的底部高程；为所述积水区的平均积水水位。

所述采煤沉陷区的人工取水量：

采煤沉陷区作为平原蓄滞区，具有水量的拦蓄功能，可作为水源使用。在计算中采煤沉陷区的人工用水包括农业灌溉用水、工业用水、生活用水及生态用水等，这些用水量将作为输入数据直接给出。采煤沉陷区当天的总人工取用水量为这四种用途的人工取水量之和，其计算公式为：

W＝w_irri+w_indu+w_life+w_eco；

其中，W为所述采煤沉陷区的人工取水量；w_irri为采煤沉陷区的农业灌溉取水量；w_indu为所述采煤沉陷区的工业取水量；w_life为所述采煤沉陷区的生活取水量；w_eco为所述采煤沉陷区的生态环境取水量。

所述采煤沉陷区的排泄流量：

采煤沉陷区的下游排泄是采煤沉陷区水量平衡过程中受人工控制程度最高的一个环节。出于防洪、兴利拦蓄、挡潮、冲砂、下游用水等目的，采煤沉陷区一般通过建设泄水闸调控采煤沉陷洼区的积水水量或积水水位以满足生产需要。在本发明中，以平原区平底堰流公式作为采煤沉陷区下泄计算的基础方法，同时考虑水量拦蓄规则以及外水位条件等因素的影响。采煤沉陷区的排泄流量的计算公式为：

$Q_{\mathrm{so}}=\mathrm{\&sigma;\&epsiv;m}{B}_0\sqrt{2}g{H}_0^{3/2},\mathrm{\&sigma;}=2.31\frac{h_s}{H_0}{\{1-\frac{h_s}{H_0}\}}^{0.4};$

其中，Q_so为所述采煤沉陷区的排泄流量；σ为堰流淹没系数；ε为堰流侧收缩系数；m为堰流流量系数；B₀为闸孔总净宽；H₀计入行近流速水头的堰上水深，可以用采煤沉陷区中积水区在时间步长内的平均积水水位代替；h_s为下游淹没水深，如果采煤沉陷区以外的水位对采煤沉陷区的排泄有顶托作用，则以外水位作为下游淹没水深参数计算淹没系数，若无顶托，淹没系数σ为1；g为重力加速度；

另外，还要计算采煤沉陷区中未积水区的降水入渗量与潜水蒸发量。其目的是用于对采煤沉陷区中未积水区部分的水循环影响作用进行评价。采煤沉陷区中未积水区的面积是采煤沉陷区总面积的组成部分，尤其在非汛期采煤沉陷区的未积水区的面积有可能占有较大比例。作为采煤沉陷区中积水区的外部环境，其上的降水入渗过程与潜水蒸发过程会对积水区的水平衡产生影响，因此需要在计算过程中予以考虑。在本发明中，采煤沉陷区中未积水区的降水入渗量通过降水入渗补给系数法计算，未积水区的潜水蒸发量通过阿维里扬诺夫公式计算。

其中，未积水区降水入渗通过入渗补给系数法的计算公式为：

$G_{\mathrm{rec}}=A_{\mathrm{nw}}\&CenterDot;(P_{\mathrm{day}}^{\mathrm{nw}}-\mathrm{rnf})\&CenterDot;{\mathrm{coef}}_{\mathrm{rec}};$

A_nw为采煤沉陷区中未积水区在时间步长内的平均面积；为采煤沉陷区中未积水区在时间步长内的陆面降水强度；coef_rec为采煤沉陷区中未积水区的降水入渗补给系数；G_rec为采煤沉陷区中未积水区的降水入渗量，计算过程中会根据未积水区的采煤沉陷区差分单元分布情况将其进行空间展布。

单个未积水区的采煤沉陷区差分单元的潜水蒸发通过阿维里扬诺夫公式计算，其公式为：

$\left\{\begin{array}{cc}{G}_{\mathrm{evp},j}={\mathrm{ET}}_{\mathrm{day}}^w\&CenterDot;A_{\mathrm{nw},j}& h_{a,j}\&GreaterEqual;{\mathrm{Btm}}_j\\ G_{\mathrm{evp},j}={\mathrm{ET}}_{\mathrm{day}}^w\&CenterDot;A_{\mathrm{nw},j}\&CenterDot;{\left(\frac{{\mathrm{btm}}_j-h_{a,j}}{{\mathrm{Depth}}_{\max }}\right)}^2& ({\mathrm{Btm}}_j>h_{a,j})\&({\mathrm{Btm}}_j-h_{a,j}<{\mathrm{Depth}}_{\max })\\ G_{\mathrm{evp},j}=0& {\mathrm{Btm}}_j-h_{a,j}\&GreaterEqual;{\mathrm{Depth}}_{\max }\end{array}\right.;$

式中，为在时间步长内潜在的潜水蒸发强度，视为等于水面蒸发强度；Depth_max为未积水区的潜水蒸发极限埋深；p_evp为未积水区的潜水蒸发指数；btm_j为第j个未积水区的采煤沉陷区差分单元处的底部高程；G_evp,j为第j个未积水区的采煤沉陷区差分单元的潜水蒸发量；A_nw,j为第j个未积水区的采煤沉陷区差分单元的面积；h_a,j为第j个未积水区的采煤沉陷区差分单元的地下水水位。

采煤沉陷区未积水区的潜水蒸发量为各未积水区的采煤沉陷区差分单元潜水蒸发量的总和，其计算公式为：

$G_{\mathrm{evp}}=\underset{j=1}{\overset{J}{\mathrm{\&Sigma;}}}{G}_{\mathrm{evp},j};$

式中，J为未积水区的采煤沉陷区差分单元的总个数。

本实施例中，通过建立采煤沉陷区的水文模拟原型，对采煤沉陷区进行空间离散，形成多个采煤沉陷区差分单元，对采煤沉陷区进行时间离散，形成多个时间片段单元；建立采煤沉陷区的积水水量-积水水位-积水面积之间数据关系；计算采煤沉陷区在时间步长内平均积水面积和平均积水水位；判断各采煤沉陷区差分单元的积水状态；根据上述获得的平均积水面积、平均积水水位、地下水水位，将这些数据代入采煤沉陷区的水量平衡关系中，计算出采煤沉陷区水资源效应，即各项水量分项，从而对采煤沉陷区的水资源的构成和水循环过程进行定量评价。实现了仅需采用一般性的气象、降水及用水统计数据和水文地质数据，便能够对采煤沉陷区复杂多变的地表水与地下水转化关系进行自动识别和计算，并对采煤沉陷区的积水过程进行数值反演。不需要长时间大量开展各种观测实验才能推理出采煤沉陷区的水资源构成和水循环过程，节省了研究成本，降低了研究难度。同时，该方法的物理仿真程度高，计算分析速度快，能够综合采煤沉陷区自身和周边汇流区的水量过程进行影响关系分析，能够进行长期的应用。

为了让上述步骤更容易理解，在此以实际某矿区的研究来举例说明，图14为本发明采煤沉陷区水资源效应的定量方法中某矿区及采煤沉陷区的子流域空间离散图。该矿区含有6个主要的采煤沉陷区，如图14所示，其中标记141为第一采煤沉陷区，标记142为第二采煤沉陷区，标记143为第三采煤沉陷区，标记144为第四采煤沉陷区，标记145为第五采煤沉陷区，标记146为第六采煤沉陷区，与这些采煤沉陷区有水位过程和水量联系的研究区总面积为4208km²。为研究地表河道向采煤沉陷区的汇入过程，地表水系统空间离散方面，全研究区根据DEM汇流分析，共划分为1265个子流域，对应每个子流域都有自己的主河道，其中6个采煤沉陷区作为特殊的子流域单独进行了划分。地下水系统空间离散方面，全研究区以500m为间距进行了有限差分格式剖分，共剖分132行276列，分浅层和深层两层，单层有限差分单元的数量为36432个，研究区范围内的有效单元为单层16166个，6个采煤沉陷区均根据其沉陷面积范围识别了其所包含的地下水有限差分网格。

图15为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区中积水区的水面降水量模拟计算数据图。图16为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区的上游河道汇入流量模拟计算数据图。图17为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区中未积水区的产流流量模拟计算数据图。图18为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区中地下水向积水区的渗出补给流量模拟计算数据图。图19为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区中从未积水区渗出、并流入积水区的地下水流量模拟计算数据图。图20为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区中积水区的水面蒸发量模拟计算数据图。图21为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区中积水区的渗漏量模拟计算数据图。图22为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区的人工取水量模拟计算数据图。图23为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区的排泄流量模拟计算数据图。图24为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区的未积水区降水入渗量模拟计算数据图。图25为本发明采煤沉陷区水资源效应的定量方法中在2000-2010年间某矿区的第一采煤沉陷区的未积水区潜水蒸发量模拟计算数据图。

通过本实施例上述步骤的应用，以某矿区的第一采煤沉陷区为例，如图15-图25所示，详细给出了该采煤沉陷区从2001年到2010年长达10年(3652天)期间各项水量的逐日反演过程，包括该采煤沉陷区中积水区的水面降水量、该采煤沉陷区的上游河道汇入流量、该采煤沉陷区中未积水区的产流流量、该采煤沉陷区中地下水向积水区的渗出补给流量、该采煤沉陷区中从未积水区渗出、并流入积水区的地下水流量、该采煤沉陷区中积水区的水面蒸发量、该采煤沉陷区中积水区的渗漏量、该采煤沉陷区的人工取水量、该采煤沉陷区的排泄流量、该采煤沉陷区的未积水区降水入渗量、该采煤沉陷区的未积水区潜水蒸发量等共11项水量过程项。这些水量分项是本发明对研究区复杂的地表水与地下水转化关系进行自动识别和计算后产生的，是通过本发明的技术手段，对采煤沉陷区的水量平衡过程进行的符合自然规律的解析和反映。体现出了本发明具有对采煤沉陷区的积水过程相关的各项水量过程进行完整数值反演的能力，且能够在多年尺度进行长期应用。

本发明在研究该矿区的采煤沉陷区水资源效应的过程中，所采用的数据主要为3个国家气象站的气象观测资料，当地11个降水站的降雨资料，以及当地的水文地质数据和用水统计数据等。这些数据都是一些较为普通的水文数据，可通过文献收集、网络查询等方式取得，相对于针对单个采煤沉陷区进行专门实验观测的难度和所需成本，本发明所需数据资料的获取难度低、成本低，且能够适用于多个采煤沉陷区集合研究的情况，如本发明一次对该矿区的6个采煤沉陷区开展了研究应用，进一步降低了应用成本。此外本发明研究效率高，提出的技术方法经过编程可在普通台式机上运行。体现出本了应用成本低和研究效率高的特点。

图26为本发明采煤沉陷区水资源效应的定量方法在水文模拟完成后对该矿区各个采煤沉陷区水资源效应的评价数据成果图。经由各个采煤沉陷区的逐日水量分项统计而来，如图26所示，图中各个采煤沉陷区的水资源来水总量，地表及地下水来水组成及比例大小，水资源各项排泄量、初始蓄水量及模拟完成后的蓄水量等均一目了然，可为将来对采煤沉陷区的水资源开发利用提供工程参考数据。另外，图27为本发明采煤沉陷区水资源效应的定量方法中某矿区的第一采煤沉陷区在2007年汛期积水水位实测结果与本发明的模拟计算结果之间的对比图。如图27所示，其中标记271为实测的积水水位，标记272为模拟计算的积水水位。从对比结果来看，实测结果与模拟计算结果之间的吻合度很高，可见本发明的模拟精度是令人满意的，实践应用性很强。

图28为本发明采煤沉陷区水资源效应的定量方法实施例三的流程图。如图28所示，在上述实施例的基础上，本方法实施例的步骤可以包括：

S301、获取计算所需的数据。

如各个采煤沉陷区经空间离散后的积水水位-积水水量、积水面积-积水水量数据点，水文气象、水文地质数据等；根据输入数据的信息，设置计算模拟的起始日期，结束日期等。

S302、以日为时间步长单位，进行水文模拟计算。

从模拟计算期的第1天开始，以日为时间步长单位，逐日进行采煤沉陷区的水文模拟计算工作。

S303、地表水分布式模拟计算，计算时间步长的河道系统流量。

在每日的水文模拟计算过程中，分布式水文模拟计算先运行，计算研究区内的河道系统中各条河道的流量。

S304、地下水在时间步长内迭代循环计算。

地下水分布式模拟基于数值算法，同时采煤沉陷区中积水区的平均积水面积、平均积水水位、积水与地下水之间的交换通量等均需试算，因此需要在每个模拟计算时间步长内进行迭代循环，以逐步逼近满足精度要求的近似解。

S305、计算各采煤沉陷区中积水区的平均积水面积和平均积水水位。

每次迭代，都先需要根据采煤沉陷区中积水区的积水水量确定模拟计算时间步长内的平均积水水位和平均积水面积。首次迭代计算时，采煤沉陷区中积水区的平均积水水位和平均积水面积为上个时间步长末(本时间步长初)的计算结果。在首次迭代计算后，由于在时间步长末的采煤沉陷区的积水水量数据已经通过试算得到更新，因此可由时间步长初(上个时间步长末)和该次迭代计算的积水水位和积水面积确定本时间步长内的平均积水水位和平均积水面积。

S306、根据平均积水水位判断各采煤沉陷区差分单元所处积水状态。

根据采煤沉陷区中积水区的平均积水水位状况，采煤沉陷区差分单元积水状态包括三种，一是完全未积水状态的差分单元，二是半积水状态的差分单元、三是完全积水状态的差分单元。

S307、计算时间步长内所有与地下水无关的采煤沉陷区源汇项。

包括积水区的水面降水量，积水区的水面蒸发量，未积水区的产流流量，积水区的渗漏量，采煤沉陷区的上游河道汇入流量，采煤沉陷区的人工取水量和采煤沉陷区的排泄流量。

S308、建立三维地下水数值计算仿真矩阵方程。

先仅从地下水有限差分单元间渗流的角度构建基本的地下水数值计算仿真矩阵方程。

S309、处理除采煤沉陷区差分单元之外的源汇项和边界项。

这里除采煤沉陷区差分单元之外的源汇项主要是其他一般地下水差分单元的降水入渗量、潜水蒸发量、河道渗漏量等。这些都是发生在采煤沉陷区之外水量，但由于地下水沟通了包含采煤沉陷区在内的整个研究区，这些水量对采煤沉陷区的积水过程有直接和间接的影响，需要进行处理。

S310、处理采煤沉陷区差分单元的地下水源汇项。

根据该采煤沉陷区差分单元的积水状态，首先确定该采煤沉陷区差分单元的积水面积和未积水面积。完全未积水状态的差分单元只有未积水面积，完全积水状态的差分单元只有积水面积，半积水状态的差分单元则既有积水面积，又有未积水面积。确定该数据后，则可对所有采煤沉陷区差分单元进行逐个循环，依次处理有关地下水源汇项。

S311、处理积水状态面积部分的采煤沉陷区差分单元的地下水源汇项。

对于某个采煤沉陷区差分单元，积水面积部分的地下水源汇项包括该面积上的渗漏量(含与地下水水位无关稳定渗漏量和与地下水位有关的渗漏量)、该面积上的地下水渗出量。

S312、处理未积水状态面积部分的采煤沉陷区差分单元的地下水源汇项。

对于某个采煤沉陷区差分单元，未积水面积部分的地下水源汇项包括该面积上的潜水蒸发量、该面积上的地下水渗出量。

S313、求解矩阵方程，得到当前迭代下的地下水水位。

将采煤沉陷区各差分单元上的源汇项整合进步骤8中建立的地下水水位矩阵方程，并用求解大型矩阵方程的算法对方程进行近似求解。这些算法包括强隐式法、共轭梯度法等。在矩阵方程得到求解后，当前迭代循环下的地下水水位即得到更新。

S314、根据采煤沉陷区的水量平衡关系公式计算采煤沉陷区在时间步长末的积水水量。

在地下水水位值更新后，采煤沉陷区的源汇项中与地下水水位有关的分项，如地下水向积水区的渗出补给流量、从未积水区渗出、并流入积水区的地下水流量、潜水蒸发量等均能够得到确定。综合之前计算的与地下水水位无关的采煤沉陷区源汇项分项，可根据采煤沉陷区的水量平衡关系公式对当前迭代循环下的积水水量进行计算，并作为本时间步长末的积水水量试算结果。

S315、更新本次迭代循环后的采煤沉陷区的积水水位和积水面积。

根据时间步长末的采煤沉陷区的积水水量的试算结果，通过积水水量-积水水位关系式和积水水量-积水面积关系式更新本次迭代循环后的采煤沉陷区的积水水位和积水面积。

S316、判断迭代是否收敛。

由于是近似求解，有可能地下水水位的求解结果并不能满足事先规定的精度要求，因此需要进行判断。如果求解结果满足精度要求，则在本次时间步长的模拟计算完成，可以进入下一个时间步长进行模拟计算；若不满足精度要求，则返回S304这一步重新反复迭代试算，直至地下水水位结果满足精度为止。

S317、判断模拟期是否结束。

判断当前水文模拟计算的时间步长是否是模拟期的最后一个时间步长，如果是，则终止计算，否则返回到S302这一步，继续下一个时间步长的水文模拟计算，并逐时间步长推进，直至规定的模拟期结束为止。

S318、结束计算。

本实施例，通过建立采煤沉陷区的水文模拟原型，对采煤沉陷区进行空间离散，形成多个采煤沉陷区差分单元，对采煤沉陷区进行时间离散，形成多个时间片段单元；建立采煤沉陷区的积水水量-积水水位-积水面积之间数据关系；计算采煤沉陷区在时间步长内平均积水面积和平均积水水位；判断各采煤沉陷区差分单元的积水状态；根据上述获得的平均积水面积、平均积水水位、地下水水位，将这些数据代入采煤沉陷区的水量平衡关系中，计算出采煤沉陷区水资源效应，即各项水量分项，从而对采煤沉陷区的水资源的构成和水循环过程进行定量评价。实现了仅需采用一般性的气象、降水及用水统计数据和水文地质数据，便能够对采煤沉陷区复杂多变的地表水与地下水转化关系进行自动识别和计算，并对采煤沉陷区的积水过程进行数值反演。不需要长时间大量开展各种观测实验才能推理出采煤沉陷区的水资源构成和水循环过程，节省了研究成本，降低了研究难度。同时，该方法的物理仿真程度高，计算分析速度快，能够综合采煤沉陷区自身和周边汇流区的水量过程进行影响关系分析，能够进行长期的应用。

本领域普通技术人员可以理解：实现上述各方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成。前述的程序可以存储于一计算机可读取存储介质中。该程序在执行时，执行包括上述各方法实施例的步骤；而前述的存储介质包括：ROM、RAM、磁碟或者光盘等各种可以存储程序代码的介质。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (10)

1.一种采煤沉陷区水资源效应的定量方法，其特征在于，包括：

获取所述采煤沉陷区中积水区的第一积水面积、第一积水水位；所述第一积水面积为时间步长初的积水面积，所述第一积水水位为时间步长初的积水水位；

获取所述采煤沉陷区中积水区的第二积水面积、第二积水水位；所述第二积水面积为时间步长末的积水面积，所述第二积水水位为时间步长末的积水水位；

根据所述第一积水面积与所述第二积水面积获得所述采煤沉陷区中积水区的平均积水面积，根据所述第一积水水位与所述第二积水水位获得所述采煤沉陷区中积水区的平均积水水位；

根据所述平均积水面积、所述平均积水水位、地下水水位通过水量平衡关系定量所述采煤沉陷区水资源效应。

2.根据权利要求1所述的方法，其特征在于，所述采煤沉陷区还包括：未积水区；

所述水资源效应包括：所述积水区的水面降水量，所述采煤沉陷区的上游河道汇入流量，所述未积水区的产流流量，地下水向所述积水区的渗出补给流量，从所述未积水区渗出、并流入所述积水区的地下水流量，所述积水区的水面蒸发量，所述积水区的渗漏量，所述采煤沉陷区的人工取水量，所述采煤沉陷区的排泄流量。

3.根据权利要求2所述的方法，其特征在于，所述水量平衡关系为：

其中，Vⁿ为第二积水水量，所述第二积水水量为所述第二积水面积与所述第二积水水位对应的积水量；V^n-1为第一积水水量，所述第一积水水量为所述第一积水面积与所述第一积水水位对应的积水量；Δt为时间步长；P为所述积水区的水面降水量；Q_si为所述采煤沉陷区的上游河道汇入流量；rnf为所述未积水区的产流流量；为地下水向所述积水区的渗出补给流量；为从所述未积水区渗出，并流入所述积水区的地下水流量；E为所述积水区的水面蒸发量；为所述积水区的渗漏量；W为所述采煤沉陷区的人工取水 量；Q_so为所述采煤沉陷区的排泄流量。

4.根据权利要求3所述的方法，其特征在于，所述第二积水面积、所述第二积水水位为通过所述第一积水面积、所述第一积水水位和时间步长内所述采煤沉陷区的水量平衡关系迭代计算后的收敛值。

5.根据权利要求4所述的方法，其特征在于，所述获取所述采煤沉陷区中积水区的第一积水面积、第一积水水位；所述第一积水面积为时间步长初的积水面积，所述第一积水水位为时间步长初的积水水位之前，还包括：

对所述采煤沉陷区进行空间离散，形成多个采煤沉陷区差分单元；对所述采煤沉陷区进行时间离散；

所述获取所述采煤沉陷区中积水区的第一积水面积；所述第一积水面积为时间步长初的积水面积，具体为：

获取所述每个采煤沉陷区差分单元的第一积水面积；所述第一积水面积为时间步长初的所述每个采煤沉陷区差分单元的积水面积之和；

所述获取所述采煤沉陷区中积水区的第二积水面积；所述第二积水面积为时间步长末的积水面积，具体为：

获取所述每个采煤沉陷区差分单元的第二积水面积；所述第二积水面积为时间步长末的所述每个采煤沉陷区差分单元的积水面积之和。

6.根据权利要求5所述的方法，其特征在于，所述对所述采煤沉陷区进行空间离散，包括：

对所述采煤沉陷区的地表水系统空间离散和对所述采煤沉陷区的地下水系统空间离散。

7.根据权利要求6所述的方法，其特征在于，所述根据所述第一积水面积与所述第二积水面积获得所述采煤沉陷区中积水区的平均积水面积，根据所述第一积水水位与所述第二积水水位获得所述采煤沉陷区中积水区的平均积水水位，具体为：

根据获得所述积水区的平均积水水位，其中，为所述积水区的平均积水水位；为所述第一积水水位；为所述第二积水水位，β为隐式加权因子，取值为0-1之间；

根据获得所述积水区的平均积水面积，其中，为所述积水区的平均积水面积；为所述第一积水面积；为所述第二积水面 积；β为隐式加权因子，取值为0-1之间。

8.根据权利要求7所述的方法，其特征在于，还包括：根据所述平均积水水位与所述采煤沉陷区的底部高程确定所述多个采煤沉陷区差分单元的积水状态；所述积水状态包括：完全积水状态、完全未积水状态、半积水状态；

所述半积水状态差分单元的积水面积比为：

其中，为所述积水区的平均积水面积；K为处于完全积水状态差分单元的总个数；AFCell_k为第k个完全积水状态差分单元的面积；P为处于半积水状态差分单元的总个数；AHCell_p为第p个半积水状态差分单元的面积。

9.根据权利要求8所述的方法，其特征在于，所述地下水水位通过三维地下水动力学方程获得；

所述三维地下水动力学方程为：

其中，h_a为所述地下水水位；K_xx，K_yy和K_zz为渗透系数在X，Y和Z方向上的分量；W为单位体积流量；S_s为孔隙介质的贮水率；t为时间。

10.根据权利要求9所述的方法，其特征在于，所述积水区的水面降水量，所述采煤沉陷区的上游河道汇入流量，所述未积水区的产流流量，地下水向所述积水区的渗出补给流量，从所述未积水区渗出、并流入所述积水区的地下水流量，所述积水区的水面蒸发量，所述积水区的渗漏量，所述采煤沉陷区的人工取水量，所述采煤沉陷区的排泄流量具体为：

所述积水区的水面降水量为：

其中，P为所述积水区的水面降水量；为所述积水区的平均积水面积； 为所述积水区单位面积内的水面降水强度；

所述采煤沉陷区的上游河道汇入流量为：

其中，Q_si为所述采煤沉陷区的上游河道汇入流量；第i个上游河道的汇入流量；

所述未积水区的产流流量为：

其中，rnf为所述未积水区的产流流量；为所述未积水区的平均面积； 为所述未积水区单位面积内的陆面降水强度；S为土壤的最大可能滞留量；A_T为所述采煤沉陷区的总面积；为所述积水区的平均积水面积；

地下水向所述积水区的渗出补给流量为：

其中，为地下水向所述积水区的渗出补给流量；K为所述积水区的地下水向采煤沉陷区渗出的差分单元总数；为第k个积水区的差分单元处地下水向所述采煤沉陷区的渗出补给量；C_s,k为第k个积水区的差分单元处采煤沉陷区与地下水含水层间的导水系数；h_a,k为第k个积水区的差分单元处的地下水水位；为所述积水区的平均积水水位；

从所述未积水区渗出，并流入所述积水区的地下水流量为：

其中，为从所述未积水区渗出，并流入所述积水区的地下水流量；M为所述未积水区的地下水渗出的差分单元总数；为第m个未积水区的差分单元处地下水向采煤沉陷区的渗出补给量；C_s,m为第m个未积水区的差分单元处采煤沉陷区与地下水含水层间的导水系数；h_a,m为第m个积水区的差分单元处的地下水水位；Btm_m为第m个未积水区的差分单元处的底部高程；

所述积水区的水面蒸发量为：

其中，E为所述积水区的水面蒸发量；为所述积水区的平均积水面积； 为所述积水区单位面积内的水面蒸发强度；

所述积水区的渗漏量为：

其中，

为地下水水位高于所述采煤沉陷区的底部高程时所述积水区向所述地下水含水层的渗漏量；N为所述差分单元的总数；为第n个积水区的差分单元处采煤沉陷区的渗漏量；C_s,n为第n个积水区的差分单元处采煤沉陷区与地下水含水层间的导水系数；h_a,n为第n个积水区的差分单元处的地下水水位；为所述积水区的平均积水水位；

其中，

为地下水水位低于所述采煤沉陷区的底部高程时所述积水区向所述地下水含水层的渗漏量；S为所述差分单元的总数；为第s个积水区的差分单元处采煤沉陷区的渗漏量；C_s,s为第s个积水区的差分单元处采煤沉陷区与地下水含水层间的导水系数；Btm_s为第s个差分单元处的底部高程；为所述积水区的平均积水水位；

所述采煤沉陷区的人工取水量为：

W＝w_irri+w_indu+w_life+w_eco；

其中，W为所述采煤沉陷区的人工取水量；w_irri为采煤沉陷区的农业灌溉取水量；w_indu为所述采煤沉陷区的工业取水量；w_life为所述采煤沉陷区的生活取水量；w_eco为所述采煤沉陷区的生态环境取水量；

所述采煤沉陷区的排泄流量为：

其中，Q_so为所述采煤沉陷区的排泄流量；σ为堰流淹没系数；ε为堰流侧收缩系数；m为堰流流量系数；B₀为闸孔总净宽；H₀计入行近流速水头的堰上水深；h_s为下游淹没水深；g为重力加速度。