cumcm0508

摘要:

本文建立了基于O'Conner修正式的一维动力学模型，对长江水域的水质评价、预测和控制问题进行讨论与分析。

通过对大量数据的统计分析，我们给出了长江两年多来水质污染状况的定量综合评价。

为了预测未来十年长江的污染情况，我们构造了动态的马尔可夫状态转移矩阵以描述三大类水质（可饮用水，中度污染水，劣质水）的内在联系以及与排污总量之间的联系，并结合总排污量的指数增长模型对未来污染发展趋势做出预测。

为了从纷繁芜杂的数据中寻找规律，论文中采用了具有反馈的Elman人工神经网络，同时使用“时间窗滑动法”对长江水质未来的发展趋势给出了评价与预测。

在论文拓展中，我们讨论了电路模拟模型和四房室模型的建模思想，并对题目中的问题进行了描述。遗憾的是没有进行相关的求解。

论文的主要求解结果如下：

长江在近两年内，第二类水质所占总量比例为61%左右。各地区水质也大多属于第二类。题目中的两种污染物的主要污染源均位于湖北宜昌南津关和湖南岳阳城陵矶之间。全流域内年均排污总量是呈逐年递增的；可饮用水比例逐年降低，水域污染日益严重。若要保证可饮用水比例不低于80%需对大量污水进行处理，2005年为112亿吨，2014年为325.66亿吨。

关键字：

水质评价预测 动态的马尔可夫预测模型 一维动力学模型 nj阶时间窗滑动法 电路模拟模型 Elman神经网络模型 房室模型

1

一．问题重述：

长江水质的评价和预测

水是人类赖以生存的资源，保护水资源就是保护我们自己，对于我国大江大河水资源的保护和治理应是重中之重。专家们呼吁：“以人为本，建设文明和谐社会，改善人与自然的环境，减少污染。”

长江是我国第一、世界第三大河流，长江水质的污染程度日趋严重，已引起了相关部门和专家们的高度重视。2004年10月，由全国政协与中国发展研究院联合组成“保护长江万里行”考察团，从长江上游宜宾到下游上海，对沿线21个重点城市做了实地考察，揭示了一幅长江污染的真实画面，其污染程度让人触目惊心。为此，专家们提出“若不及时拯救，长江生态10年内将濒临崩溃”，并发出了“拿什么拯救癌变长江”的呼唤。

附件3给出了长江沿线17个观测站（地区）近两年多主要水质指标的检测数据，以及干流上７个观测站近一年多的基本数据（站点距离、水流量和水流速）。通常认为一个观测站（地区）的水质污染主要来自于本地区的排污和上游的污水。一般说来，江河自身对污染物都有一定的自然净化能力，即污染物在水环境中通过物理降解、化学降解和生物降解等使水中污染物的浓度降低。反映江河自然净化能力的指标称为降解系数。事实上，长江干流的自然净化能力可以认为是近似均匀的，根据检测可知，主要污染物高锰酸盐指数和氨氮的降解系数通常介于0.1~0.5之间，比如可以考虑取0.2 (单位：1/天)。附件4是“1995~2004年长江流域水质报告”给出的主要统计数据。下面的附表是国标(GB3838-2002) 给出的《地表水环境质量标准》中4个主要项目标准限值，其中Ⅰ、Ⅱ、Ⅲ类为可饮用水。 请你们研究下列问题：

（1）对长江近两年多的水质情况做出定量的综合评价，并分析各地区水质的污染状况。 （2）研究、分析长江干流近一年多主要污染物高锰酸盐指数和氨氮的污染源主要在哪

些地区?

（3）假如不采取更有效的治理措施，依照过去10年的主要统计数据，对长江未来水质

污染的发展趋势做出预测分析，比如研究未来10年的情况。

（4）根据你的预测分析，如果未来10年内每年都要求长江干流的Ⅳ类和Ⅴ类水的比例

控制在20%以内，且没有劣Ⅴ类水,那么每年需要处理多少污水？ （5）你对解决长江水质污染问题有什么切实可行的建议和意见。

二．问题分析：

污染物的移流扩散大致分为三个阶段：第一阶段为初始稀释阶段。该阶段主要发生在污染源附近区域（污染带），运动主要为沿水深的垂向浓度逐渐均匀化。第二阶段为污染扩散阶段。该阶段中，由于存在浓度梯度，污染由垂向均匀化向过水断面均匀化发展。第三阶段为纵向离散阶段。该阶段中，由于沿水流方向的浓度梯度作用以及断面上

2

流速分布，出现了污染沿纵向的移流扩散。

具体分析题目及数据，我们认为由于本题考虑的是整个长江流域，河域范围十分广泛，故我们认为无需考虑污染物扩散的前两个阶段。又由题意可知评价水质状况主要有溶解氧、高锰酸盐指数、氨氮、PH值四个指标，经分析我们认为这四个指标主要受到污染源、江河自身对污染物的自然净化能力及人为因素的影响，故我们应该从各因素对指标的影响的角度去解决问题。

三．一维动力学水质模型的建立

¾ 模型假设：

1．对于水体中污染物变化规律的假设：

溶解氧、高锰酸盐指数、氨氮、PH值四个指标相互，互不影响。 2．对于污染源排污方式的假设：

[1]污染源皆为点源，即在集中点排放工业废水和城镇生活污水等污染物的污染源。 [2]污染源均匀连续排放污染物。

[3]污染源排放的污染物瞬间在该断面均匀分布，即不考虑污染物在河流垂向上的扩散过程。

[4]单位距离(1公里)内污染物的浓度不变。

¾ 建模准备： 1．有关名词说明：

周期：从所给数据的最早月份开始，每十二个月为一个周期。对于本题而言,对周期

做如下标号：

周期 2003.6至2004.52004.6至2005.52005.6至2005.9

标号1 2 3

水期：每一个周期分为枯水期、丰水期、平水期三个水期，分别标号为水期1,2,3。 2．观测点标号说明：

