长江水质的评价和预测

摘要：

本文是关于长江水质评价和对未来趋势进行预测分析的问题，通过对长江近两年来相关数据的分析，建立了以下评价和预测模型。

对于问题一，首先对两年来长江的水质进行了综合评价：对数据进行归一化处理，利用熵值确定各类污染物对水质影响的权重，得到水质的综合评测指数Q，再利用两年来各月的Q值与设定的综合测评指数衡量标准进行对比的结果，定义出反映水质优劣程度的程度系数，再对程度系数进行区间划分，作为水质分级的指标。最后计算得出程度系数0.7109，按照水质分级的指标，得到结论：两年来长江的综合水质为?良?类，但属于“良”类中较差的水平。然后分析两年来各地区的水质状况：采用模糊集对模型，以集对分析为基础,重视信息中的相对性、模糊性,综合评价两年来各观测站污染情况及其变化。

对于问题二，分析可知某站点的污染物质量主要来自该站点的排放量与上游站点流经的变化量之和，通过一维河流的稳态水质模型，确定干流上污染物的变化量，计算出各地区两种污染物的排放质量，确定高锰酸盐指数（CODMn）的主要污染源在湖南岳阳、湖北宜昌、江苏南京和江西九江等地区，氨氮(NH3-N)的主要污染源在湖南岳阳江西九江和湖北宜昌等地区，与实际情况比较相符。

对于问题三，运用灰色模型预测法，以前十年的废水排放量的数据为依据进行预测，经比较验证了模型合格。随后，对未来十年水文年内全流域各类水质所占流域河长百分比进行了趋势预测，预测结果表明：若不对污染采取有效措施，长江水质将会不断恶化，十年之后，长江全流域内将有60%以上的水不可饮用。

对于问题四，水质所占河长百分比之间有密切的联系，我们通过分析，利用多元线性回归模型来构造两者之间的关系，利用spss 软件可以快速做出预测模型并能对异常点（如第III 类水质）进行剔除，使得模型更加合理。年排污量与总流量之比和各类水质的线性关系也就明确，利用问题三预测的数据，在满足条件下，可以预测出2005-2014 年的污水处理量。

对于问题五，本文以客观的态度，将长江的现状结合文中所建模型的结论分析，给出了治理长江污染问题的若干点建议。

关键词：熵值法 模糊集对模型 灰色预测 多元线性回归

1.问题重述

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

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

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

（1）对长江近两年多的水质情况做出定量的综合评价，并分析各地区水质的污染状况。

（2）研究、分析长江干流近一年多主要污染物高锰酸盐指数和氨氮的污染

源主要在哪些地区?

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

江未来水质污染的发展趋势做出预测分析，比如研究未来10年的情况。 （4）根据你的预测分析，如果未来10年内每年都要求长江干流的Ⅳ类和Ⅴ

类水的比例控制在20%以内，且没有劣Ⅴ类水,那么每年需要处理多少污水？

（5）你对解决长江水质污染问题有什么切实可行的建议和意见。

2.模型假设

（1）长江的降解指数是近似均匀的，取整个流域的降解系数恒为0.2（单位：1/

天）。

（2）长江流水水质是一维稳定河流水质。 （3）假设未来十年不采取任何排污措施。

（4）支流均是在位于长江干流上的6处观测点处汇入长江主干流的。 （5）一个观测站的水质污染主要来自于本地区的排污和上游的污水。 （6）长江干流的自净能力是近似均匀的。

3.符号说明

n m R 观测站的个数 评价指标的个数 建n个观测站的m个评价指标的判断矩阵 在同种评价指标下不同方案中的最满意者 在同种评价指标下不同方案中的最不满意者 评价指标的熵 评价指标的熵权矩阵 相对差 污染物在各个观测站的相对差之和 水质的综合评测指数 定义衡量标准 xmax xmin Hj W dij Dj Q Qs e Qq 水质某时期内长江总体污染的程度 第q个月的综合评测指数 测点集 指标集 决策矩阵 最优测点集 最劣测点集 F E D U V

akr ckr 同一度 对立隶属度 平均同一隶属度 平均对立隶属度 相对贴近度 n观测站污染源的排污量 n观测站观测到的i污染物的质量 来自n观测站上一个观测站的排污量 污染物的降解系数 两站点间的距离 河段的平均流速 污染物的质量变化 n站点i污染物的浓度 观测站的水流量 两站点之间的平均距离 流域段长度 江水完全流过水区域段 各段各污染物的残存系数 发展灰数 内生控制灰数 残差 相对差 回归常数 ak ck rk min Min Mi*(n?1) k x v m(x) cin Vn vn L Sj K a u E(k) e(k) ?0 ?Y

?长江年处理污水量与长江的年总流量

4 问题分析

本题针对长江水的污染问题，提出对长江水质现状的评测以及对长江水质未来状况的预测问题。我们认为对长江水质的综合评价应该从各类污染物对水质的影响出发，从客观的实测数据出发，通过数据的熵处理，给定各类物质各一个权值反映其对水质影响大小，利用相对比较的原理进行评测，对一个时期内的水质的综合分析采用把大时期划分成小时期的方式考虑；对各地水质污染状况的评测，可以采用平均监测值的方法去判断一个时期内此地区的水质状况；考虑污染源的定位问题可以考虑上游观测站与下游观测站的污染物浓度存在差值的原因是污染源的污水排放量与江水自净能力，在这些数据中只要知道上游观测站浓度、下游观测站浓度、江水自净系数即可求出在观测站之间江水段上注入的污水量，此污水量相对越大就说明此江水段上存在主要的污染源；对于预测，可以考虑到将来的值与过去、现在值的联系，因为江水自身的特性，在枯水期与丰水期污染情况相差较大，因此，我们对江水进行了枯水期、丰水期、水文年预测，采用灰色系统预测法找到数据间的联系，并做出误差分析；预测处理污水量的问题，水质所占河长百分比之间有密切的联系，我们通过分析，利用多元线性回归模型来构造两者之间的关系，利用spss 软件可以快速做出预测模型并能对异常点（如第III 类水质）进行剔除，使得模型更加合理。年排污量与总流量之比和各类水质的线性关系也就明确，利用问题三预测的数据，在满足条件下，从而可以预测出2005-2014 年的污水处理量。

5.模型的建立及求解

5.1问题一的求解 5.1.1综合评价情况

本题要求对长江近两年多的水质情况做出定量的综合评价，然后再分析各地区水质的污染状况。附表3记录了近两年来，长江流域上的17个观测站中4类物质(DO，CODMn,NH3-N，PH)浓度的监测值。通过分析2003年3月至2005年6月长江流域主要城市水质检测报告数据知，在28个月中17个观测站点的PH值几乎均在6~9之间，而国标(GB3838-2002)给出的《地表水环境质量标准》中Ⅰ类、Ⅱ类、Ⅲ类、Ⅳ类、Ⅴ类PH值的标准限值都是6~9，因而在评价各站点的水质类别时未考虑PH值这一指标。

a．首先用17个观测站对3项污染物的监测值构成的决策矩阵。

这里可以类比信息论中的熵概念。在信息论中，信息熵反映了信息无序化的程度，熵越小，系统无序度越小，信息效用值越大反之亦然。对于所讨论的n

