《数学建模优秀论文北京水资源短缺风险综合评价.doc》由会员分享，可在线阅读，更多相关《数学建模优秀论文北京水资源短缺风险综合评价.doc（18页珍藏版）》请在三一办公上搜索。
1、北京市水资源短缺风险综合评价摘要本文探讨的是北京市水资源短缺风险的综合评价及预测问题。水资源短缺已成为目前大多数城市都面临的严峻问题，如何对水资源风险的主要因子进行识别，对风险造成的危害等级进行划分，对不同风险因子采取相应的有效措施规避风险或减少其造成的危害，这对社会经济的稳定、可持续发展战略的实施具有重要的意义。首先，我们利用主成分分析法定量分析和筛选水资源短缺风险的影响因子；然后，建立熵权模糊综合评价模型，计算水资源短缺风险值并作出风险等级划分；最后，用灰色系统模型预测出北京市2011年和2012年的水资源短缺风险值并判断所属风险等级。第一，我们根据1979年2010年9个北京市水资源短缺
2、风险影响因子：降水量、地下水、污水处理率、GDP、城市生活用水、工业用水、农业用水、人口和森林覆盖率的相关数据，利用主成分分析法，通过使用SPSS软件最终筛选出人口、污水处理率、GDP、森林覆盖率和降水量5个风险指标。第二，利用第一问中筛选出的5个风险指标，通过对数据的分析类比，将信息论中的熵值理论运用到确定风险指标的权重，建立了基于熵权的水资源短缺风险模糊综合评价模型,通过matlab,excel等软件工具计算出了2001年2010年北京市水资源短缺风险值，并进行了风险等级划分。具体各年份水资源短缺风险值及风险等级见表5。第三，根据第二问中解出的2001年2010年的水资源短缺风险值，建立了
3、灰色系统预测模型，预测出了2011年和2012年北京水资源短缺的风险值分别为0.7465、0.7839，所属风险等级分别为中等、较低。第四，根据第一问中确定的五个影响力较大的风险因子，通过第二问的模型定量求解出每年的水资源短缺风险值和第三问中预测出的未来两年水资源短缺风险值以及相应风险等级，向北京市水行政主管部门写了一份关于规避主要风险因子的建议书。本文亮点：水资源短缺风险评价中各指标具有模糊性和不确定性，将信息论中熵值理论运用于水资源短缺风险评价中，建立了基于熵权的水资源短缺风险模糊综合评价模型，有效地解决这个问题。关键词：水资源短缺风险评价 主成分分析法 熵权模糊综合评价模型 灰色系统预测
4、模型1问题重述近年来，我国、特别是北方地区水资源短缺问题日趋严重，水资源成为焦点话题。北京就是我国严重缺水的都市之一，北京市水资源短缺已经成为影响和制约首都社会和经济发展的主要因素。虽然政府采取了措施,但是，气候变化和经济社会不断发展，水资源短缺风险始终存在。如何对水资源风险的主要因子进行识别，对风险造成的危害等级进行划分，对不同风险因子采取相应的有效措施规避风险或减少其造成的危害，这对社会经济的稳定、可持续发展战略的实施具有重要的意义。北京2009统计年鉴及市政统计资料提供了北京市水资源的有关信息。利用这些资料和你自己可获得的其他资料，讨论以下问题：1 评价判定北京市水资源缺风险的主要风险因
5、子是什么？影响水资源的因素很多,例如：气候条件、水利工程设施、工业污染、农业用水、管理制度，人口规模等。2建立一个数学模型对北京市水资源短缺风险进行综合评价，作出风险等级划分并陈述理由。对主要风险因子,如何进行调控，使得风险降低？ 3 对北京市未来两年水资源的短缺风险进行预测，并提出应对措施。4 以北京市水行政主管部门为报告对象，写一份建议报告。2.模型的基本假设1、假设筛选出的风险因子之间是相互独立的；2、假设每一个风险因子都是可以等价量化；3、不考虑其他微弱风险因子对北京水资源短缺风险的影响。3.主要变量符号说明 特征根； 累计贡献率； 水资源短缺风险因素论域； 等级论域； 隶属度； 熵权
6、； 综合评价风险值； 相关矩阵； 模糊矩阵。注：以上变量量纲均为1。4.问题分析北京人均水资源占有量不足300m3，为全国人均的1/8，世界人均的1/30，属重度缺水地区。本文是定量评估北京市水资源短缺风险的问题。首先，本文通过系统定性分析，列出可能加大水资源风险的各方面因素。然后定量分析各风险因素对水资源短缺风险的影响程度，从而确定主要风险因子。首先，运用主成分分析对高维变量系统进行最佳综合与简化，同时客观确定各个指标的权重，从而筛选出权重较大的指标，确定敏感因子。（1）降水量（2）地下水（3）污水处理率（4）GDP风险指标体系 （5）城市生活用水（6）工业用水（7）农业用水（8）人口（9）
7、森林覆盖率然后，第二个问题的概念是相对比较模糊的，所以我们选择基于熵权的模糊综合评价法，对每个风险因子对每年缺水造成的风险度进行定量运算，然后根据具体数值大小确定本年的风险等级。最后，要根据第二问的结果对未来两年的缺水风险度进行预测，选择灰色预测模型，定量预测风险值。5.模型的建立与求解5.1 水资源短缺主要风险因子的确定5.1.1主成分分析模型的建立与求解数据的标准化处理其中为第i个分区第j个指标的值，为第j个指标的样本均值和样本标准差。计算数据表的相关矩阵S。求S 的J 个特征值:,以及对应的特征向量,它们标准正交称为主轴。求主成分: 精度分析:通过求累计贡献率E 来判断, ,一般要求取E
8、 85 %的最小m值,则可得主平面的维数m,从而可对m个主成分进行综合分析。在获得特征向量与特征值,并确定主超平面的维数之后,可以计算主因子载荷矩阵,其计算公式为: 其中根据北京历年的9个风险因子的主成分分析，得到如下结果：水资源短缺风险：计算所得的前两个特征根所以取m=2。根据附录一数据，利用spss软件求的风险因子之间的相关系数矩阵以及初始因子载荷矩阵：表1 相关系数矩阵（Correlation Matrix）降水量地下水污水处理率（%）GDP降水量10.298324198-0.45515-0.391129413地下水0.2983241981-0.83512-0.836378634污水处理
9、率（%）-0.455149554-0.83512148310.979858457GDP-0.391129413-0.8363786340.9798581城市生活用水-0.363007159-0.6914272330.7840880.7996124工业用水0.4608652480.826516422-0.92516-0.903280316农业用水0.0227758480.612722184-0.63506-0.643773483人口-0.389440799-0.8211613570.9494090.958830515森林覆盖率-0.371723509-0.7740657240.8627670.8
10、68394066续表1城市生活用水工业用水农业用水人口森林覆盖率降水量-0.3630071590.4608650.022776-0.38944-0.37172地下水-0.6914272330.8265160.612722-0.82116-0.77407污水处理率（%）0.784087852-0.92516-0.635060.9494090.862767GDP0.7996124-0.90328-0.643770.9588310.868394城市生活用水1-0.73369-0.688460.9105260.959874工业用水-0.73368766510.493078-0.88701-0.8092
11、2农业用水-0.6884560470.4930781-0.72184-0.70412人口0.910525659-0.88701-0.7218410.960056森林覆盖率0.959874436-0.80922-0.704120.9600561表2 方差分解主成分提取分析表（Total Variance Explained）ComponentInitial Eigenvalues% ofVarianceCumulative %Extraction Sums of Squared LoadingsTotalTotal% of VarianceCumulative %16.85402576.1558
12、376.155836.854076.1558376.1558321.02304911.3672087.52301.023011.3672087.523030.4894215.43801092.9610440.2960753.28972496.2507750.1993262.21473698.4655060.0847050.94116599.4066770.0304320.33813199.7448080.0159040.17671699.9215290.0070630.078477100在获得特征值和特征向量后，计算正交旋转后的主因子载荷矩阵,结果如表3所示。表3 初始因子载荷矩阵（Compo
13、nent Matrix）Component12降水量-0.440240.831401221地下水-0.87244-0.054546952污水处理率（%）0.963287-0.082689493GDP0.961511-0.01886892城市生活用水0.8997670.090408486工业用水-0.913260.177236703农业用水-0.72227-0.52348974人口0.9886820.046216029森林覆盖率0.950710.076731512由初始因子载荷矩阵可知：主成分（1）中各因子载荷值，从正方向看，比较大的是人口，污水处理率，GDP,森林覆盖率，分别为0.988682
14、，0.963287，0.961511和0.95071。主成分（2）中各因子载荷值，从正方向看，唯一比较大的就是降水量，为0.831401221。根据以上分析可知：在水资源短缺风险中载荷较大的指标为人口，污水处理率，GDP,森林覆盖率和降水量。5.2问题二模型的建立与求解5.2.1熵权模糊综合模型的建立（1）建立影响水资源短缺的因素域（2）建立等级域（3）在水资源短缺的因素域U与等级域V之间进行隶属度分析，建立模糊关系矩阵 （1）矩阵中：表示因素域U中第i个因素对于等级域中第j个等级的隶属度。（4）模糊综合评价模型及其值的计算 水资源短缺风险模糊综合评价B等于W与R两个矩阵的乘积，即 （2）其中
15、：为两个因素对水资源短缺风险指标的权重，且满足，因此B的水资源风险评价结果等于,就是最终的评价结果，也就是此对象具体属于那个风险等级。5.2.2建立评价集水资源的短缺风险大小是相对且模糊的，不可能定性描述，是属于模糊集理论。那么根据模糊数学理论1，可以直接把主成分分析中提取的五个影响因素，根据指标大小不同分为五个等级,，分别为低风险，较低风险，中等风险，较高风险，高风险。根据文献3,4的分级标准，我们可以把个指标的分级情况表示如下：表4：等级划分指标级级级级级森林覆盖率指数0.430.420.40.380.3GDP指数0.30.40.50.70.9污水处理指数0.70.60.50.50.4人口
16、指数0.140.150.160.170.18降雨量指数0.60.40.40.350.35.2.3建立隶属函数水资源缺乏风险是一个模糊概念，它分级标准也是模糊的，用隶属度来刻画分级界线较为合理。以各指标的五级标准，作个级别的隶属函数，如第一个指标的隶属函数： （3） 以此类推，可得GDP，污水处理，人口和降雨量的隶属函数。这样，通过5个指标，5级标准，将组成5行5列的隶属函数矩阵。 （4） 5.2.4用熵权法确定权重系数 建立10个评价对象，5个评价指标的判断矩阵R 把判断矩阵R通过归一化公式进行归一，得到矩阵B （5）公式中： 和 分别为同一评价指标下的对象评价值中的最大和最小值。参考文献4可
17、以有发现关于熵权法的公式,根据其中熵的定义公式，计算10个评价对象5个评价指标的熵值 （6）公式中： ，其中但是当时无意义，因此需要对加以修正，定义为 （7）根据熵值计算评价指标的熵权 （8） （9）公式中，i=1,2,3,4,5，且满足 至此就算出每个指标的权重，虽然计算权重的方法很多，但使用熵权法避免了主观权重的误差性，计算的权重值更具有代表意义。5.2.5水资源短缺风险综合评价（1）2001-2010年水资源短缺指标的取值 表5：水资源短缺指标年份降雨量污水处理率GDP人口森林覆盖率20010.340.420.370.140.3920020.370.450.430.140.4120030
18、.440.500.500.150.4120040.480.540.610.150.4220050.410.620.690.150.4220060.320.740.790.160.4320070.480.760.940.160.4320080.630.791.050.170.4420090.340.771.060.170.4420100.460.791.100.180.44（2）权重系数的确定 用excel根据公式（5）归一化判断矩阵B 由公式（6）和（7）可得评价指标的熵值：H=（0.9780,0.9455,0.9882,0.9451,0.9303）再通过公式（9）计算出各评价指标的权重W=（
19、0.2328,0.2238,0.1996,0.1490,0.1948）。（3）每年综合评价值的求解模糊矩阵的求解由表五的数据，用excel根据公式（3）中的隶属函数求出隶属度，得到模糊矩阵 通过matlab由式（2）可得北京市水资源风险综合评价为 同理可求出2002至2009年的综合评价值，如下表：表6 2001-2010水资源短缺风险评价结果年份v1(低）v2(较低）v3(中等）v4(较高）v5(高）等级20010.34250.66520.45760.56640.7425高20020.24790.39780.44320.74450.6934较高20030.35420.37560.47820.
20、78760.6543较高20040.34230.45760.78950.48760.6787中等20050.46740.65470.78990.69870.5647中等20060.82210.89340.67880.65640.2341较低20070.47630.32460.89340.37540.3212中等20080.76450.54260.47340.73340.4354低20090.67580.78340.63430.65330.2136较低20100.87460.79870.43230.38560.6576低（4）结果分析 由上面得出的表三中的结果可以看出，2000年以后缺水风险值基
21、本处于降低趋势，特别是2008年以后，由于有奥运会的影响，很明显缺水风险很低。与此同时，对于分析出来的主要风险因子，主要是降雨量和各种用水量在影响，那么我们的调控措施就主要在用水控制与循环利用上，可适当提高水价和提高污水处理率，这样可使得风险降低。5.3水资源短缺风险的预测5.3.1灰色预测模型的建立我们假设一组原始序列 引入一阶弱化算子D，令 其中 那么的1AGO为对作紧邻均值生成从而可以得到矩阵和矩阵 对参数进行最小二乘估计，得到估计参数 和。则GM(1,1)白化方程为 相应的时间式为（预测模型） 5.3.2预测模型的求解首先我们可以得到北京市2001年2010年水资源短缺风险值的原始序列
22、引入一阶弱化算子D，令其中从而得到 那么的1AGO为 对作紧邻均值生成 在此，我们使用MATLAB编程实现对的求解（具体程序见附录二）。于是可以得到于是我们可以得到矩阵和矩阵 对参数进行最小二乘估计，采用matlab编程完成得到结果（具体程序见附录二）于是我们就可以得到估计参数：。则GM(1,1)白化方程为：相应的时间式为（预测模型）那么通过这个预测模型，借助MATLAB软件就可以重新得到2001年2010年北京市水资源短缺风险值序列（具体程序见附录二）：残差序列为相对误差序列平均相对误差 精度为二级 精度为一级所以 , 均方差比值为一级。计算小误差概率：所以，小误差概率为一级。通过上面的检验
23、分析，该模型所有的检验都合格，故可用来进行预测，2011、2012年的预测值为：得出结论：2011年北京市水资源短缺风险等级为中等2012年北京市水资源短缺风险等级为较低5.4问题四的解决方案给北京市水行政主管部门的建议报告本文通过系统分析,结合定性与定量研究结果,得出北京水资源短缺的主要风险敏感因子为: 人口，污水处理率，GDP,森林覆盖率和降水量。找到敏感风险因子后,也就找到了风险防范的突破口。我们把以上敏感因子分为两类:一类为不可改变因素:降水量;另一类为可改变因素:人口，污水处理率，GDP，森林覆盖率。不可变因素是生态系统长期作用的结果,非短期内人类所能改变,而可改变因素却完全可以通过
24、人为的自身调控在短期内取得效果。因此我们应重点对可改变的敏感因子进行调控,以期最大限度的降低水资源短缺风险。由此依据所得出的风险敏感因子,我们提出以下四条降低流域水资源短缺风险的措施和建议:(1) 均衡人口流动方向，在有越来越多人涌向北京趋势下，而环境承载力有限。北京水利部门可以与有关部门协商，向中央提交有关报告，大力支持其他地区建设，均衡人口密度。(2) 大幅度地提高污水处理率能减少水环境污染,从而改善水质,从质上降低水资源短缺风险。努力实现对工业和城市生活排放的污水进行100 %的处理,处理后的污水应至少能达到农业灌溉要求,能够大量回用。在建立污水厂的同时,就应该对污水回用工程进行规划布置
25、,尽量将水的重复利用率提高到最高水平。(3) 可持续发展，在追求GDP的增长的同时，考虑资源承受能力，适当提高水价，从而限制不必要的水资源的浪费，提高水资源利用率。(4) 以水养林，以林护水。从水资源收益中的一部分拿来支持植树造林，森林覆盖率高了，对于水土保持、降低水资源短缺风险有着显著效果。6.模型的改进与评价6.1模型优点（1）本文建立了主成分分析模型，消除了降水量、GDP、工业用水量等多个指标间的相关影响，减少了指标选择的工作量，更加客观的反映了风险因子与缺水量之间的显示关系。（2）本文建立了基于熵权的水资源综合评价模型，避免传统层次分析法中两两对比构造判断矩阵易于发生目标先后不一致性问
26、题。引入信息论中的熵值理论从数据本身所反映的信息的无序化效应值来计算权重系数使得权 重的确定有了一定的理论依据减少了评价中人为主观因素的影响和评价结果的主观性。（3）本建立了灰色预测模型，进行预测北京2011、2012两年水资源短缺风险，切实量化了北京未来两年的水资源短缺风险。（4）本文所有模型，整个过程思路清晰、计算简单、评价结果客观。6.2模型缺点(1)主成份分析模型的不可避免的要损失一部分信息量。（2）灰色模型进行预测，鉴于序列的光滑性不好，背景值的选取比较困难。（3）由于水资源受环境以及人为影响较大，所以预测结果难免有误差。6.3模型改进方向由于数据的光滑度不是很好，灰色模型在进行预测
27、时难免会带来一定的误差，基于正弦函数变换的改进GM(1.1)模型，在已建立的模型上提高光滑度，从而提高预测精度。设非负原始数据虚序列为递增序列： 将序列进行标准处理，使得对任意的，标准化后的均在区间内。 展开得到： 其中，为背景值为：记则得待辨识向量的最小乘解为：白化微分方程的离散解为： 还原到原始数据为： 最后将标准化后的还原为。6.4模型的推广（1）基于熵权的水资源短缺风险模糊评价模型可以对北京以外的城市进行水资源风险预测，或者对其他能源缺乏的风险预测。也可以应用到国家战略水资源调配、储备。（2）灰色模型也可以应用到其他自然资源的风险预测。7.参考文献1 衷平，沈珍瑶等，石羊河流域水资源短
28、缺风险敏感因子的确定，干旱区资源与环境，19（2）：81-86，2005。2 韩宇平，阮本清，解建仓，水资源系统风险评估研究，西安理工大学学报，19（1）：41-45，2003。3 刘涛，邵东国，水资源系统风险评估方法研究，武汉大学学报工学版，38（6）：66-77，2005。4 黄明聪，解建仓，阮本清等；基于支持向量机的水资源短缺风险评价模型及其应用，水利学报，33（3）：255-299，2007。5 阮本清，韩宇平，王浩等，水资源短缺的模糊综合评价，水利学报，36（8）：906-912，2005。6 北京统计局，北京统计年鉴2009，/2009-tjnj/，2011年5月7日7 姜启源，数
29、学模型（第三版） M，北京：高等教育出版社， 2003。8.附录8.1附录一 各风险因子原始数据表7 各风险因子原始数据年份降雨量（毫米）地下水污水处理率（%）GDP城市生活用水工业用水农业用水（亿立方米）人口（万人）森林覆盖率1979718.426.910.2120.14.3714.3724.18897.122.31980380.727.19.4139.14.9413.7731.83904.320.11981393.227.910.8139.24.312.2131.6919.220.11982544.427.610.9154.94.5213.8928.8193520.11983489.927
30、.410.2183.14.7211.2431.695020.11984488.827.810216.64.01714.37621.8496520.1198572126.810257.14.3917.210.1298122.11986665.325.98.9284.97.189.9119.46102822.861987683.927.67.7326.87.2614.019.68104722.91988673.3277.4410.26.414.0421.991061251989442.226.86.64566.4513.7724.421075261990697.326.97.3500.87.041
31、2.3421.741086281991747.927.36.6598.97.4311.922.7109428.431992541.527.41.2709.110.9815.5119.94110230.331993506.727.23.1886.29.5915.2820.35111231.331994813.227.19.61145.310.3714.5720.93112532.391995572.526.919.41507.711.7713.7819.331251.132.681996700.926.721.21789.29.311.7618.951259.433.241997430.925.
32、8222075.611.111.118.12124034.221998731.725.922.5237612.210.8417.391245.635.61999266.926.2252677.612.710.5618.451257.236.32000371.127.139.4316113.3910.5216.491363.636.52001338.927.2423710.5129.217.41385.138.782002370.424.2454330.410.87.515.51423.240.572003444.925.450.15023.8138.413.81456.440.87200448
33、3.526.853.96060.312.87.713.51492.741.912005410.724.962.46886.313.46.813.2153842200631824.373.8786113.76.212.8158142.52007483.924.376.29353.313.95.812.41633432008626.324.278.910488145.212169543.52009338.924.577.21056614.25.113169943.62010456.424.978.51100914.14.912.6175543.78.2附录二 灰色预测模型程序代码（1）clc%求对
34、作紧邻均值生成代码%x=;z(1)=x(1);for i=25z(i)=0.5*(x(i)+x(i-1);endformat long gz（2）B=;X=ones(4,1) B;Y=;b,bint,r,rint,stats=regress(Y,X)8.3附录三 归一化程序代码d=338.9,42,3710.5,1385.1,38.78;370.4,45,4330.4,1423.2,40.57; 444.9,50.1,5023.8,1456.4,40.87;483.5,53.9,6060.3,1492.7,41.91;410.7,62.4,6886.3,1538,42,318;73.8,7861,1581,42.5,483.9;76.2,9353.3,1633,43,626.3;78.9,10488,1695,43.5,338.9;77.2,10566,1699,43.6,456.4;78.5,11009,1755,43.7,436.4; N,M=size(a);b=ones(1,M); %1矩阵c=b;d=c;for i=1:M c(i)=prod(a(i,:);%每行元素乘积 d(i)=c(i)(1/M);%求M次方根 b(i)=c(i)/sum(c);%归一化处理endb%特征向量lamda=sum(a*b)./(N.*b)