标号 观测点名称 标号观测点名称 1 四川攀枝花 10 四川泸州沱江二桥2 重庆朱沱 11 湖北丹江口胡家岭3 湖北宜昌南津关 12 湖南长沙新港 4 湖南岳阳城陵矶 13 湖南岳阳岳阳楼 5 江西九江河西水厂14 湖北武汉宗关 6 安徽安庆皖河口 15 江西南昌滁槎 7 江苏南京林山 16 江西九江石 8 四川乐山岷江大桥17 江苏扬州三江营 9 四川宜宾凉姜沟

3．水质类别标号说明：

Ⅰ类、Ⅱ类、Ⅲ类、Ⅳ类、Ⅴ类、劣Ⅴ类分别对应类别标号1至6。 4．坐标系的建立：

查阅资料得到题目中所给的各个观测点的地理位置如下图所示：

3

长江流域水系简图

为了便于对问题的描述，我们建立以长江源头为原点沿干流至入海口的一维坐标系。

¾ 符号约定： D ：水深

E ：各种污染物及溶解氧的综合扩散系数u ：水体中有机物的沉降速率Cs ：水体中饱和溶解氧浓度RDO ：大气溶解氧系数SDO ：底泥沉积物耗氧速率M ：周期总数I ：观测点总数J ：污染物种类总数m ： 周期标号(m=1,...,M)n ：水期标号(n=1,2,3)l ：月份标号(l=1,...,12)i ：观测点序号(i=1,...I)j ：污染物种类标号(j=1,...,J)k ：水质类别标号(k=1,...6)xi ：第i个观测点的坐标Vi ：第i个观测点的水流量dj ：水体对第j种污染物的降解系数

λj ：第j种污染物溶解量占总量的比例

wij ：第i个观测点单位时间内第j种污染物的排放量

4

v(x,t) ：时间t位置x处的水流速CDO(x,t) ：时间t位置x处的溶解氧浓度Cj(x,t) ：时间t位置x处第j种污染物的浓度

cilj ：第i个观察点在第l个月测得的第j种污染指标浓度Ωkj(t) ：时间t以第j种污染物为指标时对应于第k类别的浓度区间Lkj(t) ：时间t以第j种污染物为指标时第k类水质的总河长Mi ：第i个水文年污水的总量

∆Mi ：第i个水文年相对于第i-1个水文年污水总量的变化量

¾ 建模思想：

如下图所示，我们在模型中综合考虑以下几个方面：

(1) 分析研究河段内污染物总量的变化：

总量增大的原因：该段污染源、上游干流及支流污染源的注入。 总量减小的原因：研究河段内污染物浓度减小的原因主要是在水环境中的物理降解、

化学降解和生物降解等自然降解，可综合考虑有机物在生物降解的过程中对水体溶解氧的消耗。

(2)可能污染源的讨论：

经过查阅资料，我们发现与干流的水流量相比，出口支流的水流量较小，故可认为出口之流对干流污染物浓度的变化没有影响，因此我们不考虑出口支流作为污染源的情况，即认为污染源只有可能位于干流上任一点及入口支流处。

¾ 模型建立：

经过上述分析，我们建立了基于O’Conner修正式的长江水质的评价、预测和控制一维动力学模型：

5

J

∂2CDO(x,t)∂CDO(x,t)∂CDO(x,t)SDO

1)+v⋅=E⋅−R⋅(C−C)−d⋅C(x,t)−∑DOsDOjj

∂t∂xD∂x2j=1

2)

∂Cj(x,t)∂t

Jj=1

+v⋅

∂Cj(x,t)∂x

=E⋅

∂2Cj(x,t)∂x2

+zj(x,t)−dj⋅Cj(x,t)−

u

⋅(1−λj)⋅Cj(x,t)(j=1,...,J)D

3)M=∑∫∫ρj⋅zj(x,t)dxdt

4)Lkj(t)=L[Cj(x,t)] Cj(x,t)∈Ωkj (j=1,...,J;k=1,...,K)5)ηkj(t)=Lkj(t)

J

∑L

a=1

kj

(t)×100% (j=1,...,J;k=1,...,K)

(1)模型中各方程的说明：

1)式描述研究区域内溶解氧的浓度变化规律。

2)式描述研究区域内第j种污染物的浓度变化规律。 3)式描述全流域各种污染物的总量。

4)式描述各个水质类别的河长的抽象函数。

5)式描述各个水质类别的河长在总河长中所占的比重。 (2)不同水质类别河长Lkj(t)的计算方法：

在求解各个水质类别的河长时，我们使用“分类一维搜索”的方法：按照所给的第j种污染物的分类标准，从源头开始逐点搜索，记录满足Cj(x,t)∈Ωkj的第一个坐标点，继续搜索直到找到属于另一水质类别的第一个坐标点，记录该点坐标值，则两个坐标值之差即为该类别河水河长的一部分，对整个流域进行遍历搜索，将符合该类别浓度要求的河长进行累加，最终可得到该类别水的总河长，以同样方法搜索可得其他各类水质的河长。

(3)单位时间排污口污染物浓度的变化zj(x,t) 的构造：

由于污染源和支流注入的水具有相似的特征，我们构造zj(x,t)来表示单位时间

排污口（包括污染源和入口支流）的污染物浓度的变化，这样做一方面可以方便求解每个研究区域内的污染物浓度，另一方面可以有效地给出全流域内污染物的总量。 3．模型的化简： (1)Ex⋅

∂2Cj(x,t)∂x2

∂2CDO(x,t)，Ex⋅ ：这两项描述河流的弥散作用对污染物、溶解氧浓2

∂x

度的影响，由于在一般的情况下河流的弥散作用对污染物、溶解氧浓度的影响远小于河水的纵向迁移对污染物、溶解氧浓度的影响，因此可以忽略不考虑。 u

(2)⋅(1−αj)⋅Cj(x,t) ：该项描述河流对非溶解有机物的沉降作用，对于本题而言，由