个方案的m个评价指标的初始矩阵，决策矩阵显然是一种信息的载体，故可用信息熵所获系统信息的有效度及效用，由于信息熵可用来度量m个指标的信息效用值，因此，利用熵值法计算每个指标的权重，其本质就是利用该指标信息的效用值来计算的，效用值越高，其对评价的重要性越大。步骤如下：

1）构建n个观测站的m个评价指标的判断矩阵

R=(xij)mn（i=1,2…n, j=1,2…m）

2）将判断矩阵R归一化处理，得到c矩阵

ijcij?xij?xminxmax?xmin

式中，xmax、xmin分别表示在同种评价指标下不同方案中的最满意者或最不满意者（根据不同的情况）。例如，在水质评价中，DO的含量越大越好，而其他的评价指标则基本都是越小越好。

根据熵的定义, n个方案的m个评价指标,可以确定评价指标的熵为:

1n Hj?? j=1,2,…m） (?fijlnfij(i=1,2,3…n,)lnni?1在上式中，有

fij?1?cijn

ij?(1?c)i?1计算评价指标的熵权矩阵公式W：

W?(?j)1?m

其中矩阵元素的计算公式为?j?1?Hjm??Hjj?1m ，且满足??j?1

j?1mb．下面定义的是标准值I3(水质的标准值):

因为Ⅰ、Ⅱ、Ⅲ类水为饮用水，因此以Ⅲ类水的标准值作为衡量限。定义

相对差dij：

??dij?xij?I3 ?d?I?x?3ij?ij

上式中，第一个式子用来描述DO这类随评价指标值的递增的污染物； 第二个式子用来描述CODMn,NH3-N这类随评价指标值的递减的污染物。

显然dij值越大意味i评测站的j类污染物对水质的污染越小，满意度越高。 对相对差的求和：

n Dj??di j

i?1Dj指j污染物在各个观测站的相对差之和，Dj可以整体反映j污染物对长

江水质的影响。定义水质的综合评测指数Q:

Q??Dj?j

j?1mQ是从各类污染物所占权重与求和相对差的值的角度全面考虑长江每月水

质，通过进一步计算，可以求出近两年来长江的综合水质情况。

定义衡量标准Qs：

(Ⅰ?j?Ⅱj)?QS??Ⅱj???n??j 2??

其中Ⅰj,Ⅱj是指第j类污染物的一类水和二类水的标准值,n指观测站的个数,?j指j污染物两年来各月权重的平均值。

衡量标准Qs合理性分析:通过在Ⅰ,Ⅱ类水的标准值之间找到一个点的值与

Ⅲ类水的标准值的差与相对差做比较可以大体评判水质的好坏,我们希望水质情况的目标是Ⅰ,Ⅱ,Ⅲ类水,因此在Ⅰ,Ⅱ类水的标准值之间找到一个点的值与Ⅲ类水的标准值的差来考虑具有合理性。

定义程度系数：

k?(Q e?q?1kq?Qs) (6)

q?Qq?1?Qs

上式中，e表示水质某时期内长江总体污染的程度；

k表示要调查的月份数(本题k=28)；

Qq表示第q个月的综合评测指数；

分子的意义：各个月的实际综合评测指数Qq与衡量标准Qs的差求和,因为每月Qq与衡量标准Qs的差有正有负,求和反映了k个月水质达标(大于0,为正数)与不达标(小于0,为负数)在量上的相对差。

分母的意义：是综合评测指数Qq与衡量标准QS的差的绝对值求和,这个值为正,反映了水质达标(大于0,为正)与不达标(小于0,取绝对值,为正)在量上的总和。

比值的意义：水质达标与水质不达标的相对优势。

从表达式中可以看出:

e?[?1,1]，按照我们的模型用e的均匀分配把水质综合评价定义为以下几

类比较合理:

e?[?1,0] 定义为崩溃,根据上述模型可知Ⅰ，Ⅱ，Ⅲ类水基本上不存

在

e?（0，1/6] 定义为?差? ； e?（1/6,1/3] 定义为?较差?； e?（1/3,1/2] 定义为?中?； e?（1/2,2/3] 定义为?较良?； e?（2/3,5/6] 定义为?良?； e?（5/6,1] 定义为?优?。

c.模型应用与评价：

可以看出相邻水质级别在量上的差值不大，比如从?优?变为?良?，e上变化量并不太大，这与实际情况比较相符。

长江两年水质综合评价图:

长江两年的水质综合评价504540综合评价指标3530252015051015年月份202530

横坐标1代表2003.6,2代表2003.7,3代表2003.8依次类推； 纵坐标代表Q值；直线是Qs值

结果与分析：

?1=0.2707

?2=0.3327 ?3=0.3965

Qs=0.2707*17*1.75+0.3327*17*3+0.3965*17*0.675=29.5733

e= 0.7105

所以在这两年来看长江的水质是‘良’类水，但是e的值更加接近2/3，即在程度上是‘良’类水中较差的水平，意味着很快进入‘较良’类水的范围。 1.1.2．各地区水质污染的状况及其变化评价

对于长江各地区水质污染状况的分析与评价，现采用模糊集对模型，综合评价各地区的水质情况及其两年多的变化趋势。

评价原理及其步骤：

设多属性决策问题Q?{F,E,D}，其中：

，fk 为第k个测点； F?{fk} (k=1,2,…,m) 为测点集（月份）

E?{er} ( r=1,2,…,n ) 为指标集，E有不同的类型指标，记E1为收益

型指标（即越大越好，比如溶解氧），E2为成本型指标（即越小越好，比如高锰酸盐指数等）；

D?{dkr} 为决策矩阵，dkr为测点fk关于指标er的属性值。

记最优测点集U={u1,u2,…,ur}，最劣测点集V={v1,v2,…,vr}，其中：ur，

vr分别为指标er的最优值和最劣值。

对于er?E1 ，比较区间为{vr,ur}，定义集对{dkr,ur}同一度akr和对立隶属度ckr如下：

akr?dkr ur?vrurvr

(ur?vr)dkrckr?

同理，对于er?E2，在比较区间[ur,vr]可得：

akr?urvr

(ur?vr)dkrckr?dkr ur?vr记平均同一隶属度为ak，表示对fk接近最优测点U的肯定程度；平均对立隶属度为ck，表示对fk接近最优测点U的否定程度，可得：

1nak??akr

nr?1

1nck??ckr

nr?1n为指标集E的元素个数

由于ak、ck是相对确定的，那么在相对确定条件下可定义fk与U的相对贴近度为：

rk?ak ak?ck由上式可得：rk越大则ak越大或Ck越小，对fk接近最优测点U的肯定程度就越高，或对fk接近最优测点U的否定程度就越低，故rk越大表明水质综合质量最好，可根据rk的大小，排序各月份综合水质。

根据以上模型，以湖北宜昌南津关为例，选取溶解氧、高锰酸盐指数、氨氮三项指标计算该点2003年6月至2005年9月这28个月中rk值的变化情况，并依此分析其水质变化。

表1湖北宜昌南津关各指标测定值统计结果

指标 月份 2003-06 2003-07 2003-08 2003-09 2003-10 2003-11 2003-12 2004-01 2004-02 2004-03 2004-04 2004-05 2004-06 2004-07 2004-08 溶解氧(DO) 7.81 7.89 6.65 10.6 11.9 9 9.63 8.17 8.38 11.6 7.6 9.76 7.96 7.82 6.94 高锰酸盐指数(CODMn) 5.8 5.2 2.8 3.6 2.4 3.2 1.8 1.9 3.5 3.6 2.2 3 3.8 3.2 3.1 氨氮(NH3-N) 0.55 0.22 0.31 0.36 0.26 0.22 0.2 0.33 0.25 0.43 0.27 0.3 0.34 0.16 0.19

