本发明涉及一种基于城市降雨径流污染负荷量模拟的河道黑臭预警方法，适用于对城市河道降雨后可能产生黑臭现象的预警，属于黑臭水体治理的。
背景技术：
1、随着点源污染得到有效治理，城市河道黑臭现象已得到有效遏制，然而城市的不透水面积不断扩大，降雨径流冲刷汇流所形成的面源污染已成为城市水环境质量下降的主要原因。
2、针对降雨后产生的面源污染通过模型模拟其污染负荷量，是准确预测场次降雨是否导致河道黑臭的关键。面源污染模型能够定量模拟污染物的迁移转化过程，现有的城市面源污染模型以子汇水区作为地表汇流模拟的计算单元，忽略子汇水区内部汇流过程，导致模拟结果产生较大偏差，进而无法准确预测降雨是否会导致河道黑臭。
技术实现思路
1、为了解决上述存在的问题，本发明公开了一种基于城市降雨径流污染负荷量模拟的河道黑臭预警方法，基于元胞自动机模拟城市降雨径流污染过程，从而改善现有城市面源污染模型中忽略子汇水区内部汇流过程所导致的较大偏差，实现城市降雨径流污染负荷准确核算的基础上对河道黑臭现象进行预警。其具体技术方案如下：
2、一种基于城市降雨径流污染负荷量模拟的河道黑臭预警方法，包括如下步骤：
3、步骤(1)：选定需要进行河道黑臭预警的城市小流域河段作为目标河段，搜集目标河段n场次黑臭水体对应的场次降雨产生的降雨径流数据，n为≥5的整数，降雨径流数据包括降雨径流中各污染物的浓度和总流量；根据管网数据及实地调研情况，确定目标河段流域内的雨水排口位置，于雨水排口处采集这n场次黑臭事件对应的场次降雨的降雨径流数据；
4、步骤(2)：对步骤(1)中n场次黑臭水体对应的场次降雨产生的降雨径流数据进行主成分分析，得到降雨径流中核心污染物，对筛选出来的核心污染物的n场次黑臭水体对应的场次降雨产生的降雨径流数据求其污染负荷平均值，记作即为发出河道黑臭现象预警的临界城市降雨径流核心污染物负荷指标；
5、步骤(3)：确定雨水排口收水范围，作为城市降雨径流污染的研究区，获取研究区的土地利用类型、数字高程数据(dem)和随时间变化的降雨量，据此设置与栅格数据相匹配的元胞空间，元胞空间包含元胞大小、元胞形状和边界元胞，设置基于元胞自动机的城市降雨径流污染负荷核算模型运行时间步长，建立基于元胞自动机的城市降雨径流污染负荷核算模型的基础地理场景；所述的元胞空间是研究区的地理空间，即为雨水排口的收水范围，由元胞组成；元胞形状为正方形网格，将研究区栅格化后，每个栅格对应一个元胞；利用雨量计测定或通过气象站点获得研究区的降雨量数据，使用arcgis软件对降雨数据进行空间插值处理，得到降雨数据的栅格表面，即获得每个元胞不同时刻的降雨量数据；
6、步骤(4)：根据质量守恒原理通过降雨量、植被截留降雨量、下渗量和洼地蓄积量计算每个元胞的净雨深，水流向下游多个方向流动，整体考虑中心元胞和邻域元胞的水位状态，采用多流向法中常用的最小差异法初步进行水量分配计算，使得水量分配之后的系统达到水位平衡的状态，在最小差异法计算出流量的基础上，采用曼宁公式对流量进行修正；
7、步骤(5)：采用指数函数冲刷模型计算各元胞内的非点源污染冲刷负荷，采用饱和函数累计模型计算模拟初始的干燥期地表污染物累积量，并综合考虑清扫管理措施以及径流冲刷，冲刷产生的污染物沿径流传输方向在各邻元胞之间传输交换，计算各元胞的非点源污染负荷净输入量，综合考虑因下渗而产生的非点源污染损失后，即可实现单个时间步长内各元胞的非点源污染负荷和浓度更新；
8、步骤(6)：采用一维元胞空间模拟城市雨水管网的运行情况，各元胞采用统一的坡度和断面尺度，整个管网的径流及非点源污染输入仅来源于地表中于管网入口连接的元胞，在较小的时间步长内将城市雨水管网中的雨水视作静态，故采用曼宁公式计算城市雨水管网中的径流和非点源污染的传输过程；
9、步骤(7)：当目标河段发生降雨时，通过步骤(3)至(6)对其产生的非点源污染负荷进行模拟计算，当核心污染物浓度达到发出河道黑臭现象预警的临界城市降雨径流核心污染物负荷指标立即发出河道黑臭预警。
10、进一步的，所述步骤(2)通过主成分分析得到降雨径流中核心污染物的具体过程为：
11、将各污染物浓度xij标准化，即对同一变量减去其均值再除以标准差，以消除量纲影响；
12、
13、式中，为标准化后的第i场次降雨的第j种污染物；为第j种污染物的样本均值；sj为第j种污染物的标准差；
14、在标准化数据阵x＝(xij)的基础上计算原始指标相关系数阵r＝(rij)，其中：
15、
16、式中，rij为第i种污染物与第j种污染物的相关系数；
17、求r阵的p个特征值λ1≥λ2≥…≥λp及相应的正则化单位特征向量cj＝(c1j，c2j，…cpj)′，其中每一特征值为对应主成分的方差，方差越大，则对总变异的贡献越大；
18、
19、式中yj为第j个主成分；是第k个标准化指标变量；
20、通过计算m个特征值的信息贡献率和累积贡献率合理选择主成分，计算式如下：
21、
22、式中，dj为第j个主成分yj的信息贡献率；aj为第j主成分yj的累计贡献率,通过设置阈值，确定主成分，即为降雨径流中核心污染物。
23、进一步的，所述步骤(4)中根据质量守恒原理将每个元胞的径流深度h通过降雨量r、植被截留降雨量p、下渗量i和蓄积量ds计算得到，计算公式如下：
24、h＝r-p-i-ds
25、式中，h为元胞的径流深度，mm；r为总降雨量，mm；p为植被截留雨量，mm；i为下渗量，mm；ds为蓄积量，mm；
26、下渗量i通过horton模型计算，计算公式如下：
27、
28、式中，i为下渗量，mm；fu为稳态或最小下渗速率，mm/s；f0为初始或最大下渗速率，mm/s；kd为下渗能力衰减系数，s-1；t为降雨持续时间，s-1；ts为时间步长；e为自然常数；
29、植被截留雨量p是指植被覆盖的地表在降雨时截留雨水的能力，计算式如下：
30、
31、pmax＝0.935+0.498lai-0.00575lai2
32、式中，p为植被截留雨量，mm；rcum为累积降雨量，mm；lai为叶面积指数；pmax为最大存储降雨量，该参数为lai的函数。
33、进一步的，所述步骤(4)中进行各方向元胞的分配水量qi的计算，具体为：计算中心元胞和其领域元胞水深的平均值，排除领域元胞中水深超过该平均值的元胞，不予分配水量，继续执行该过程直到没有领域元胞需要被排除，此时按照分配后剩余领域元胞水深相同的原则将中心元胞水量分配给剩余领域元胞中，计算公式如下：
34、
35、式中，qi为未修正的中心元胞向邻域元胞i转移的水量；hi为第i个邻域元胞的径流深度，mm；m为剩余邻域元胞数量。
36、进一步的，所述步骤(4)中对流量修正具体为：在产流和流向确定的基础上，对于部分元胞在一个时间步长内不能完成完整的水量分配的情况，采用曼宁公式对最小差异法进行时间上的修正，计算公式如下：
37、
38、式中，t为中心元胞水深转移至邻域元胞所需时间，既水量/污染物转移时间步长；v为流速，m/s；h为径流深度，m；j为中心元胞与第i个邻域元胞间的水面坡度；n为曼宁糙率系数；l为需要进行时间修正这个方向上的元胞间的距离；t为模型时间步长；qci表示修正后中心元胞向邻域元胞i转移的水量；
39、计算完每个邻域元胞转移水量qci求和可以得出转移总水量，将其与中心元胞水量进行对比，若中心元胞水量小于转移总水量，即中心元胞没有足够的水进行转移，每个邻域元胞水量转移量需按照水深进行缩小，重新计算qci；若中心元胞水深大于转移总水量，转移水量不变，将此规则应用于整个元胞空间，计算出每一时刻整个元胞空间内所有元胞的转移水量和剩余水量，并根据下一时刻步长的降雨和下渗量，获得下一时刻末的元胞水量，循环计算，直到整个模拟过程结束。
40、进一步的，所述步骤(5)中非点源污染冲刷负荷计算：采用指数函数冲刷模型公式计算每个元胞的非点源污染冲刷产生的污染负荷，计算公式如下：
41、
42、式中，w为冲刷产生的污染负荷，mg；kw为冲刷系数；nw为冲刷指数；mb为各元胞内各时刻的剩余污染物累积量，mg；
43、所述步骤(5)中地表污染物累积量计算：各元胞内各时刻的剩余污染物累积量按照饱和函数累计模型并综合考虑清扫管理措施以及径流冲刷计算，其计算公式如下：
44、mb＝bmaxad/(kb+ad)×a×(1-sse)-wt-1
45、式中：bmax为单位面积地表最大可能的污染物累积量，mg/mm2；kb为半饱和累计常数，即达到最大污染负荷累积量一半所需要的天数，d；ad为降雨前干燥天数，d；sse为清扫措施对累积污染物的去除率；wt-1为上一时刻径流冲刷去除的累积污染物负荷，mg；a为单个元胞面积，m2。
46、进一步的，所述步骤(5)中元胞的非点源污染负荷和浓度计算：各邻域元胞之间污染物沿径流方向进行传输交换，由此得到各元胞的非点源污染负荷净输入量，综合考虑因下渗而产生的非点源污染损失后即可实现单个时间步长内各元胞的非点源污染负荷和浓度更新，其计算公式如下：
47、l＝qci×ts×a×10-6×c+w
48、l(i,j)＝l(i+1,j+1)+l(i+1,j)+l(i+1,j-1)+l(i,j+1)+l(i,j-1)+l(i-1,j+1)+l(i-1,j)+l(i-1,j-1)+∑l(i,j)
49、式中：l代表各时间步长内各元胞沿径流方向向邻域元胞传输的非点源污染负荷，mg；c代表当前时刻各元胞内的非点源污染浓度，mg/l；l(i,j)代表位于坐标(i，j)的元胞内非点源污染负荷的净输入量，mg；l(i+1,j+1)+l(i+1,j)+l(i+1,j-1)+l(i,j+1)+l(i,j-1)+l(i-1,j+1)+l(i-1,j)+l(i-1,j-1)分别为来自各邻域元胞输入的非点源污染负荷，mg；∑l(i,j)代表(i,j)元胞内向各方向邻域元胞输出的非点源污染负荷总和，mg；
50、根据计算获得每个元胞当前时刻末的非点源污染负荷，循环计算，直到整个模拟过程结束。
51、进一步的，所述步骤(5)中综合考虑因下渗而产生的非点源污染损失后即可实现单个时间步长内各元胞的非点源污染负荷和浓度更新：
52、mt+1＝qrcr+l(i,j)-qic
53、
54、qi＝(i0+ftc)×a/106
55、式中，mt+1代表各元胞内更新的地表非点源污染物负荷，mg；c代表当前时刻各元胞内的非点源污染浓度，mg/l；qr代表各元胞输入的降雨量，l；cr为降雨中含有的污染物浓度，mg/l；qi代表下渗径流量，l。
56、进一步的，所述步骤(6)城市雨水管网中径流传输过程的计算：在较小的时间步长内可以将城市雨水管网中的雨水视作静态，故采用曼宁公式计算城市雨水管网中的径流传输过程，计算公式如下：
57、
58、式中：qp代表中心元胞向下游元胞输出的径流量，mm3/s；hp代表中心元胞径流深，mm；sp代表管道坡度；np代表管道曼宁系数；dp代表断面尺寸，mm；
59、城市雨水管网中非点源污染传输过程的计算：在较小的时间步长内可以将城市雨水管网中的雨水视作静态，故采用曼宁公式计算城市雨水管网中的非点源污染的传输过程，计算公式如下：
60、lp＝qp×ts×10-6×cp
61、lp代表中心元胞向下游输出的非点源污染负荷，mg；cp代表各管网元胞中的非点源污染浓度，mg/l。
62、进一步的，所述步骤(6)城市雨水管网中元胞状态的全局模拟：每个时间步长内的管道元胞内的径流和非点源污染负荷仅受两个邻域元胞和上一个时间步长自身的元胞状态决定，计算各时刻管道元胞的径流和非点源污染负荷的公式如下：
63、
64、
65、式中：qp_t+1和qp_t分别代表当前时刻和下一时刻管道元胞内的径流量，mm3；pp_t和pp_t+1分别代表当前时刻和下一时刻管道元胞内的非点源污染负荷，mg；k为一阶衰减系数，d-1；和分别为各元胞的径流入量和出流量，mm3；和分别代表各元胞内的非点源污染负荷输入和输出量，mg；
66、得出整个元胞空间每个元胞当前时刻末的非点源污染负荷，循环计算，直到整个模拟过程结束。本发明的有益效果是：
67、与基于子汇水区的城市面源污染模型相比，元胞自动机具备与输入栅格数据分辨率相匹配的网格结构，可以完全表现城市地表特征，因此元胞自动机算法具有较好的城市地表降雨径流过程模拟能力。
68、与现有技术相比，本发明在核算降雨产生的非点源污染负荷时考虑了城市降雨径流地表子汇水区内部汇流过程。该模型考虑了地表子汇水区内部汇流过程对城市降雨径流的影响，使用了真实的时间步长，相较于传统模型，本发明提出的模型计算单元不是仅具有单一地理特征的子汇水区，不仅能够准确模拟城市地表降雨径流过程，而且能够模拟不同降雨强度下径流污染分布规律和研究区的降雨径流污染负荷，进而提高了预警的准确性。
1.一种基于城市降雨径流污染负荷量模拟的河道黑臭预警方法，其特征在于，包括如下步骤：
2.根据权利要求1所述的基于城市降雨径流污染负荷量模拟的河道黑臭预警方法，其特征在于，所述步骤(2)通过主成分分析得到降雨径流中核心污染物的具体过程为：
3.根据权利要求1所述的基于城市降雨径流污染负荷量模拟的河道黑臭预警方法，其特征在于，所述步骤(4)中根据质量守恒原理将每个元胞的径流深度h通过降雨量r、植被截留降雨量p、下渗量i和蓄积量ds计算得到，计算公式如下：
4.根据权利要求1所述的基于城市降雨径流污染负荷量模拟的河道黑臭预警方法，其特征在于，所述步骤(4)中进行各方向元胞的分配水量qi的计算，具体为：计算中心元胞和其领域元胞水深的平均值，排除领域元胞中水深超过该平均值的元胞，不予分配水量，继续执行该过程直到没有领域元胞需要被排除，此时按照分配后剩余领域元胞水深相同的原则将中心元胞水量分配给剩余领域元胞中，计算公式如下：
5.根据权利要求1所述的基于城市降雨径流污染负荷量模拟的河道黑臭预警方法，其特征在于，所述步骤(4)中对流量修正具体为：在产流和流向确定的基础上，对于部分元胞在一个时间步长内不能完成完整的水量分配的情况，采用曼宁公式对最小差异法进行时间上的修正，计算公式如下：
6.根据权利要求1所述的基于城市降雨径流污染负荷量模拟的河道黑臭预警方法，其特征在于，所述步骤(5)中非点源污染冲刷负荷计算：采用指数函数冲刷模型公式计算每个元胞的非点源污染冲刷产生的污染负荷，计算公式如下：
7.根据权利要求1所述的基于城市降雨径流污染负荷量模拟的河道黑臭预警方法，其特征在于，所述步骤(5)中元胞的非点源污染负荷和浓度计算：各邻域元胞之间污染物沿径流方向进行传输交换，由此得到各元胞的非点源污染负荷净输入量，综合考虑因下渗而产生的非点源污染损失后即可实现单个时间步长内各元胞的非点源污染负荷和浓度更新，其计算公式如下：
8.根据权利要求1所述的基于城市降雨径流污染负荷量模拟的河道黑臭预警方法，其特征在于，所述步骤(5)中综合考虑因下渗而产生的非点源污染损失后即可实现单个时间步长内各元胞的非点源污染负荷和浓度更新：
9.根据权利要求1所述的基于城市降雨径流污染负荷量模拟的河道黑臭预警方法，其特征在于，所述步骤(6)城市雨水管网中径流传输过程的计算：在较小的时间步长内可以将城市雨水管网中的雨水视作静态，故采用曼宁公式计算城市雨水管网中的径流传输过程，计算公式如下：
10.根据权利要求1所述的基于城市降雨径流污染负荷量模拟的河道黑臭预警方法，其特征在于，所述步骤(6)城市雨水管网中元胞状态的全局模拟：每个时间步长内的管道元胞内的径流和非点源污染负荷仅受两个邻域元胞和上一个时间步长自身的元胞状态决定，计算各时刻管道元胞的径流和非点源污染负荷的公式如下：