D

于所给的数据没有包含此方面内容，所以可以忽略此项。 (3)

SDO

：描述底泥沉积物的耗氧，同上在本题中也不予考虑。 D

6

4．模型在本题的应用 (1)问题一的解决：

由模型中的(5)式可以直接得到整个长江流域各个水质类别的河长在总河长中所占的比重，进而可对长江的水质进行综合评价，若想通过属于各个类别的断面个数比例进行评价，只需应用模型求得各个断面所属类别进行统计即可；由模型亦可得到各地区的污染物浓度，由之可对各个地区的水质污染状况进行分析。 (2)问题二的解决：

由模型求解得到各个分段地区的各类污染物浓度之后，即可以在全流域内进行主要污染源的寻找，认为污染物浓度大的地方可能存在污染源。 (3)问题三、四的解决：

对于问题三、四，模型中的(3)式可以直接给出全流域总的污染量，(5)式可以直接给出各个水质类别的百分比。这样以全流域各个点的各类污染物浓度作为连接，可以直接给出各个类别河水的百分比和总污染量，这样就可以得到问题3，4的处理方法。

四．长江水质的综合评价：

¾ 长江近两年多水质情况的综合评价：

由于题目中给出了长江17个观测点的相应数据，首先我们认为如果可以得到各观测点在各个周期的水质类别，则可以用各周期长江分别属于六个类别的观测点百分比对长江水质情况进行综合评价。解决问题分如下三个步骤： 步骤一：评价各周期各个观测点在各个水期的水质情况。

数据给出的为各个月份17个观测点的四个指标数值，由四个指标值分别可以对应得到一个水质等级，在确定各观测点在该月份的水质类别时，我们认为水质类别反映的是水质最差的方面，故应该以上面得到的四个类别中最差的一个作为各观测点在该月份的水质类别。而在确定各观测点在各个水期的水质类别时，考虑到观测点在每一水期的各个月份的水流量的不同，我们用一个水期内水流量对浓度的加权平均作为该水期的污染物浓度：

l∈水期n

∑c

ijl

⋅Vl

l∈水期n

∑V

l

(i=1...17,j=1...4)

进而又取四个类别中最差的一个作 为各观测点在该月份的水质类别。经如上分析我们可以得到17个观测点在各周期各个水期的水质类别。

步骤二：评价各个观测点各个周期的水质情况。

表一给出了17个观测点在各周期各个水期的水质类别，采用与上一步骤相同的处理方法，我们以17个观测点在各周期各个水期的水质类别中最差的一个作为各观测点在各周期的水质类别。 即：

观测点i第m个周期的水质等级= max（观测点i第m周期各个水期的水质类别） 得到如下表二：

时间＼观测点

2003年6月至2004年5月2004年6月至2005年5月

122

223

322

422

522

622

722

8 9 54

42

7

2005年6月至2005年9月2211221

212443

313232

214663

2156

216342

22

时间＼观测点 102003年6月至2004年5月2004年6月至2005年5月2005年6月至2005年9月

622

17 2 3 2

表二

步骤三：统计各个周期中属于各个类别的观测点总数。

周期＼类别

2003年6月至2004年5月 2004年6月至2005年5月 2005年6月至2005年9月

Ⅰ类 Ⅱ类 Ⅲ类 Ⅳ类Ⅴ类 劣Ⅴ类

0 10 1 2 1 3 0 9 3 3 0 2 1 12 3 1 0 0 61.66%

表三

总计 1.67%11.67%10%1.67% 8.33% 由于2005年6月至2005年9月仅仅为丰水期对应的月份，因此我们仅对2003

年6月至2004年5月与2004年6月至2005年5月进行比较讨论，由上表可见：2004年6月至2005年5月中劣于Ⅲ类的观测点数较2003年6月至2004年5月多，说明长江的水质变差，同时可见整个长江水质属于类别Ⅱ的观测点个数最多，占总数的61.66%，故可认为长江的整体水质为Ⅱ类。 ¾ 各地区水质污染状况分析：

为了分析各地区水质污染状况，我们对所给数据进行统计，计算各个地区在所给数据时间范围内水质状况分别属于各个类别的时间比例，分别考虑干流观测点和支流观测点，得到如下表格：

类别＼观测点 1 2 3 4 5 6 7 Ⅰ类 21.43% 0.00%Ⅱ类 Ⅳ类 Ⅴ类

78.57% 85.71%0.00% 0.00%0.00% 0.00%Ⅲ类 0.00% 14.29%

0.00% 0.00%100.00%

85.71%

0.00% 14.29%0.00% 0.00%0.00% 0.00%0.00% 0.00%

0.00% 7.14% 7.14% 100.00%

92.86% 92.86% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 0.00% 劣Ⅴ类 0.00% 0.00%表四：干流观测点水质评价结果表

由上表可见，干流中观测点1、3、5、6、7的水质较好，水质一直为Ⅰ类或Ⅱ类；观测点2、4的水质次之，水质属于Ⅱ类的时间比例为85.71%，其余时间属于Ⅲ类。

类别＼观测点 8 9 10 11 12 13 14 15 16 17 Ⅰ类 0.00% Ⅱ类 Ⅲ类

0.00%

0.00%

.29%

0.00%0.00%

0.00% 57.14%42.86%0.00% 0.00% 0.00% 100.00%0.00%

0.00% 78.57%21.43%0.00% 0.00% 0.00% 100.00%0.00%

0.00% 0.00% 0.00%

0.00%

0.00%

21.43% 78.57% 57.14% 35.71%71.43% 78.57%21.43% 21.43%21.43% 0.00% 28.57% 0.00% 71.43%0.00% 0.00%

28.57%0.00%0.00%71.43%28.57%

Ⅳ类 42.86% 21.43% 0.00% Ⅴ类 14.28% 0.00% 劣Ⅴ类 0.00%

0.00%

0.00%

14.29% 7.14% 0.00% 0.00%