2004-09 2004-10 2004-11 2004-12 2005-01 2005-02 2005-03 2005-04 2005-05 2005-06 2005-07 2005-08 2005-09 7.56 8.4 9 9.07 9.2 10.4 9.6 8.22 7.88 7.1 7 6.5 6.51 3.4 3.7 1.9 2.4 1.9 2 2.1 2.4 2.4 2.1 2.6 1.3 3.2 0.29 0.23 0.24 0.16 0.22 0.18 0.19 0.13 0.22 0.17 0.15 0.63 0.2

根据上述模型，通过Matlab程序，确定U?(11.90, 1.30, 0.13)，

V?(6.50, 5.80, 0.63)，得Matlab程序Rk，代入数据计算得下表：

表 2 湖北宜昌南津关水质变化趋势的平均同一隶属度、相对贴近度及其排序等

平均同一隶属度平均对立隶属度相对贴近度rk 0.55 0.22 0.31 0.36 0.26 0.22 0.2 0.33 0.25 0.43 0.27 0.3 0.34 根据rk排序 1 13 2 3 16 27 10 4 9 28 17 15 12 月份 2003-06 2003-07 2003-08 2003-09 2003-10 2003-11 2003-12 2004-01 2004-02 2004-03 2004-04 2004-05 2004-06

ak 0.803488 1.122862 1.088312 1.170421 1.503701 1.310829 1.65217 1.329509 1.189908 1.176039 1.294881 1.243636 1.029025 ck 5.8 5.2 2.8 3.6 2.4 3.2 1.8 1.9 3.5 3.6 2.2 3 3.8

2004-07 2004-08 2004-09 2004-10 2004-11 2004-12 2005-01 2005-02 2005-03 2005-04 2005-05 2005-06 2005-07 2005-08 2005-09 1.430386 1.28692 1.094811 1.212077 1.497076 1.608943 1.548765 1.694888 1.594615 1.718175 1.360582 1.525471 1.507307 1.341215 1.224486 3.2 3.1 3.4 3.7 1.9 2.4 1.9 2 2.1 2.4 2.4 2.1 2.6 1.3 3.2 0.16 0.19 0.29 0.23 0.24 0.16 0.22 0.18 0.19 0.13 0.22 0.17 0.15 0.63 0.2 11 6 8 24 14 26 25 18 5 20 19 22 23 7 21 通过表2数据得下图：

0.650.60.550.5相对贴近度0.450.40.350.30.250.20.15051015第i个月202530

图1 湖北宜昌南津关相对贴近度随时间变化图

图1及其表2反映了湖北宜昌南津关28个月内的水质变化情况，其结果表

明，2003年6月份相对贴近度最大，即rk值为0.6827，排在第一位，表明2003年6月份水源水质综合质量最好，其次为：2004年10月，2005年1月，2004年6月，2003年11月。2004年8月的水质综合质量最差，其次为：2003年7月，2004年1月，2004年9月， 2003年10月，可看出湖北宜昌夏季水质综合质量较差，而春季较好。

同理可求出干流其余6个点位及其支流上的10个点位的相对贴近度随月份的变化情况，见附表。 5.2问题二的求解

本问题主要研究和分析长江干流近一年多主要污染物高锰酸盐指数和氨氮，最终确定其污染源在哪些地区，为解决此问题，建立了以下模型。 5.2.1 污染物的污染源评价模型的建立

通过对原始数据的分析可知，观测站的水质污染的来源，主要来自本地区的排污部分以及上游流经过的污水。那么，在观测站n的污染源i污染物每秒排放量即目标函数为：

* min?Mi?nM( ) ?i1n (i=1,2;n=1,2,3…7)

其中min为n观测站污染源的排污量；Min为n观测站观测到的i污染物的质量；Mi*(n?1)为来自n观测站上一个观测站的排污量。

对于长江干流上的七个主要观测站：攀枝花、重庆、 宜昌、岳阳、九江、安庆、南京，可分别称其为A1、A2、A3、A4、A5、A6、A7。

对于攀枝花至重庆、重庆至宜昌、宜昌至岳阳、岳阳至九江、九江至安庆、安庆至南京的长江干流流域段(River way Section)，可分别定义为：S1、S2、S3、

S4、S5、S6。

为了确定Mi*(n?1)，即来自上游的污水在流动过程中污染物的质量变化情况，我们对各情况进行以下分析来确定各约束条件：

河流对污染物有一定的净化能力，那么水流从上游流到下游的过程中污染物的含量会发生变化，当流量不变时，由一维河流的稳态水质模型，在下游n处污染物的浓度为：

ci(x)?c?e?k?t(x)?c?e?k?x/v

其中c为上游起点的污染物浓度；k为污染物的降解系数；x为两站点间的

距离；v为该河段的平均流速。

由质量与浓度、体积的关系式：m?cV，可得下游x处的污染物质量：

m(x)?c(x)?V?c?V?e?k?x/v

考虑到各段流量变化的情况，但经分析发现河流流量的变化对污染物的质量变化产生的影响极小，可以忽略不计。

从而确定污染物的质量变化为m(x)，它的大小只与位置有关。 所以上游的观测站污染物到达下游相邻的观测站后的质量变为：

Mi*(n?1)?Mi(n?1)?e?k?x/v

其中x为两观测站的距离，v为该河段的流速。

同时，观测站污染物的质量也可以通过该观测站的水流量与污染物的浓度来确定：

Min?ci?nV n其中cin为n站点i污染物的浓度，Vn为观测站的水流量。

简化上面公示，可得本地区的污染物源产量为：

min?Min?Mi*(n?1)?Min?Mi(n?1)?e?k?x/vn?cin?Vn?ci(n?1)?Vn?1?e?k?x/vn 其中 vn为两站点之间的平均距离； 式中i=1,2；n=1,2,3…,7 。 5.2.2 模型的求解

由题中假设，长江干流的江水做匀加速运动(a≤0或a≥0)。我们定义S1、S2、

S3、S4、S5、S6这6段干流流域段的平均流速矩阵Vij'(13?6)如下：

Vij'?Vij?V(i?j1)，13; j?1,，…6 ? ? i= 1,…2注：1）其中Vij为流域段Sj上游的观测站Aj在i时段的水流速。

2）其中Vi(j?1)为流域段Sj下游的观测站Aj?1在i时段的水流速。

3）其中i为从2004年4月至2005年4月的13个时段。

由题目中所给的原始数据经分析处理，结合公式，运用MATLAB编程（程序见附录）得到速度矩阵Vij'如下表：

表3 速度矩阵表 单位：m/s

流域段 时段 2004．04 2004．05 2004．06 2004．07 2004．08 2004．09 2004．10 2004．11 2004．12 2005．01 2005．02 2005．03 2005．04 S1 S2 S3 S4 S5 S6 2.9 2.8 3 3.2 2.95 4.95 2.7 2.3 2.3 1.8 1.5 1.4 1.65 1.5 1.35 1.65 1.85 1.75 3.25 1.9 1.3 1.1 1 0.7 0.65 0.8 0.9 0.85 1.25 1.45 1.4 1.8 1.55 0.75 0.75 0.55 0.5 0.5 0.45