0.00%

0.00% 0.00%

14.29% 0.00% 85.71% 0.00% 0.00%

不低于Ⅲ类 42.86% 78.57% 85.71% 100.00%劣于Ⅲ类

57.14% 21.43% 14.29% 0.00%

92.86% 100.00%100.00% 7.14% 0.00%

表五：支流观测点水质评价结果表

8

由上表可见，支流中观测点11、13、14、17的水质较好，水质一直不低于Ⅲ类；观测点9、10、12、16的水质较好，水质不低于Ⅲ类的时间比例大于90%；观测点8的水质较差，水质不低于Ⅲ类的时间比例为42.86%；观测点15的水质最差，水质一直劣于Ⅲ类。

综合考虑上述二表亦可知，长江干流的水质整体上较支流水质好的多。

五．长江干流主要污染源的确定：

应用已经建立的一维动力学模型，我们完全可以对该问题进行处理，但是在实际的求解过程中，为了方便计算，我们提出了相应的假设对上述模型进行适当的简化。 ¾ 模型求解的简化：

由于题目给定的支流和相关全流域内存在的污染源的数据太少，并且通过查阅相关资料，我们发现题目中所给的支流并未包括了长江的所有主要支流，且支流的流量相对于干流流量较小（如下表六所示），故在模型求解中，我们对该问题进行了适当的简化，认为污染物排放口只可能是干流上的7 个观测点，故只需对这7个观测点进行分析，从中得到主要污染源。这样问题简化为：在污染源已定的情况下求各个污染源的排污量，并由排污量给出干流上各个观测点所属类别，划分各污染源的等级。

约定：等级越大，污染能力越强。

支流

雅砻江 大渡河 岷 江

1570

2850

嘉陵江2120

乌 江1650

沅 江2170

湘 江 汉 江 赣 江2370

1710

2130

年径流量 1810

表六：长江主要支流的多年平均水流量

¾ 建模准备：污染物剩余比例函数的构造

步骤一：水流由初始位置x0到达任意点x所用时间的计算推导。

根据物体运动规律可得：

t

∫

t0

v[x(t)]dt=x(t)−x0

求解上式得到水流时间计算公式：

t−t0=∫

1

dx x0v(x)

x

其中：t0,x0分别为水流的初始时间和初始位置。 公式的物理意义：

dx

表示水流流速为v(x)时通过dx的路程所需时间，则对通过的路程进行积分v(x)即得到水流通过该路程所需的时间。

步骤二：构造污染物剩余比例函数。

9

由步骤一可以得到xq处的污染物达到xi所需时间为：

t=∫

iq

dx

xqv(x)

xi

又由题意可知，长江干流的自然净化能力可以认为是近似均匀的，即污染物每天被降解的程度相同，所以污染物在水环境中的时间越长，剩余的质量越少，因此我们构造如下的污染物剩余比例表示自然降解对于污染物质量的减少程度：

r=(1−K

iq

i

(j)tq

)

¾ 简化模型的建立：

基于以上分析准备，我们建立如下的简化模型：

∑wqj⋅rqi=Cj(xi,t)⋅Vi (i=1,...,7)

q=17

说明：各个方程描述各个干流观测点污染物总量等于上游各观测点的注入与本地注入的污染物质量之和。

¾ 简化模型求解：

在求解的过程中，首先得到各个月份各观测点水流速的线形插值，并利用Matlab求解最终得到两种主要污染物高锰酸盐和氨氮在各干流观测点单位时间内的注入量，如下表所示：

干流 1 2 3 4 5 6 7

高锰酸盐（kg/S）

单位时间流经总量 单位时间本地注出量

8.7729 29.971 49.732 68.959 58.113 59.574 57.048

8.7729 25.961 38.990 52.452 38.151 19.259 33.210

0.67494 3.321 3.6518 6.5146 4.6824 4.9708 2.840

氨氮(kg/S)

单位时间流经总量

单位时间本地注出量

0.67494 3.0126 2.4615 5.3025 2.7966 1.7225 0.85098

分析上表数据：

对于整个长江流域而言，评价污染源等级的指标是各污染源单位时间的本地注入量，而由于水体自然降解作用的影响，各污染源对其他观测点污染物总量的影响不同，具体而言即某污染源的污染物流至其他观测点时可能会由于水体的降解作用而对该观测点的污染物总量没有影响，也就是说各个观测点有其相应的主要污染源。 分别研究各观测点两种污染物的情况得条形图如下所示：

10

对整个长江流域的污染情况进行研究，分析各观测点单位时间内高锰酸盐的本地注入量，对各观测点进行分类：观测点4即湖南岳阳城陵矶为高锰酸盐重污染源，观测点2、3、5、7为二级污染源，观测点1、6的污染能力最低。

对整个长江流域的污染情况进行研究，分析各观测点单位时间内氨氮的本地注入

观测点2、3、量，对各观测点进行分类：观测点4即湖南岳阳城陵矶为氨氮重污染源，

5、6为二级污染源，观测点1、7的污染能力最低。 问题二结果的综合分析：

从单位时间内本地的注入量和单位时间内的流经总量，比较可以得到:在某些地

11

区，单位时间的流经总量大，但是单位时间内本地的注入量不大，由此可得:对于该地区而言，它的污染主要来自上游，而自身并不是一个主要污染源。例如本题中的观测点6。

另外，为了方便求解，在计算当中我们考虑只有观测点处才有污染物的注入，如果考虑实际情况，我们可以通过已求得的解来大概推测污染源的所处范围。例如：本题中，4 处的流经污染量较大，我们就可以认为在3 与4 之间，4与5之间都有可能存在着主要的污染源。当然，此结论是在假设污染源的单位时间注入量符合实际的情况下得到的。

六．长江水质预测与污水控制：

分析问题三和问题四，可知我们要做的即是给出长江未来水质发展的趋势和相应处理污水的方案。在解决过程中，我们以一个水文年为单位时间，对长江未来水质进行预测和控制，并且按照用途对六个类别的水进行合并，得到适合于分析和预测的三个水质类别：可饮用水A类（包括原Ⅰ类、Ⅱ类和Ⅲ类水）、中度污水B类（包括原Ⅳ类、Ⅴ类水）和劣质水C类（原劣Ⅴ类）。同时我们认为每个水文年的总污水量和各个类别水资源所占的比例之间存在着一定的函数关系，而不是存在的，所以在模型建立和求解中，需对二者统一考虑。通过分析，我们建立了以下两个模型：

¾ 动态马尔可夫预测模型的建立： 1．建模准备：

如图所示，我们认为A类、B类和C类水质可看作是在一个物质交换系统中的三个类，每个类中的元素均有机会互相转换，并且在不考虑河流中的污染物积累时，可近似认为这种相互转换具有马尔可夫性，因此我们就可由这种经过简化的结构来推测马尔可夫转移矩阵的结构。

然而，两个具有不同排污总量的年份所具有的三类水质的结构应该是不同的，两个具有相同排污总量的年份所具有的三种水质的比例却应该大致相同，这说明这个矩阵不应简单地构造成静态的转移矩阵，而应该是随时间或者说是随排污总量的不同而不同的。

如上图，我们先假设三类水之间可以相互转换，且先不考虑归一化的问题。

hAA

设A,B,C三类水之间的转换流矩阵P0=hBA

hCA

hABhBBhCB

hAC

hBC，其中hAA, hBB, hCC不hCC

随时间变化，它们只与水体自身的性质有关，物理含义为表征一种物质留在原有状

12

hAB表示A类水流向B类水的流量，hBA表示B类水流向A类水的流量。态的能力。hAC

与hCA表示可饮用水类别与劣质水类别之间的交换，由于我们考虑的是一个大的系统，可饮用水直接变成劣质水的概率很小，同理劣质水直接变成可饮用水的概率也很小，所以可作近似假设认为hAC=hCA=0，这样我们需求hAB, hBA, hCB, hBC四个动态量。

2．马氏转移矩阵的构造：

我们认为hAB, hBA, hCB, hBC四个动态变量是∆Mi的函数，由此我们构造了以下函数形式：

∆Mi∆Mi∆Mi∆Mi

=,=,=,=ααααhhhhAB1BA2BC3CB4

MiMMMiii



∆Mi∆Mi∆Mi∆Mi

=−,=−,=−,=−hαhαhαhα1AB2CB3BC4BA

MMMMiiii

∆Mi>0

∆Mi<0

考虑∆Mi的正负，我们得到如下两个矩阵：





∆Mi>0

∆Mi<0

Pi(1)

Pi(2)

∆Mi

hAA0α1

Mi

∆Mi∆Mi

hBBα3=α2

MMii

∆Mi

h0αCC4Mi

∆MihAA0−α2

Mi

∆Mi∆Mi

hBB=−α1−α4

MiMi



∆Mi

0αhCC−3Mi





注：∆Mi<0表示假设的流向刚好与实际流向相反。 对两个矩阵进行归一化处理如下：

13

∆Mi>0时

hAA⋅Miα1⋅∆Mi0

hMαMhMαM⋅+⋅∆⋅+⋅∆11AAiiAAii

αMhMαM⋅∆⋅⋅∆23iBBii~(1)P=α2⋅∆Mi+hBB⋅Mi+α3⋅∆Miα2⋅∆Mi+hBB⋅Mi+α3⋅∆Miα2⋅∆Mi+hBB⋅Mi+α3⋅∆Mi

i



αMhM⋅∆⋅4iCCi

0α4⋅∆Mi+hCC⋅Miα4⋅∆Mi+hCC⋅Mi



∆Mi<0时hAA⋅Mi−α2⋅∆Mi0

hMαMhMαM⋅−⋅∆⋅−⋅∆AAi2iAAi2i



αMhMαM−⋅∆⋅−⋅∆1iBBi4i~

(2)

Pi=−α1⋅∆Mi+hBB⋅Mi−α4⋅∆Mi−α1⋅∆Mi+hBB⋅Mi−α4⋅∆Mi−α1⋅∆Mi+hBB⋅Mi−α4⋅∆Mi

αMhM−⋅∆⋅3iCCi

0−α3⋅∆Mi+hCC⋅Mi−α3⋅∆Mi+hCC⋅Mi



~~

则Pi(1)与Pi(2)就是所需的动态马尔可夫矩阵。用µi={ai , bi , ci}表示第i各水文年

的水质向量（ai , bi , ci分别表示第i个水文年A、B、C三类水的百分比），则有：

µi+1

~

µiPi(1) ∆Mi≥0

=~(2)

µP M0 ∆3．马氏转移矩阵中参数的确定：

用题目中提供的十年的数据构造最小二乘误差函数来确定模型中的参数，当然也要注意转移概率作为概率非负且不大于1的约束。通过观察易知各个参数若满足比例就是等价的，在此不赘述。 4．各水文年废水总量的确定：

考虑不受阻滞的经济增长，我们认为表征沿岸经济发展规模的因子Q(t)符合指数增长规律,因此我们建立如下的模型：

Q(t)=a⋅eb⋅t

同时考虑废水排放的总量是与长江沿岸经济的增长是呈线性的关系的：

M(t)=γ⋅Q(t)

记c=γ⋅a则有： M(t)=c⋅eb⋅t

14

由此去拟和前十个水文年的废水总量，确定其参数，在去预测未来十年的废水排放的总量

指数模型拟和污水排放总量600函数拟和与预测55050045040035030025020015010095% 致信区间年污水排放量实际数据年污水排放量19951996199719981999200020012002200320042005200620072008200920102011201220132014年份 由图可以看出，拟和结果与实际符合得相当好的，回归检验结果R2=0.9583, 系数c=156置信区间为(144.7, 167.2)，参数b=0.05965置信区间为(0.04939, 0.06991)。图中虚线给出置信度为0.95的置信区间估计。