0.95 1 1.4 1.5 1.45 2 1.6 0.85 0.8 0.65 0.65 0.7 0.6 1.05 1.1 1.5 1.55 1.6 2.75 1.65 0.9 0.8 0.7 0.7 0.8 0.75 1.15 1.15 1.55 1.65 1.7 3.4 1.8 0.95 0.85 0.75 0.75 0.85 0.8 定义S1、S2、S3、S4、S5、S6的流域段长度(Length)分别为：L1、L2、L3、L4、

L5、L6，

Li??Ln??Ln

n?1n?1ii?1由问题原始数据累积差可得：

L??950,778,395,500,164,464?

通过L与矩阵Vij'的计算( 过程见附录MATLAB程序)可以得到时间矩阵

Tij(13?6)：

?t11?t1n???Tij??????(m?13,n?6) ?t??m1?tmn?其中tij为i时刻，江水完全流过水区域段Sj所需的时间(单位：天)。

表4 时间矩阵表

流域段 时段 2004．04 2004．05 2004．06 2004．07 2004．08 2004．09 2004．10 2004．11 2004．12 2005．01 2005．02 2005．03 2005．04

S1 S2 S3 S4 S5 S6 3.7915 3.9269 3.6651 3.4361 3.7272 2.2213 4.0724 4.7806 4.7806 6.1085 7.3302 7.8538 6.6639 6.0031 6.6701 5.4574 4.8674 5.1455 2.7707 4.7393 6.9266 8.1860 9.0046 12.8638 13.8533 11.2558 5.0797 5.3785 3.6574 3.1529 3.2655 2.5399 2.9495 6.0957 6.0957 8.3123 9.1435 9.1435 10.1595 6.0916 5.7870 4.1336 3.8580 3.9910 2.8935 3.6169 6.8083 7.2338 8.9031 8.9031 8.2672 9.6451 1.8078 1.7256 1.2654 1.2246 1.1863 0.6902 1.1504 2.1091 2.3727 2.7116 2.7116 2.3727 2.5309 4.6699 4.6699 3.4648 3.2548 3.1590 1.5795 2.9835 5.6530 6.3181 7.1605 7.1605 6.3181 6.7130 由于河流对污染物有一定的净化能力，下面求解污染物经过一段水流后的残存浓度，首先获得各段各污染物的残存系数K＝e?k?x/vn：

由相关资料，我们确定长江对CODMn的降解系数k＝0.25，对NH3－N的降解系数k＝0.40，（单位1/天）。

对以上数据运用MATLAB编程可获得各段各污染物的残存系数K（i,j）,程序及结果均见附录。

根据残存浓度，很容易算出各时期各流段的污染物质量变化量： Mi*(n?1)?Mi(n?1)?e?k?x/v

运用MATLAB程序（见附录），分别得到CODMn和NH3-N两种污染物的质量变化量，绘制成表格如下：

表5 CODMn污染物质量变化量 单位：kg/s

时段 2004.04 2004.05 2004.06

四川攀枝花 3.2892 5.9931 4.0100 重庆朱沱 湖北宜昌 湖南岳阳 江西九江 安徽安庆 10.7689 5.4387 10.8864 12.9753 15.4817 30.9159 18.4232 51.8596 17.3675 60.0101 28.1434 55.8988 25.7012 33.8534 21.5995

2004.07 2004.08 2004.09 2004.10 2004.11 2004.12 2005.01 2005.02 2005.03 2005.04

4.7373 8.5432 21.9845 0.9422 1.2712 0.4605 0.1855 0.0881 0.0962 0.1335 16.0284 5.8569 7.9264 2.8920 1.1000 0.5502 0.2602 0.2821 0.4378 33.0255 32.8868 33.8061 4.3876 3.8691 1.0869 0.9172 1.1061 1.0222 38.5822 51.6864 38.1976 54.8817 31.5990 35.3430 5.6879 6.1384 4.5103 2.7574 3.6539 1.8834 18.9578 16.7761 13.2808 16.2101 14.2232 16.0406 24.0357 28.8079 30.8779 9.1013 4.9396 5.3387 6.1751 11.9069 8.6742 104.7706 96.3985 101.7858 116.3975 79.9889 表6 NH3-N污染物质量变化量 单位：kg/s

时段 2004.04 2004.05 2004.06 2004.07 2004.08 2004.09 2004.10 2004.11 2004.12 2005.01 2005.02 2005.03 2005.04

四川攀枝花 0.1215 0.0541 0.0370 0.0472 0.8421 0.2324 0.0512 0.0133 0.0112 0.0043 0.0049 0.0070 0.0045 重庆朱沱 湖北宜昌 湖南岳阳 江西九江 安徽安庆 0.2626 0.2454 0.2881 0.4213 0.2301 2.5143 0.5597 0.1381 0.1264 0.0592 0.0115 0.0100 0.0255 0.7433 0.6909 1.5982 1.0290 1.2350 5.6174 1.3501 0.2221 0.1034 0.0362 0.0209 0.0254 0.0121 0.5821 0.6480 1.5139 1.8540 1.7318 5.0728 1.8892 0.2836 0.1837 0.0977 0.0975 0.1083 0.0642 3.6815 3.2876 2.6674 2.6470 4.7932 13.8091 2.9741 1.3816 0.6132 0.7737 0.4526 0.6089 0.7682 1.0023 1.3127 1.3648 1.9090 1.5104 5.5226 1.5978 0.2126 0.1577 0.1403 0.1406 0.3932 0.2067 欲求此问题必须求出每个站点的排污量，即要用到公式：

min?Min?Mi*(n?1)

然后对每个站点13个月的排污量求平均数，即：

113min??mikn （i=1,2;n=1,2,3…,7）

13k?1

其中mikn表示i污染物在n站点第k个月的排污量。

运用MATLAB编程，获得结果绘制成表格如下：

表7 各站点污染物排放量 单位：kg/s

四川攀 CODMn 排放量 NH3-N 排放量

0.4816 8.9861 枝花 重庆朱沱 湖北宜昌 湖南岳阳 江西九江 安徽安庆 江苏南京 33.1952 38.0456 49.8466 35.2169 19.7911 35.8746 2.9356 3.7655 5.4107 3.9606 2.0693 1.2893 为了比较直观的分析数据，分别对各站点的CODMn 排放量和NH3-N 排放量利用Excel绘制折线图如下：

CODMn 各站点排放量60排放量(Kg/s)504030201001234站点

图5-2-1 CODMn 排放量图

567

NH3-N 各站点排放量65排放量(Kg/s)432101234站点567

图5-2-2 NH3-N 排放量图

5.2.3 结果分析

观察图5-2-1 CODMn排放量图，可知高锰酸盐指数的污染源主要集中在湖南岳阳、湖北宜昌、江苏南京和江西九江等地区，其中湖南岳阳排放量最大，居各排放点污染源之首。

观察图5-2-2 NH3-N排放量图，可知氨氮的污染源主要集中在湖南岳阳、江西九江和湖北宜昌等地区，其中湖南岳阳排放量最大。

从以上结果不难发现两种污染物的污染源都主要集中在湖南岳阳、江西九江和湖北宜昌。而实际情况是：宜昌至岳阳河段属于长江中游地区，该地区沿江工业以耗能大、污染大的重化工业为主。该地区目前已形成了以钢铁为主的治金工业、以汽车为主的机械工业、以石油化工和磷、盐为主的化学工业、以水泥为主的建材工业、以纺织和食品为主的轻纺工业等五大支柱产业，岳阳地区还是我国中南地区最大的石油化工生产基地和全国最大的客车空调机、抗菌素生产基地，这些沿江的重工业给长江带来了巨大的污染。而处于长江中下游地区的岳阳至九江河段也属于重污染河段，这些地方所排放的工业废水是导致水质污染的主要原因。