5．模型求解：

(1)问题三的求解：

步骤一:根据已经拟和好的指数增长模型，可以求得未来十年污水排放总量Mi

~~

步骤二:根据Mi和∆Mi构造Pi(1)与Pi(2)

步骤三:用迭代公式µi+1最终求解结果如图：

~

µiPi(1) ∆Mi≥0

求得µi+1 =~(2)

µiPi ∆Mi<0

15

2014年长江干流三种水质水域长度的百分比分别是：46.69% 25.77% 27.54%

2014年长江干流三种水质水域长度的百分比分别是：46.28% 25.69% 28.03%。

年份 废水总量（亿吨）年份废水总量（亿吨） 2005 300.67 2010405.15 2006 319.15 2011430.05 2007 338.76 2012456.49 2008 359.59 2013484.54 2009 381.69 2014514.33 2014年的排污总量已经达到514.33亿吨。可见按照现有的发展趋势，不超过十年，长江全流域包括干流就会有超过一半流域的水质达不到饮用水标准，近三分之一的流域被严重污染。

(2)问题四的求解：

~~

由Pi(1)与Pi(2)的结构易知，实际上由于矩阵第三列不全为零，我们无法做出关于劣V类在一年内被净化为0的事件的预测，该事件只可能经过无穷长的时间才有可能发生。于是我们不是一般性地假设劣V类小于3%时就达到了要求。而且，我们假设污水被处理就意味着不曾也将不会污染长江。

&使用我们在求解时运用了这样的求解方法：我们从高到低搜索满足条件的Mi

~(1)~(2)

&=M&−M&构造的P&是实际未经处∆M求得µi+1使其满足上述要求。其中Mi与Piiii−1i

&=M， M−M&=C理就排入长江流域的污水，M1010iii(i=11K20)，Ci就是当年

需要处理的污水量。具体结果如下：

年份 需处理污水（亿吨）年份需处理污水（亿吨） 2005 112 2010216.48

16

2006 130.48 2011241.38 2007 150.1 2012267.82 2008 170.92 2013295.88 2009 193.02 2014325.66

由上表容易看出，要达到水质的要求，必须处理以指数增长的污水，而 这也是不现实的，只有从根本控制污染源增长，才是解决问题的根本。

¾ 基于反馈型神经网络的水质预测与控制模型： 1.建模准备：

反馈型神经网络又称为递归网络或回归网络，是一种反馈动力学系统，比前向神经网络具有更强的计算能力。我们选用该网络对神经元进行训练，分析和预测题目中水质污染的发展趋势，并且在满足长江干流中各类水的比例要求下，给出相应的污水处理方案。

。该网络的特点是隐含层针对本题，我们选用 Elman 神经网络（结构图见下）

的输出通过承接层的延迟与存储，自联到隐含层的输入，这种自联方式使其对历史状态的数据具有很强的敏感性，内部反馈网络的加入增加了网络本身处理动态信息的能力，从而达到了建立动态模型的目的，同时该网络对非线性很强的系统也能得到很好的结果。

Elman 神经网络结构

2．模型求解：

(1)问题三的求解：

首先我们考虑以第一年的数据作为网络的训练样本，给出下一年总污水量和各个类别所占的比例的预测，即按照y(t)=f(y(t−1),u(t−1))依次类推得到未来十年的水质污染的发展趋势。实际计算中由上述方法得到的结果不理想，所以我们考虑采用时间窗移动法。

ny阶时间窗移动法：

每次选择前ny年的数据作为网络的训练样本，即将前ny个水文年的总污水量和

17

各个类别水所占的比例作为的输入向量，对第(ny+1)年的总污水量和各个类别水所占的比例作出预测。 对应的系统函数为：

y(t)=f(y(t−1),L,y(t−ny),u(t−1),L,u(t−nu))

注：

u(t)表示影响整个系统的控制因素。

y(t−1),L,y(t−ny)为系统前ny次输出对系统新的输出的反馈影响。

在具体解决问题三时，我们选择3阶时间窗移动法，中间层神经元的数目确定是8，得到了以下的结果： 长江全流域水质变化预测：

长江干流内水质变化的预测：

18

对于总量的预测：

年份 废水总量（亿吨）年份废水总量（亿吨） 2005 306.85 2010463.20 2006 328.62 2011512.45 2007 355.38 2012567.50 2008 385.29 2013625.07 2009 421.26 2014678.55 简单分析上述图表：我们可以得到，若不采取有效的治理措施。在未来的十年内，由于废水总量的不断增加，江河自然降解能力的，长江的水质在全流域范围将面临严重的污染。从我们对水质的简单分类上，也可以的到，引用水所占比例逐年下降，劣质水的比例逐年上升。

(2)问题四的求解：

同样采用3阶时间窗移动法，利用“问题三”中得到的20个水文年的污水总量，对问题四按以下步骤进行求解：

步骤一：对于前10个水文年，以每4年的数据作为网络训练的样本，即以4年水文年的

污水总量和前3个水文年中各个类别的水所占的比例作为网络的输入向量，以第4个水文年各个类别的水所占的比例作为目标向量，对神经网络进行训练。

步骤二：为了使未来10个水文年每一年各个类别的水所占的比例满足问题的要求，我们

构造了准则一：

准则一：按下图所示规则，判断每一个水文年内各个类别的水所占的比例，若满足要求，

则求解下一个水文年的各个类别的水所占比例；若不满足要求，则人为控制排放的污水总量，直到各个类别的水所占比例均满足要求。此处我们利用穷举的方法给出需要降低的污水量。

污水控制流程图

步骤三：记录步骤二中每年需要人为控制的污水量，则为未来十年内每个水文年所要处

理的污水量。（列表如下）

年份 需处理污水（亿吨）

年份

需处理污水（亿吨）

251 280 313 357 386 2005 125 20102006 147 20112007 173 20122008 203 20132009 229 201419

结果的简要分析：所需要处理的污水量在一年之内是恒定的，但随着年份的

推移，由于排污水量的不断增加，所需要处理的污水总量也是呈现出一种增加的趋势。

八．模型扩展：

¾ 模型扩展一：水质预测与控制的电路模拟模型

根据题目所给数据和规律，可将分析污染物的问题转化为分析电流电压的问题。 1．由河流拓扑结构提取电路图： (1)节点的设置：

每个观测点作为电路的节点。

(2)电流源的设置：

电流源代表污染源，由于各个干流观测点及入口支流皆可为污染源，故将各

个干流观测点及入口支流作为电流源。电流源的大小代表该观测点本地污染物注入量，电流源军非负值，为零表示该观测点不是污染源，各电流源总和等于长江流域的总污染排放量。 (3)电阻和导线的设置：

根据河流的拓扑结构在可以直接连通的节点之间用导线连接，代表各干流

河段和支流河段。在每个导线上串联一个电阻，大小基本与该河段上水体对污染物的降解比例成正比(我们在问题二的模型中已经得到这个比例的构造方法)。源头与入海口之间亦用导线连接，在该导线上串联一个无穷大的电阻，以表示二者之间没有直接联系。 2．合理性论证：

在用电路模型分析河流中污染物总量问题，最关键的问题是要保证每个河段上实际流过的污染物总量与电路模型中相应导线上的电流量对应相等。可如下论证这一点：

由于在一定流速条件下，污染物由一个观测点到另一个观测点过程中总量减少的比例是一定的，那麽通过调整各个河段的电阻大小，总可以使各导线上的污染物量与河段上实际流过的污染物量相同。这样对水质的分类即相当于对导线上电流值大小的分类。 3．模型应用：

对于问题二，相当于已知各导线上的电流值求各电流源的值；问题三则在求得各电流源值总和（即河流总流量）的变化趋势的条件下，讨论各节点的电流值；问题四则为调解各电流源的值，使各类电流值所占比例满足条件。由于时间的原因，我们并未进行实际求解。 4．模型评价：

优点：应用专业软件求解，电路分析模型可大大简化运算过程，同时可以真实

的反应污染物的分布。

缺点：各导线上电阻的取值对于电路模型的求解影响很大。

¾ 模型扩展二：房室模型

1．模型假设：

20

只有相邻类别的水可以相互转化。

2．符号约定：

在原有模型的基础上作如下的符号约定：

pi：第i个水文年的废水排放总量 pij：第i个水文年排放的废水转化为第j类水的量 Vi：第i个水文年河流总流量Vi：第i个水文年第j类水的流量Cj：第j类水的污染物浓度

j

θi：第i个水文年排放的废水中污染物的浓度

3．模型建立：

考虑各个水文年河流中A、B、C三类水之间的转化及废水排放量对于三类水所占比例的影响，如下图所示：

建立如下房室模型：

N

(1)pi=∑pij

j=1N

(2)Vi=∑Vij

j=1

jj

(3)Cj⋅(Vi+1−Vi)

= Cj⋅Vij⋅hjj+Cj+1⋅Vij+1⋅hjj +1+θi⋅Pji+1−Cj⋅Vij⋅hjj+1 (j=1,N)

jj(4)Cj⋅(Vi+1−Vi)

=Cj⋅Vij⋅hjj+Cj−1⋅Vij−1⋅hjj -1+Cj+1⋅Vij+1⋅hjj +1+θi⋅Pji+1−Cj⋅Vij⋅hjj+1 (j=2,...,N-1)

说明：

(1)式描述转化为各类水的废水之和等于该水文年的废水排放总量。 (2)式描述各类水的流量之和等于该年河流的总流量。 (3)式描述最优与最差水的变化规律。 (4)式描述中间类别水的变化规律。 4．模型应用：

21

对于本题，由于我们将水质分为三类，故对应为四房室模型，同时在此模型中我们讨论的皆为各类别水的流量所占的比例，而题目中给出的比例为河段比例，对于这个问题可按问题二中的方法将流量比例转化为河段比例。但由于数据量巨大及时间有限，我们未对该模型进行求解。 5．模型评价：

优点：模型适用于对水质进行任意多种分类的情况，并且考虑到了排污量与各类水

之间的内在联系，较为深入地对问题进行了描述。 缺点：模型对于参数的依赖性较大，求解复杂。

八．模型评价：

¾ 问题一中水质评价方法优缺点分析： 1．评价方法的优点：

(1)所建立的评价方法可方便地应用于河流的水质现状评价分析，同时也适用于对其他对象的综合评价。

(2)方法简单，评价符合客观实际,结果可信度高。

(3)不仅能够确定评价对象所属的级别，还能进行不同评价对象间的优劣比较。 2．评价方法的缺点：

该方法在确定各周期各个观测点在各个水期的水质情况、各个观测点各个周期的水质情况时皆是选取最差的一个水质类别，若某观测点在一个丰水期的5个月内水质皆为Ⅰ类，只有一个月水质为Ⅲ类，则由我们的方法，该观测点在该水期的类别为Ⅲ类，显然这样评价该观测点是不科学的，但是鉴于水质状况不会突变，此种情况出现的几率很小。

¾ 问题二中模型简化的评价：

在求解问题二的时候我们作了只考虑干流上观测点处作为污染源的简化假设。该假设可以使得问题的分析变得简单、明了，在求解污染源的过程中，也不需要再去考虑支流和流域上的其他污染源对干流水域的影响。在该假设的基础上，我们对全文的总模型进行了简化，能够较为快捷的给出主要的污染源位置。也可以大致的确定污染源所在的区域。

但是，在实际的生活当中，污染源是会分布在整个河流域上的，支流对干流上的影响也是不容忽视的。问题二的回答也就归属于对全文总模型的求解。这样可以更准确的给出各个主要的污染源的所在位置。

¾ 长江水质预测和治理模型的综合评价：