由此我们判定建模及计算结果的准确可靠性。 5.3问题三的求解

本题主要研究在不采取任何措施的情况下，根据1995~2004年长江流域水质报告统计数据，对未来10年水质污染的发展趋势做出预测分析。由于所给样本数据较少，故采用灰色预测模型进行预测。又水文年是指与水文情况相适应的一

种专用年度，其数据是一年各时期数据的平均数，可以更好地反映一年的整体水质情况，因此，在本题中我们采用的原始数据只有水文年全流域各水质类别的河长百分比，然后以这些数据为初始数列，对未来10年内不同水质的河长的比例进行预测。

首先建立灰色预测模型：

1）设处理后的1995年至2004年在不同时期的不同水质的河长比例的数据列为：

x(0)?{x(0)(1),x(0)(2),?,x(0)(10)}

?{?x(0)(j)i?1,2?,10}； 经一次累加，其数学表达式为：x（i）(1)j?1)),x得到： x(1)?{x((11(1)i(?2),x,(1)(2 0)}；

dx(1)?ax(1)?u； 设x满足以下白化形式的一阶常微分方程： dt(1)其中：a为发展灰数；u称为内生控制灰数； 白化形式的微分方程的解为

uux(1)(t)?[x(0)(1)?]e?at?

aauu2）将时间响应函数离散化，故有：x(1)(k?1)?[x(1)(1)?]e?ak?；

aa3）通过最小二乘法来估计常数a与u，因x(1)(1)已知，?t?(t?1)?t?1,则有：

?x(1)(2)?x(1)(n)(1)(1)(1)(0)??x(2)?x(2)?x(1)?x(2),类推得：?x(0)(n)； ?t?t1因x(i)(i)=[x(i)(i)?x(i)(i?1)],(i?2,3,...,20)，则有

2(1)(1)?x(0)(2)???11?2[x(2)?x(1)]??(0)??1(1)(1)?[x(3)?x(2)]1x(3)??a?.； ????2???????1??u???(0)??1(1)(1)?[x(20)?x(19)]1x(20)?????2?4）令y?(x(0)(2),x(0)(3),?,x(0)(20))?，?为转置符号

(1)(1)??11?2[x(2)?x(1)]?1(1)?(1)?[x(3)?x(2)]1?a?2?? 令B=,U???, 故可得：y?BU； ????u???(1)(1)1?[x(20)?x(19)]1???2?

????a??(B?B)?1B?y 而y?BU的最小二乘估计为：U?u????预测值由累减生成处理求得：

?(k)?x(k)?x(k?1) x?(1?ea)(x(0)(1)?u/a)e(?a)(k?1)，其中k?2,3,??????,n

(0)?(1)?(1)由上式，采用matlab编程，可以预测出未来十年的废水排放总量数值，程序见附录，结果如下表所示。

表未来十年废水总量（亿吨）

年份 废水总量 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 300.3 319.5 340.0 361.8 385.0 409.7 436.0 463.9 493.7 525.3 为了进一步看出未来十年废水排放总量的趋势，我们用EXCEL做出废水排放总量折线图，如图所示。

未来10年废水排放总量预测图600500排放量(亿吨）40030020010002005200620072008200920102011年份(2005-2014)201220132014

图 未来十年废水排放总量预测图

从上面的可以清楚地看到，如果不采取任何措施，未来十年废水排放总量将成直线上升趋势。这必将会导致水污染的情况越来越严重，可饮用水量将会越来越少。

同样道理，按照这个模型，我们来预测未来十年水文年内全流域各类水质所占总流域的百分比，预测结果如下表

表 未来十年水文年内全流域各类水质所占百分比 类别 年份 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 Ⅰ 1.925866 1.476568 1.122411 0.845758 0.631637 0.467478 0.342844 0.249157 0.179442 0.128092 Ⅱ 38.77851 39.13572 39.15857 38.83972 38.18143 37.19641 35.90803 34.34971 32.56336 30.59721 Ⅲ 27.43371 24.84769 22.31305 19.86221 17.52358 15.32114 13.27397 11.39598 9.695657 8.176157 Ⅳ 11.21941 10.8821 10.4647 9.975552 9.424833 8.824358 8.18718 7.527078 6.857933 6.193078 Ⅴ 6.407494 6.829613 7.217309 7.560497 7.849683 8.076564 8.234608 8.319553 8.329749 8.266282 劣Ⅴ 14.23502 16.82831 19.72396 22.91626 26.38883 30.11405 34.05336 38.15852 42.37386 46.63918 由上表数据绘制折线图如下：

Ⅰ类水质所占百分比2.5Ⅱ类水质所占百分比4540235301.525120150.51050200520062007200820092010201120122013201402005200620072008200920102011201220132014

Ⅲ类水质所占百分比302520151050200520062007200820092010201120122013201412Ⅳ类水质所占百分比10864202005200620072008200920102011201220132014

Ⅴ类水质所占百分比劣Ⅴ类水质所占百分比95087654321020052006200720082009201020112012201320144540353025201510502005200620072008200920102011201220132014

结果分析：观测水文年的十年预测数据，发现全流域中，前三类可饮用水质的河长百分比都随时间呈下降趋势，而Ⅴ类和劣Ⅴ类水所占百分比在逐年上升。十年后整个长江全流域大约60%的河长被污染，不可饮用，因此，加强排污、防污工作力度刻不容缓。 5.3.3模型三的结果检验 1）预测值的检验 ①残差检验

残差：E(k)?x(0)?(0)(k)?x(k),?(0)k?1,2,3,?,n;(k)]/x(0)(k),k?1,2,3,?,n;

相对残差：e(k)?[x(0)(k)?x如果e(k)?0.2，则认为达到一般要求；如果e(k)?0.1，则认为达到较高要求。

②级比偏差值的检验

首先由参考数据x(0)(k)，x(0)(k?1)计算出级比?(k)，结合发展系数求出响应的级比偏差：

?1?0.5a???(k)

1?0.5a???(k)?1??如果?(k)?0.2，则认为达到一般要求；如果?(k)?0.1，则认为达到较高要求。

③检验预测结果

这里对1995年至2004年的废水排放总量的预测结果进行检验,结果如下表：

表：1995年至2004年废水排放总量的预测结果

年份 项目 实际值 预测值 1995 174 1996 179 1997 183 182.6 1998 189 1999 207 2000 234 2001 205 2002 256 2003 270 2004 285 282.2 174.0 171.6 194.4 206.8 220.1 234.2 249.2 265.2

残差 相对残差 极比偏差 0.0 7.4 0.4 -5.4 0.2 13.9 -29.2 6.8 4.8 2.8 0.000 0.041 0.002 0.028 0.001 0.060 0.142 0.027 0.018 0.010 0.000 -0.034 -0.041 -0.030 0.028 0.059 -0.215 0.148 -0.009 -0.008 由上表数据知，相对残差、级比偏差值均小于0.1，说明预测的结果达到了较高的要求，也说明了模型GM（1,1）对未来十年废水排放总量的预测结果是合理的。进而也说明模型GM（1,1）对未来十年水文年内全流域各类水质所占总流域的百分比的预测结果是合理的。

5.4问题四的求解

根据附录所提供的数据进行分析，可以知道，长江年废水排放量和年总流量会影响各类别的废水占长江河长的百分比，当废水排放量越多时，超标的长江河长就会越长；而当长江的年总流量增大时，超标的长江河长会变短。由此可以判断，各类别水质河长的长短能够反映出长江污水年排放量和长江的年总流量的大小。因此，我们考虑将各个类别的水质河长占总长的百分比作为对长江污水年排放量与长江年总流量之比的影响因子。通过分析，我们建立了多元线性回归模型。 a．理论预测模型

预测对象：长江污水年排放量与长江的年总流量的比作为因变量Y

诸影响因素：各级别水质的所占河长百分比为自变量Xj(j＝1，2，…，k)，(X4 代表第IV类和第V 类的污水对应的河长)

Y与各个X之间存在的线性关系，从理论上可表述为

Y??0??1X1??2X2??3X3??????kXk??

式中：?0是回归常数；?j(j＝1，2，…，k)虽未知但均为某一固定数值，称为回归系数；?是除Xj(j＝1，2，…，k)以外的，可以忽略的随机因素，称为随机干扰项，也称为残差项。 b．预测模型

回归余项?的变化往往无法预测，而预测实际上是通过主体部分来完成的，即

?Y?b0?b1X1?b2X2?b3X3?????bkXk

我们利用spss 软件的线性回归分析功能对1995-2004 年的每年长江干流水文年各类水质的所占河长的百分比做为历史数据进行回归分析，我们得到了各个因子的常量、系数，分别为：

b0?2.67389,b1?0.02312,b2?0.01181,b4?0.00384,b5?0.02323

代入预测模型，得Y?2.67389?0.02312X1?0.01181X2?0.00384X4?0.02323X5 模型的相关说明：

1．由于我们在处理第IV和第V类水的时候要求它们的量之和不能超过20%，所以我们在处理这两个数据的时候把它们并为一类进行处理。

2. 在对进行数据的分析中，由于第Ⅲ类污水的残差统计量是一个奇异值，系统自动剔除了它，所以我们的模型中剩下四个变量。 模型的求解：

题目中要求我们估计需要处理的污水量，且第IV第V类的河长百分比之和小于20%，劣V类为零；那么我们对式中的因子进行控制，得到减少的各类水质河长的百分比差值，代入式

??Y?2.67389?0.02312X1?0.01181X2?0.00384X4?0.02323X5

中就可以得到长江年处理污水量与长江的年总流量之比?Y。

??Y?0.00384X4?0.02323X5

?0,X4?0.2?X4??上式，?X4?0.2,X4?0.2

?X5?X5?考虑长江年总流量的平稳趋势，可以用前十年的长江年总流量的算术平均值来作为后十年的平均长江年总流量。 根据?Y的定义，推出以下关系表达式： 未来的年处理污水量=期望长江年总流量×?Y

我们以第三题用灰色模型预测出的未来十年的水文年内总流域所占河长的百分比作为数据，用于我们的多元线性回归模型中，计算出未来十年每年需要处理多少污水才能使长江的水质满足要求。结果见下表：

年份 污水处理量（亿吨） 2005 2006 2007 2008 2009

??3933.9 4533.1 5197.6 5925.7 6713.9

2010 2011 2012 2013 2014

5.5问题五

7555.9 8443.1 9364.8 10309 11261 长江流域是我国经济发展最快的地区之一，长江流域在我国经济建设中的地位更加重要。保护好长江水资源，不仅对长江流域的经济、社会发展有重要意义，而且随着南水北调工程的实施，也将对华北、京津地区乃至全国的经济与社会发展产生重要作用。

根据问题一所建立的模型分析表明1994－2004年长江水质基本上属于良好水质，但是问题三中灰色预测数据已经表明长江水质正逐渐恶化并将继续下滑。这一情况应引起我们足够的警惕，防微杜渐，我们应尽早采取相应措施，因此我们提出了以下的建议和意见：

（1）长江的污水处理是刻不容缓的，根据预测情况，每拖一年，长江的污染状

况将有很大的恶化，超标的河长将有大幅度的提升，所以必须尽快地采取保护措施。

（2）虽然长江近两年的水质情况比较合格，但是长江污染的发展情况还是令人

担忧的，不要看到现状，就忽略了将来。所以，我们必须增强长江沿岸居民的环保意识。

（3）从问题中我们可以知道长江水有自净能力，影响长江的降解系数主要是长

江水系生态，我们不能人为破坏长江的水系生态，具体措施是合理的培植水生植物（对自净有利的植物），在长江沿岸植树，严禁乱砍乱伐树林，提高长江沿线森林覆盖率，保护长江流域生物的多样性，进而提高水环境自净能力。

（4）因为长江的污染深受废水排放量的影响，所以应尽快在长江流域建立水源

保护区，保护区沿岸严禁兴建排污量大的工厂，更不能将污水直接排入河道。

（5）政府加大对治污的投入。买进先进的净化仪，加大对废水的处理量，可以

使重污染河长的比例减少；采取护岸、护坡措施，防治水土流失；开发低成本、高效率、易普及的废弃物处理设施，降低污物处理成本；改进垃圾填埋技术；发展生态农业。

（6）治污工程要持之以恒，要使经济的发展与环境相适应。政府在发展经济的

同时应该要兼顾环境因素，不要以牺牲环境为代价，盲目寻求经济越快越好发展，以免以后自食苦果。

以上是根据所建模型数据结果分析所得结论并结合长江实际情况对长江水质进行治理提出的几点意见和建议。

6.模型的评价、改进及推广

6.1模型的评价 优点：

（1）在求解第一问时，采用熵值法计算评价指标权重，避免了主观确定权重的人为因素，权重有强烈的现实意义。

（2）问题三采用的灰色预测模型具有少数据性、良好的时效性、较强的系统性和关联性等特征，弱化了预测对象的随机性变化，突出了趋势性，可以合理地对数据做出预测。

（3）在求解问题二时，我们把长江流水水质假设成是一维稳定河流水质，这样可大大方便求解上游给下游地区带来的污染物质量。 缺点：

（1）采用模糊综合评价时，实际中最常用的是最大隶属度原则，但用此分级标准评价环境污染程度存在一定的主观性。

（2）在求解问题三时，为了使问题得到方便的解决，采用水文年内全流域河长所占百分比进行求解，求出的结果与真实值会存在一定的偏差。 6.2模型的改进

（1）可以考虑固体废弃物的污染，每年长江汛期，垃圾随暴雨、洪水沿江而下，这些固体废弃物及白色污染也是相当严重的，因此评价预测污染情况应将其考虑在内。

（2）考虑生态系统逐年变坏，这必将会影响到长江的自净能力，随着污染物逐年累积增多，生态越来越差，处理污染的能力有可能出现指数下降，从而影响到未来污染的预测值。 6.3模型的推广

（1）模型一中，除了用模糊算法评价外，还可以用BP神经网络来求长江水质的在综合评价；而在求解问题三和问题四时采用的灰色预测，还可以用回归分析

预测、时间序列预测和神经网络预测等代替。

（2）建立的评价模型方法既适应于水环境质量评价，对其他一些环境质量评价也同样适用，如大气环境质量评价、富营养化程度评价等，特别是对参数和等级较多的综合评价，只需在标准样本训练时改变输入节点数和隐层节点数即可。因此，该方法有较广的适用范围。

7.参考文献

[1]. 郭显光.改进的熵值法及其在经济效益评价中的应用[J]. 系统工程理论与实

践.1988（12）：98－102

[2].蔡锁章.数学建模原理与方法.北京：海洋出版社。

[3].岂兴明等.MATLAB程序设计快速入门.北京：人民邮电出版社，2009。 [4].百度百科－熵值法：http://baike.http://www.32336.cn//view/2523301.htm [5].司守奎，数学建模算法与程序，烟台：海军航空工程学院，2007

8.附录

问题一 1－1

clear all,clc

m=xlsread('data2'); %?ó?a×?o??à2a??±ê′ú?? a=zeros(1,3); a=zeros(1,3); b=10*ones(1,3); c=zeros(17,3); sumc=zeros(1,3); f=zeros(17,3); f1=zeros(17,3); f2=zeros(17,3);

H=zeros(1,3); W=zeros(1,3); for i=1:3 for j=1:17 if a(1,i)<=m(j,i) a(1,i)=m(j,i); end

if b(1,i)>=m(j,i) b(1,i)=m(j,i); end end end f=a; y=b;

z=f-y; 2? z1=y-f; for i=1:17

c(i,1)=(m(i,1)-y(1,1))/z(1,1); c(i,2)=(m(i,2)-f(1,2))/z1(1,2);

c(i,3)=(m(i,3)-f(1,3))/z1(1,3);

%?¨ò?ì????ó %?¨ò?è¨???ó %?ó?a?à????±ê?D?????óμ?áD×?′óó?×?D??μ?°??%?ù?Yê?1-11

end

c %?ó1éò??ˉ?D?????ó for i=1:3 for j=1:16

sumc(1,i)=c(j,i)+c(j+1,i); end end sumc for i=1:17 for j=1:3

f(i,j)=(1+c(i,j))/(sumc(1,j)+17); %?ù?Yê?1-13μ?????1-14 f1(i,j)=log(f(i,j))/log(exp(1)); f2(i,j)=f(i,j).*f1(i,j) end end f f1 for i=1:3 for j=1:17

H(1,i)=-log(exp(1))/log(17)*(sum(f2(:,i))); %?ù?Yê?1-12 end end

H%?ó?aì????ó for i=1:3

W(1,i)=(1-H(1,i))/(3-sum(H(1,:))) %?ù?Yê?×ó1-15 end W

T=sum(W(1,:)) %?ó?aè¨???ó for i=1:17

m(i,1)=m(i,1)-5; %?ù?Yê?×ó1-16 m(i,2)=6-m(i,2); m(i,3)=1-m(i,3); end a=m for i=1:16

a(i+1,1)=a(i,1)+a(i+1,1); %?ù?Yê?×ó1-17 a(i+1,2)=a(i,2)+a(i+1,2); a(i+1,3)=a(i,3)+a(i+1,3); end f=a(17,:)

Q=sum(W.*a(17,:)) %?ó×?o???±ê￡??ù?Yê?×ó1-18

1－2

%?ó?a3ì?è?μêye clear all,clc; Q=xlsread('Q.xls'); qs=29.5733; Qs=repmat(qs,[28,1]); q=Q-Qs; [k,r]=size(q); m=zeros(1,k); n=zeros(1,k); for i=1:k if q(i,1)>0 m(1,i)=q(i,1); else m(1,i)=0;

end if q(i,1)<0 n(1,i)=q(i,1); else n(1,i)=0;

end end

m q1=sum(m) q2=sum(n) q3=sum(abs(q)) e=(q1+q2)/q3

%??·Y?D±ê×?2?′óóú0μ?oí %??·Y?D±ê×?2?D?óú0μ?oí

1－3

clear all; clc;

m=xlsread('data0.xls'); u=zeros(1,3);

u(1,1)=max(m(:,1)); % èü?a???aê?ò?Dí??±ê u(1,2)=min(m(:,2)); % ???ì?á????êy?a3é±?Dí??±ê u(1,3)=min(m(:,3)); % °±μa?a3é±?Dí??±ê U=u; v=zeros(1,3);

v(1,1)=min(m(:,1)); % èü?a???aê?ò?Dí??±ê v(1,2)=max(m(:,2)); % ???ì?á????êy?a3é±?Dí??±ê v(1,3)=max(m(:,3)); % °±μa?a3é±?Dí??±ê V=v; for i=1:28

ak(i,1)=m(i,1)/(U(1,1)+V(1,1)); % ?ù?Yê?1-21,?ó?aakr ak(i,2)=(U(1,2)*V(1,2))/((U(1,2)+V(1,2))*m(i,2)); % ?ù?Yê?1-23,?ó?aakr ak(i,3)=(U(1,3)*V(1,3))/((U(1,3)+V(1,3))*m(i,3)); % ?ù?Yê?1-23,?ó?aakr ck(i,1)=(U(1,1)*V(1,1))/((U(1,1)+V(1,1))*m(i,1)); % ?ù?Yê?1-22,?ó?ackr ck(i,2)=m(i,2)/(U(1,2)+V(1,2)); % ?ù?Yê?1-24,?ó?ackr ck(i,3)=m(i,3)/(U(1,3)+V(1,3)); % ?ù?Yê?1-24,?ó?ackr end

for i=1:28

m(i,1)=sum(ak(i,:)); % ?ù?Yê?1-25,?ó?aak,m?′ak n(i,1)=sum(ck(i,:)); % ?ù?Yê?1-26,?ó?ack,n?′ck l(i,1)=m(i,1)/(m(i,1)+n(i,1)); % ?ù?Yê?1-27,?ó?ark,l?′ck end m;n;f=l [a,b]=sort(f)

问题二 2－1

clear all,clc syms L;format short;

v0=xlsread('data00.xls'); for i=1:1:13 for j=1:6

v(i,j)=0.5*(v0(i,j)+v0(i,j+1));%?ó?a???ùá÷?ù???ó end end

L=[950 778 395 500 164 464]'; for i=1:13 for j=1:6

t(i,j)=L(j)/(v(i,j)*3600*24/1000);%?ó?aê±?????ó end end for i=1:13 for j=1:6

k1(i,j)=exp(-0.25*t(i,j));%2D′??μêy k2(i,j)=exp(-0.40*t(i,j));%2D′??μêy end end

V=xlsread('2liuliang.xls'); C1=xlsread('gao nongdu.xls'); C2=xlsread('an nongdu.xls'); for i=1:13 for j=1:6

m1(i,j)=C1(i,j)*V(i,j)*k1(i,j)/1000;%??è????êá?±??ˉá? m2(i,j)=C2(i,j)*V(i,j)*k2(i,j)/1000;%??è????êá?±??ˉá? product1(i,1)=C1(i,1)*V(i,1)/1000; product2(i,1)=C2(i,1)*V(i,1)/1000;

product1(i,j+1)=C1(i,j+1)*V(i,j+1)/1000-m1(i,j);%??è??′2ú??á? product2(i,j+1)=C2(i,j+1)*V(i,j+1)/1000-m2(i,j);%??è??′2ú??á? end end

问题三 2－1

clear all,clc;

t0=[174,179,183, 189, 207, 234, 205, 256, 270, 285]'; n=length(t0);

lamda=t0(1:n-1)./t0(2:n) %??????±è t1=cumsum(t0); %à??ó????

B=[-0.5*(t1(1:end-1)+t1(2:end)),ones(n-1,1)];Y=t0(2:end); r=B\\Y;

y=dsolve('Dy+a*y=b','y(0)=y0'); y=subs(y,{'a','b','y0'},{r(1),r(2),t1(1)}); yuce1=subs(y,'t',[0:n+9]);

%?aìá???¤2a???è￡??è?????¤2a?μ￡??ù??ê??￠·?·?3ìμ??a y=vpa(y,6) %???Dμ?6 ±íê???ê?6 ??êy×? yuce=diff(yuce1); %×÷2?·?????￡???DDêy?Y?1?- yuce=[t0(1),yuce]; yuce_old=yuce(1:n);

yuce_new=yuce(n+1:end); %?óμ?μ?ê????¤2a?μ epsilon=t0'-yuce_old ;%????2D2? delta=abs(epsilon./t0') %?????à???ó2?

rho=1-(1-0.5*r(1))/(1+0.5*r(1))*lamda' %??????±è??2??μ yuce_new=yuce_new';

2－2

%ò?àà?? clear all,clc; t=xlsread('yin.xls'); t01=t(:,1); n=length(t01);

lamda=t01(1:n-1)./t01(2:n); %??????±è t11=cumsum(t01); %à??ó????

B1=[-0.5*(t11(1:end-1)+t11(2:end)),ones(n-1,1)];Y1=t01(2:end); r1=B1\\Y1;

y1=dsolve('Dy+a*y=b','y(0)=y0'); y1=subs(y1,{'a','b','y0'},{r1(1),r1(2),t11(1)}); yuce11=subs(y1,'t',[0:n+9]);

%?aìá???¤2a???è￡??è?????¤2a?μ￡??ù??ê??￠·?·?3ìμ??a

y1=vpa(y1,6); %???Dμ?6 ±íê???ê?6 ??êy×? yuce1=diff(yuce11); %×÷2?·?????￡???DDêy?Y?1?- yuce1=[t01(1),yuce1];

yuce_new1=yuce1(n+1:end); %?óμ?μ?ê????¤2a?μ result1=yuce_new1'; %?tàà?? t02=t(:,2); n=length(t02);

lamda=t02(1:n-1)./t02(2:n); %??????±è t12=cumsum(t02); %à??ó????

B2=[-0.5*(t12(1:end-1)+t12(2:end)),ones(n-1,1)];Y2=t02(2:end); r2=B2\\Y2;

y2=dsolve('Dy+a*y=b','y(0)=y0'); y2=subs(y2,{'a','b','y0'},{r2(1),r2(2),t12(1)}); yuce12=subs(y2,'t',[0:n+9]);

%?aìá???¤2a???è￡??è?????¤2a?μ￡??ù??ê??￠·?·?3ìμ??a y2=vpa(y2,6); %???Dμ?6 ±íê???ê?6 ??êy×? yuce2=diff(yuce12); %×÷2?·?????￡???DDêy?Y?1?- yuce2=[t02(1),yuce2];

yuce_new2=yuce2(n+1:end); %?óμ?μ?ê????¤2a?μ result2=yuce_new2'; %èyàà?? t03=t(:,3); n=length(t03);

lamda=t03(1:n-1)./t03(2:n); %??????±è t13=cumsum(t03); %à??ó????

B3=[-0.5*(t13(1:end-1)+t13(2:end)),ones(n-1,1)];Y3=t03(2:end); r3=B3\\Y3;

y3=dsolve('Dy+a*y=b','y(0)=y0'); y3=subs(y3,{'a','b','y0'},{r3(1),r3(2),t13(1)}); yuce13=subs(y3,'t',[0:n+9]);

%?aìá???¤2a???è￡??è?????¤2a?μ￡??ù??ê??￠·?·?3ìμ??a y3=vpa(y3,6);%???Dμ?6 ±íê???ê?6 ??êy×? yuce3=diff(yuce13); %×÷2?·?????￡???DDêy?Y?1?- yuce3=[t03(1),yuce3];

yuce_new3=yuce3(n+1:end); %?óμ?μ?ê????¤2a?μ result3=yuce_new3'; %??àà?? t04=t(:,4); n=length(t04);

lamda=t04(1:n-1)./t04(2:n); %??????±è t14=cumsum(t04); %à??ó????

B4=[-0.5*(t14(1:end-1)+t14(2:end)),ones(n-1,1)];Y4=t04(2:end); r4=B4\\Y4;

y4=dsolve('Dy+a*y=b','y(0)=y0'); y4=subs(y4,{'a','b','y0'},{r4(1),r4(2),t14(1)}); yuce14=subs(y4,'t',[0:n+9]);

%?aìá???¤2a???è￡??è?????¤2a?μ￡??ù??ê??￠·?·?3ìμ??a y4=vpa(y4,6); %???Dμ?6 ±íê???ê?6 ??êy×? yuce4=diff(yuce14); %×÷2?·?????￡???DDêy?Y?1?- yuce4=[t04(1),yuce4];

yuce_new4=yuce4(n+1:end); %?óμ?μ?ê????¤2a?μ result4=yuce_new4'; %??àà?? t05=t(:,5); n=length(t05);

lamda=t05(1:n-1)./t05(2:n); %??????±è t15=cumsum(t05); %à??ó????

B5=[-0.5*(t15(1:end-1)+t15(2:end)),ones(n-1,1)];Y5=t05(2:end); r5=B5\\Y5;

y5=dsolve('Dy+a*y=b','y(0)=y0'); y5=subs(y5,{'a','b','y0'},{r5(1),r5(2),t15(1)}); yuce15=subs(y5,'t',[0:n+9]);

%?aìá???¤2a???è￡??è?????¤2a?μ￡??ù??ê??￠·?·?3ìμ??a y5=vpa(y5,6); %???Dμ?6 ±íê???ê?6 ??êy×? yuce5=diff(yuce15); %×÷2?·?????￡???DDêy?Y?1?- yuce5=[t05(1),yuce5];

yuce_new5=yuce5(n+1:end); %?óμ?μ?ê????¤2a?μ result5=yuce_new5'; %áó??àà??

t06=t(:,6); n=length(t06);

lamda=t06(1:n-1)./t06(2:n); %??????±è t16=cumsum(t06); %à??ó????

B6=[-0.5*(t16(1:end-1)+t16(2:end)),ones(n-1,1)];Y6=t06(2:end); r6=B6\\Y6;

y6=dsolve('Dy+a*y=b','y(0)=y0'); y6=subs(y6,{'a','b','y0'},{r6(1),r6(2),t16(1)}); yuce16=subs(y6,'t',[0:n+9]);

%?aìá???¤2a???è￡??è?????¤2a?μ￡??ù??ê??￠·?·?3ìμ??a y6=vpa(y6,6); %???Dμ?6 ±íê???ê?6 ??êy×? yuce6=diff(yuce16); %×÷2?·?????￡???DDêy?Y?1?- yuce6=[t06(1),yuce6];

yuce_new6=yuce6(n+1:end); %?óμ?μ?ê????¤2a?μ result6=yuce_new6';

result=[result1,result2,result3,result4,result5,result6]; m=result m

2－3

clear,clc; a=xlsread('yin4.xls'); x=a'; s=sum(x); for j=1:10 for i=1:6

x(i,j)=x(i,j)/s(j); end end c=x.'*100