题目中的问题三和问题四，要我们对未来十年内，长江水质污染的发展做出预测分析。并在预测的基础上，给出相应的治理方案，以满足题目中对长江水体系中，每一类水的比例要求。

我们对问题三和问题四进行了综合考虑，建立了“动态马尔可夫链状模型”和

，两个模型各有特点： 基于Elman网络的“反馈型神经网络模型”

其中“动态马尔可夫链状模型”可以从本质意义上描述出不同类别的水比例之间的转换。但是，在计算的过程中，需要对马氏转移矩阵中的各个参数进行拟合，增加了误差的来源。

22

而“反馈型神经网络模型”可以更快捷的给出长江水质的发展趋势，可是在应用的过程中，需要设定中间层神经元的个数，而该值的确定只能从经验上去给出，无疑增大了解题的难度。另一方面，神经元在处理各类水比例预测的时候，对于神经元的训练，不能加入总比例之和为一的约束。因此，求解到的预测比例，求和会不满足一，这就需要在得到结果的时候，对同一时期的各个比例进行重新的归一化处理。

十．问题五的解决：

经过三天时间对长江水质的评价和预测问题的讨论，查阅资料时看到的长江污染情况的数字让人触目惊心，各界人士拯救长江矢志不移的决心更让我们深深感动。对长江水质的评价和预测仅仅是揭示问题的严重性，而真正解决长江水质污染问题才是根本目的。为此，我们提出以下建议：

1.保护生态环境，提高水体降解系数：

把水环境保护和治理作为环保工作的重点，逐步增加并确保水环境投入，通过植树造林提高长江流域的森林覆盖率，减少长江的泥沙含量，维持良好的生态环境，以提高水体的降解系数。

2.沿岸经济发展策略：

在确定经济项目时，要考虑水环境承受能力。严格控制工业及城市污水排放和农村化肥、农药的使用量，推广矿产资源开发废弃物最小量化和清洁生产技术，推进循环经济，减少矿产资源开发的环境代价。

可以从问题四的求解中，得到这样的结果，只是进行污水的处理不能从根本上解决问题。如果不对沿岸的经济发展进行相应得制约，不从根本上减少排污。这样，排入长江的污水是不能通过处理来实现。即使可以，这样的投入也要远大于对沿岸经济发展的规划。

同时，制定和实施有利于水环境保护的产业，调整产业结构和生产力布局，对重点水污染单位实行严格的管理措施，从根本上减少水环境污染。在一些厂矿相对集中的区域，建设污水处理设备，相对集中地处理污水，这样做所需管道线路短，处理费用较低。具体做法是：各厂矿共同投资兴建污水处理站，每年根据各厂矿的污水量、处理的难易程度收取处理费用，污水处理站负责处理污水至达标，同时供给部分中水至可利用环节（中水就是污水经过处理后达到国家规定水质标准的再生水，它可以作为公用设施和住宅冲厕、浇灌绿地、河湖景观、道路降尘、工业冷却、洗车用水等非人体接触的用水）。每年各厂矿的污水排放量、中水利用率、处理排污量等均由污水处理站负责上报水行政主管部门,由于与污水处理站的直接经济效益挂钩，数据会统计得校准确，水行政主管部门只需管住处理站一个环节，既减少了污染物排放量,降低了的治污投资,又便于管理。 3.加强立法：

加强水环境保制建设,制定和完善地方性法制法规,做到依法治理水环境污染。运用经济手段,促进水污染治理。依法征收排污费,逐步提高收费标准, 以强制单位尤其是工业污染源进行水污染治理。征收水资源保护补偿费以保护水资源,该费以对水资源影响程度的大小为收费依据，用于水体功能破坏治理和各种受害损失的补偿。另外，要立法对在长江治污过程中存在的现象进行严厉的打击。 4.提高忧患意识：

23

加强宣传和社会监督作用,通过一定的新闻媒体,设立专栏定期公开评定污染企业的环保等级、污染物排放量等,扬优贬劣,形成压力, 向民众宣传污染的危害性及治污的重要性，唤醒民众、企业、领导干部和全社会水资源保护意识的提高,不能因为长江水量较大而掩盖了问题的严重性，要使公众真正意识到问题的严峻性。 5.教育治污：

加大宣传和教育力度,加强青少年的环保意识。 6.加强水环境科学研究：

加强综合性、公益性强的重大水环境保护及示范性工程的建设,重点开展高浓度有机废水处理技术、污水资源化技术和湖泊富营养化控制技术等水资源问题和应用技术的研究，应用科学的手段有效的控制长江的污染情况。

奔腾不息的长江，孕育了古老的华夏文明，哺育着亿万中华儿女，为了让母亲河重现昔日的雄奇秀丽、气势磅礴，让我们共同努力掀起治理长江之污的风暴！

十一．参考文献：

[1]姜启源、谢金星、叶俊，数学模型，北京：高等教育出版社，2003。 [2]陆书玉，环境影响评价，北京：高等教育出版社，2004。

[3]刘昌明、岳天祥、周成虎，地理学的数学模型与应用，北京：科学出版社，2000。 [4]宋新山、邓伟，环境数学模型，北京：科学出版社，2004。

[5]李涛、贺勇军、刘志俭，Matlab工具箱应用指南—应用数学篇，北京：电子工业出版社，2000。

[6]飞思科技产品研发中心，神经网络理论与Matlab7实现，北京：电子工业出版社,2005。 [7]虞孝感，长江流域可持续发展研究，北京：科学出版社，2002。 [8]汪斌，水环境保护与管理文集，河南：黄河水利出版社,2002。

[9]邹宇庭、郑晓练、缪旭晖，SARS传播的数学原理及预测与控制，工程数学学报， 第20卷第7期：29页，2003。

[10]史戈宇、谢必克、毛燕杰，电路模拟方法在数学建模中的应用，工程数学学报， 第21卷第7期：43页，2004。

[11]话说长江，http://www.cjw.com.cn:6666/index.asp，2005年9月18日。

24