北京市水资源短缺风险综合评价89481

PAGE10北京市水资源短缺风险综合评价西南交通大学熊竣熙、周成龙、王义目录北京市水资源短缺风险综合评价……………………1摘要……………………21问题提出………………………31.1问题背景………………31.2问题的提出……………32问题分析………………………32.1概念的解释……………32.2问题的具体分析………33模型假设………………………44符号约定………………………45模型的建立与求解……………55.0数据的预处理…………55.1问题一：确定北京水资源短缺的主要影响因子……55.1.1指标的选取……………………55.1.2使用主成分分析法确定主要因子……………65.1.3主成分分析结果………………75.2问题二：水资源短缺风险等级的划分………………85.2.1等级的划分……………………85.2.2构建风险潜在函数……………95.2.3历年风险等级归类……………95.2.4结果的合理性检验…………115.2.5相关的调控措施及建议……115.3问题三：北京市水资源短缺的预测………………125.3.1利用灰色模型预测北京市水资源短缺风险………………125.3.2预测结果分析………………145.3.3利用BP神经网络检验灰色模型的预测正确性……………145.4问题四：对北京市水行政主管部门的建议………166模型的推广与评价…………176.1模型的缺陷…………176.2模型的优点…………176.3模型的推广…………177参考文献……………………178附录…………17摘要：本文采用了主成分分析法、模糊数学、灰色模型预测等方法对北京市的水资源短缺风险进行了评价与预测。问题一中，为了找出北京市水资源短缺的主要风险因子，本文首先找出了11个影响水资源短缺的因子，构建了包括来水风险、用水风险、管理制度调节在内的指标体系。利用插值、均值替换法等方法对数据进行修补与校正，接着本文采用主成分分析法，求得累计贡献率达到86.69%的三个主成分，在每一个主成分中找出影响力较大的因素为：蓄水量、降雨量、地下水埋深和人口总量。问题二中，本文先将风险暂且分为5级，分别为{风险极限，高风险，中度风险，低风险，无风险}，相应的评分为{5,4,3,2,1}。接着选取偏大型柯西分布作为隶属函数，计算出各级风险的隶属度，然后本文构建了潜在风险值函数，运用查找到的数据代入风险潜在函数计算出每一年的风险潜在值；将数据极差归一化后得到每一年的风险隶属度，从而评判出每一年的风险等级。本文划分的结果为，前十年为低风险甚至是无风险年，中间十年在向中度风险过度，后十年为中度风险年。为了验证等级划分的合理型，本文利用matlab对30年的11项指标进行聚类，30年的数据大致被分为3类，从1979—1994年的数据为一类，1995—1999年的数据为一类，2000—2008年的数据为一类，与等级划分结果很相似，从而验证了本文方法的合理型。针对本小问的分析结果，本文提出了一些控制水资源短缺的建议（见正文）。问题三中，本文构建了灰色模型GM（1,1）对11项指标的2009和2010年数据进行预测。由于较早的数据不具有太大的代表性，于是本文取近20年的数据进行预测，将预测的数据标准化后代入风险潜在函数，得到2009和2010年的风险潜在值为，；将其极差归一化后得到其风险隶属度，其隶属风险等级为中度风险。为了检验预测的合理性，本文又采用BP神经网络预测，利用近20年潜在风险值，对2009和2010年度的水资源短缺的潜在风险值进行预测，预测结果与灰色系统很相似，并且风险程度同样为中度风险，从而验证了本文预测的合理性。问题四中，本文针对主要风险因子提出一些实用建议，以通过这些方法改善目前北京市的水资源短缺状况。关键词：主成分分析风险潜在函数模糊数学灰色模型GM（1,1）BP神经网络一问题提出1.1问题背景水资源，是指可供人类直接利用，能够不断更新的天然水体。主要包括陆地上的地表水和地下水。风险，是指某一特定危险情况发生的可能性和后果的组合。水资源短缺风险，泛指在特定的时空环境条件下，由于来水和用水两方面存在不确定性，使区域水资源系统发生供水短缺的可能性以及由此产生的损失。近年，我国、特别是北方水资源短缺问题日趋严重，水资源成为焦点话题。以北京市为例，北京是世界上水资源严重缺乏的大都市之一，其人均水资源占有量不足300m31.2问题的提出《北京2009统计年鉴》及市政统计资料提供了北京市水资源的有关信息。利用这些资料和你自己可获得的其他资料，讨论以下问题：评价判定北京市水资源短缺风险的主要风险因子是什么？影响水资源的因素很多,例如：气候条件、水利工程设施、工业污染、农业用水、管理制度，人口规模等。2建立一个数学模型对北京市水资源短缺风险进行综合评价，作出风险等级划分并陈述理由。对主要风险因子,如何进行调控，使得风险降低？3对北京市未来两年水资源的短缺风险进行预测，并提出应对措施。4以北京市水行政主管部门为报告对象，写一份建议报告。二问题分析2.1概念解释2.1.1风险潜在函数由于风险大小是一个模糊的概念，于是本文定义风险潜在函数将风险的大小进行量化。风险潜在函数通过风险的指标综合求得。2.2问题的具体分析2.2.1问题一的分析问题一要求评价判定北京市水资源短缺风险的主要风险因子。在本小问中无疑会有很多指标和大量的数据，若逐步分析每一个因子，工作量势必很大。其次，指标之间有一定的信息重叠程度，采用的方法必须避开这一重叠现象。于是本文决定采用主成份分析法进行分析，确定各个因子的贡献率，以累计贡献率达到85%-90%来确定主成份。2.2.2问题2的分析问题二要求对北京的水资源短缺做出综合评价并别划分等级，由于等级的划分是一个很模糊的概念，于是本文决定采用模糊数学中的隶属函数。不妨先将风险等级划分为5级，即｛风险极限，高风险，中度风险，低风险，无风险｝，相对应的数值为｛5，4，3，2，1｝。然后根据实际情况，本文采用偏大型柯西分布确定每一级风险的隶属度，再通过计算每年的风险潜在函数，由此将每一年的风险等级归类。由于等级划分的合理性未知，本文采用模糊聚类进行验证。利用matlab编程将11个指标的30年数据进行聚类，一共可以分为4大类。再将风险等级划分出的4类年份与模糊聚类确定的4类年份进行对比即可。2.2.3问题3的分析问题三要求对北京市未来两年的水资源短缺进行预测，本文使用灰色系统模型，将各项指标的未来两年发展情况进行了预测，然后利用预测出的2009和2010的数值带入风险潜在函数，计算出2009和2010年的水资源短缺风险。为了验证结果的合理性，本文又采用神经网络模型进行预测，比较两种预测的结果从而确定方法的合理性。三模型假设1.所获得的数据真实可信2.不会出现导致指标产生巨大变动的突发状况3.水资源短缺风险不受除与之相关统计指标以外其他因素的影响4.30年中大气中的水循环正常进行四符号约定——一致化以后第个指标第年的数据——标准化以后第个指标第年的数据——第个主成分——第个主成分的贡献率——第个主成分的累计贡献率——风险潜在函数——发展灰度——内生控制灰度——灰色模型级比——灰色模型级比界限——隐层神经元数目五模型的建立与求解5.0数据的预处理5.0.1异常数据的处理在不考虑自然灾害等突发状况的影响条件下，查找到的数据应该是平缓变化的，如果发现某个数据异常，则使用均值替换法将其替换：（1）5.0.2缺失数据的处理对于缺失的数据，本文直接调用matlab中的一维插值函数，选用样条函数进行修补。5.0.3数据的一致化处理在选取的11个指标中，部分指标是极大型指标，部分指标是极小型指标。对于风险函数，某些数据越大对风险函数的促进作用越大，而某些数据则相反。所以将极小型指标极大化处理（2）5.1问题一：确定北京水资源短缺的主要影响因子5.1.1指标的选取对水资源短缺有影响的指标有很多，大致有三类：1.自然状况；2.产业用水和生活用水；3.管理体制影响。在查阅相关知识后，最终选取了11个指标。潜在风险函数体系如图5-1所示：图5-1潜在风险函数体系5.1.2使用主成份分析法确定主要因子Step1：对原始数据的标准化处理问题一涉及到的评价指标一共有8个，评价对象为1979年——2008年。在数据的预处理完成后，将第年的第个指标记为，按照标准化方法进行标准化处理。（3）其中，，得到标准化矩阵，记为Step2：计算相关系数矩阵相关系数矩阵：（4）记第个指标与第个指标的相关系数为，则相关系数矩阵为，其中，，Step3：计算特征值与特征向量计算相关系数矩阵的特征值，及其对应的特征向量，其中，有特征向量组成11个新的指标变量：（5）其中，为第主成分，Step4：选择个主成分计算特征值的信息贡献率和累计贡献率。称（6）为主成分的信息贡献率；称（7）为的累积贡献率。当时，可以选取前个指标变量作为主成分，代替原来的11个指标进行综合分析。5.1.3主成分分析结果根据以上步骤，利用matlab编程，首先求出各指标的相关性系数表，从表中可以发现，某些指标具有很强的相关性，如果直接用这些指标对影响力进行综合评估，必然造成信息的重叠，影响评价的客观性。主成分分析可以把多个指标转化成少数几个不相关的综合指标。相关系数矩阵的前几个特征值及其累计贡献率如表5-1所示5-1主成份分析结果序号特征值分项贡献率累计贡献率17.16590.65140.651421.33930.12180.773231.0310.09370.8669可以看出，前三个特征值的累计贡献率达到86.69%，效果很好，故本文中选取前三个主成份进行分析，其对应的特征向量如表5-2所示表5-2前三个特征值对应的特征向量第一第二第三-0.2520.3490.2920.5220.490.3290.1950.586-0.102-0.570.1810.370.3250.229-0.071-0.1990.377-0.527第一-0.210.236-0.002-0.0220.077第二-0.044-0.207-0.254-0.0520.009第三-0.1330.3650.4720.018-0.04由此，得出第一、二、三主成分分别为：（8）由主成分系数可以看出，第一主成分中，指标4占的比重最大；第二主成份中，指标2、4占的比重最大；第三主成份中，指标6、9占的比重最大。于是本文得出，影响水资源短缺风险的主要因子为：指标2、4、6、9。实际指标为蓄水量、降雨量、地下水埋深、人口总量。蓄水量与降雨量的的一再的减少，增加了水资源短缺的风险；人口总量的一再增加，使得用水量大量增加，在水资源来源减少的情况下，无疑大量增加了水资源供给的负担；而地下水位的一再下降，使得地下水的使用越来越难。这些都是使得潜在风险函数增加的因素。5.2问题二：水资源短缺风险等级的划分5.2.1等级的划分本文暂且将风险等级分为五类，即｛风险极限，高风险，中度风险，低风险，无风险｝，对应的数值为5，4，3，2，1.表5-3风险等级划分风险度风险极限高风险中度风险低风险无风险风险度函数（）54321注：风险极限是指达缺水程度到了城市的负荷极限根据实际情况，本文选取偏大型柯西分布隶属函数作为隶属函数其中、、、为待定系数当风险度为无风险时，取当风险度为中度风险时，隶属度为0.8，即当风险度达到极限时，隶属度为1，即计算出系数、、、分别为=1.1086、=0.8942、=0.3915、=0.3699于是得到风险度的隶属函数为（9）利用风险隶属函数计算出风险等级的隶属度如表5-4表5-4风险等级隶属度风险度风险极限高风险中度风险低风险无风险风险度函数（）10.91260.80.52450.015.2.2在第一问中，本文通过分析得到了主成份。将主成分分别代入下式，得到风险潜在函数（10）其中，为第个主成分的信息贡献率，风险潜在函数就是对水资源短缺风险的一个综合评估分别以主成分贡献量为权重，构建风险潜在函数为（11）将、、代入，得到风险潜在函数的表达式为：（12）5.2.3历年风险归类将预处理后的数据标准化处理后，带入风险潜在函数进行计算，得到历年的风险潜在值如表5-4所示：表5-4历年风险潜在值年份潜在风险年份潜在风险年份潜在风险1979-1.227921989-0.3889419991.2813251980-0.731941990-0.7361520000.9569551981-0.742271991-0.7758820011.0018131982-1.0490519920.86814520022.4767661983-0.8808219931.3310120030.6948391984-0.884431994-0.7673920040.5047531985-1.039381995-0.5553820050.8615781986-0.588681996-0.8746620061.3039731987-0.6071119970.03631420071.1701391988-0.926651998-0.5292820080.818082由表5-4所示，本文得出了历年的风险潜在值，但是仅仅这样无法将风险归类，于是本文将风险潜在值进行平移—极差化变换，得到历年风险隶属度如表5-5：表5-5历年风险隶属度年份风险隶属度年份风险隶属度年份风险隶属度1979019890.22646519990.6773161980013274120000.5897619810.1310919910.12201820010.60186819820.04828119920.5657872002119830.09369119930.69072820030.51900719840.09271819940.1243120040.46769719850.050892199505640151986009535520060.6834319870.16757419970.34125220070.64730419880.08132119980552274由表5-5可以看出，风险最小的一年是1979年，风险度为“无风险”，而风险最大的一年是2002年，该年度的风险隶属度达到了风险极限，并且概念的风险隶属度与每一年都相差较大。而本文所研究的水资源短缺风险是在不考虑突发状况的情况假设下的，数据应该是平缓变换的，于是本文采用均值替换法求该年的风险隶属度均值替换公式如下：替换后，2002年的风险隶属度为0.560438。使用matlab编程将历年划分等级，得到历年的潜在风险等级如表5-6：表5-6历年的风险等级年份风险等级年份风险等级年份风险等级197911989219993198021990220003198121991220013198221992320023198321993320032198421994220042198521995220053198621996220063198721997220073198821998220083由表5-6可以看出风险变化趋势为，前十年风险度都比较低，而中间十年开始逐渐出现中度风险的现象，到了后面十年风险度几乎都是中度风险再将风险隶属度变化做散点图得到图5-1，也可以得到相同的结论：图5-1历年风险隶属度变化趋势由此分析可以看出，如果不立即采取保护水资源的措施，水资源缺乏风险将会愈演愈烈5.2.4结果的合理性检验为了检验本文对水资源短缺划分的合理程度，本文采用聚类的方法将11个指标30年数据的变化进行聚类分析。暂不制定聚类数目，使用matlab编程进行聚类分析，得到的结果如图5-2：图5-211个指标30年数据的聚类结果从聚类图可以看出，30年的数据大致被分为3类，从1979—1994年的数据为一类，1995—1999年的数据为一类，2000—2008年的数据为一类。与本文5.2.3的分析结果相差不大，即随着城市的发展，水资源短缺的风险度由无风险、低风险逐渐在想中度风险变化。如不立即进行调控，后果将难以想象。5.2.5相关的调控措施5.2.5在问题一中，本文通过主成份分析的方法，找出对水资源短缺风险影响较大的因素为：1.蓄水量2.降雨量3.地下水位埋深4.人口总量。通过作图分析各个指标的变化趋势如图5-3：图5-3主要指标变化趋势由图5-3可以看出，北京市的蓄水量在逐年减少，降雨量变动较大，但是总体趋势也在递减。地下水位埋深不断增加为地下水的开采加大了难度。人口总量在逐渐增加，导致了生活用水的增加。有这些不利的趋势导致了水资源短缺的风险。1.蓄水量：蓄水量的大小由水利工程的数量与规模决定，在短期内无法做出有效的调整，目前北京只有官厅、密云水库。于是从长期来说，可以增加水库的容量，甚至是考虑增加水库的数量，以此增加蓄水量；其次是加强管理制度，科学合理的防水补源。2.地下水位：地下水位主要接受大气降水补给，其次为侧向径流补给、河渠渗漏补给、地表水灌溉和井灌回渗以及城市工业和生用水的排泄。随着城市规模扩大和人类活动的加剧，大量的工业废水、生活污水等通过河渠排放渗入地下，使大部分浅层地下水遭受污染。于是，采取的相应的措施为：适当的采取人工降雨，增加降雨量；提高废水处理率，减少上游污水的排放；3.降雨量：可以通过人工降雨的方法增加降雨量，增大绿化面积，维持一个正常的水循环。4.人口总量：人口总量对生活用水会产生直接的影响。由图5-3可以很直观的看到，总体上人口总量呈增长趋势，生活用水势必会大幅度增加，所以限制常住人口的数量也是一个重要的措施。5.3问题三北京市水资源短缺的预测为了对北京市未来两年水资源短缺的风险进行预测与研究，本文决定采用灰色模型GM(1，1)预测各项指标2009和2010年的数值，在将数据处理之后带入风险潜在函数计算风险潜在值，得到2009与2010年的水资源短缺情况，具体步骤如下：5.3.1灰色模型的建立分析处理后的数据，共有11个指标，30年的指标数据。由于过早的数据不具有代表性，于是本文取前后二十年的数据进行预测。将数据的数量级化为一致后，将11个指标的平均值记为（13）平均价格的级比规定范围为（14）经验证，级比均落在内。对数列进行一次累加，得到数列则记一次累加生成数列（15）取的加权均值，则，，记（16）于是可以得到GM（1,1）的白化微分方程为（17）取为灰导数，为背景值，则方程（12）所对应的灰微分方程为（18）则可得其中，用最小二乘法求得参数的估计值为（19）解得白化微分方程的特解为则可以解得（20）使用matlab编程求解，得到2009和2010的预测值如表5-7：表5-72009与2010各项指标预测值指标农业用水蓄水量入境水量降水量(毫米)污水处理能力200911.98148.52224.5662370.5357440.6512201011.58428.34954.2916358.9641494.4715指标地下水埋深工业用水第三产业人口总量生活水价工业水价200924.65236.11517.564117135.21143.2909201025.95425.826418.2521757.66.16293.73485.3.2预测结果分析将上述灰色模型预测值处理后，代入模型（12），即风险潜在函数表达式，得到2009年和2010年北京市的风险潜在值分别为：，；数据标准化后，其相应的风险等级均为3，风险程度为中度风险。由于预测的未知性，没有足够的依据证明本文预测的准确性，于是本文采用第二种预测方法，神经网络进行对比检验。.3.1神经网络的设计人工神经网络具有自组织、自适应和自学习能力,以及具有非线性、非局域性、非定性和非凸性等特点，并且在理论上可以逼近任何非线性函数。使用神经网络进行预测，不需要输入、输出之间明确的函数关系，主要通过对数据训练、学习完成模拟过程，并利用训练好的网络对新输入的数据进行预测。其中BP(BackPropagation)网络是反向传播的多层前馈式网络，是目前使用最为广泛的一种人工神经网络。国内学者利用多种预测方法进行预测研究，对比预测模型的预测精度，发现BP神经网络预测结果与实际值非常接近，效果很好，完全可以满足对预测的要求。故本文采取BP神经网络进行进一步的价格预测，然后将其预测结果与灰色模型的预测结果进行比较，从中选出预测精度较高的模型作为分析基础。对于BP神经网络,有一个非常重要的定理,即对于任何在闭区间内的一个连续函数以用单隐层的BP神经网络逼近。因而可以用一个三层BP神经网络就可以完成任意的m维到n维的映射,即一般的预测问题都可以通过单隐层的BP网络实现。所以,本文采用三层BP网络设计预测模型。首先，进行输入数据的选取。为了和灰色系统预测区分开来，本文决定采用神经网络直接对风险潜在值进行预测。由于过早的数据没有代表性，于是选取近20年的风险潜在值进行预测，网络输入层的神经元个数为20。输出层的神经元数为1。数据的预处理在进行网络预测之前，要对数据进行归一化处理，将数据处理为区间［0，1］之间的数据。于是本文采用极差化方法对数据进行归一处理。网络预测完成后，对预测结果进行分析时，应对预测数据进行反归一化，然后才能够与真实值进行比较。依本文的归一化式子，反归一化式子为（21）隐层的设计隐层神经元数目的确定是一个很复杂的问题，目前还没有科学的统一计算式。但有一些经验公式可用来大致确定隐层神经元的数目范围。本文采用下式：（22）其中：——输入层节点数——输出层节点数——常数。取值为1,2,3…10由上式知：其最佳隐层神经元数目可能在20到30之间。为了得到最佳的隐层神经元数，所以本文增大搜索范围,以10为初始值逐步增长，即先从一个简单的神经网络开始，若不符合要求，逐步增加隐层神经元数目，直到合适时为止。最终确定的最佳神经元数目以实验调试结果为准，上式确定的范围仅是一个供参考的初始范围。BP神经网络在matlab上的实现现在用matlab中的神经网络工具箱对2009年和2010年风险潜在值进行预测。先对BP神经网络结构进行初步设计。由于本文是一个用于预测，较为复杂的神经网络，输入量与输出量之间满足非线性关系，故输入层、输出层传递函数分别选取非线性函数S型正切函数tansig、对数函数logsig。再进行对训练函数的选择。由于BP神经网络有众多的训练函数，而训练函数的选取直接关系到预测所需步长和预测精度。为选取最合适的训练函数，本文取1989年至2007年的数据对2008年的数据进行预测。2008年的潜在风险值数据作为目标输出，1989年至2007年的数据作为输入，对神经网络进行训练。将各种训练函数得到的预测结果与2008年的实际风险值进行比较，得出的各组平均误差如表5-8：表5-8训练函数精度表训练函数相对误差trainrp36.95%trainlm2.22%traingd5.72%trainbfg27.25%traingdx24.61%由上表中结果容易看出，trainlm训练函数的预测结果最为精确。网络神经结构初步设计好后，用逐步增长法确定最佳隐层神经数。最后本文发现最佳隐层神经元数22。至此，BP神经网络的构建完毕。利用神经网络进行预测分析利用matlab编程实现预测，得到2009年度与2010年度的风险潜在值分别为，；该预测值与灰色系统的预测值相近，且风险等级都为3级，中度风险。由此可见，本文对于北京市2009年度与2010年度的水资源短缺风险预测具有一定的合理性与正确性。5.3.4北京市水资源短缺应对措施在本文的研究中，北京市历年水资源短缺的潜在风险值是由各个指标所决定的，所以要进行有效的调控，则必须注意11个指标的变化情况。1.北京市的人口总量有稳定增加的趋势，人口的增加必然会导致生活用水的增加，于是控制北京市常住人口的数量成为了一个主要因素。2.生活用水在逐步增加，而工业用水和农业用水在逐步减少，所以可以降低工业用水的水费，提高生活用水的水费，这样可以起到较好的限制作用而又可以不影响工业的发展。3.污水处理能力对可使用的地下水量有很大的影响，本文的预测中，污水处理能力在逐步增强，这也许是由于科技的增强所致。如果有关部门能尽力到达本文对污水处理能力的预测值，那么水资源的短缺得到有效的调控。5.4问题四：对北京市水行政主管部门的建议北京作为中国的首都，一个国际化的大都市，经济、文化等都相当发达。然而今天却面临着水资源短缺这一问题，水资源短缺成为了限制北京进一步发展的短板因素，如何调控相关影响因子，是北京的水资源短缺风险减轻至关重要。由题目中给定的数据和2009年北京统计年鉴上的数据可以看到，北京市每年的水资源总量在递减，且递减趋势较大，而用水总量却没有显著减少，如不进行调控，则水资源短缺问题将会越来越严重。根据本文对北京市水资源短缺问题的分析，找出了11个与水资源短缺联系的指标，这11个指标可以归结为3个方面：用水风险因子，来水风险因子，管理调节因子；11个指标分别为农业用水、工业用水、人口总量、第三产业及生活等用水；地下水位埋深、入境水量、降雨量；蓄水量、污水处理能力、工业水价、生活水价。而根据本文在第三文中做的预测，如果不及时调控水资源短缺的峰下因子，那么水资源短缺风险将进一步加剧。从本文的分析中，11个指标中的主要风险因子为地下水位埋深、人口总量、蓄水量、降雨量。这四个指标在影响水资源短缺风险中占的比重较大，所以应该首先考虑这4个因子的调控。本文经过分析后提出以下调控措施：1.节约用水应当从水费调控。北京市拥有将近1700万常住人口，生活用水量很大，适当提高水价是控制用水量的一个重要措施。而水价和用水量之间的关系有待进一步研究；2.北京市的蓄水量在总体趋势上有一定的减小。增加蓄水量是一个较为重要的措施。增加蓄水量可以通过跨流域调水等实现；3.地下水是一个相当重要的因素。随着每一年地下水位的逐渐下降，北京市水资源短缺的风险在逐步增加。提高污水处理能力是一个保护地下水的实用措施，这样可以让北京市拥有更多的可用地下水资源。4.在降雨量整体减小的情况下，人力对其的控制能力十分微弱，可以通过改善环境进行一些微弱的调控。除了以上针对主要风险因子采取的措施外，还可以通过采取其它方法对非主要因子进行调控，本文在此就不再赘述。六模型的推广与评价6.1模型的缺陷本文运用合理的方法按照题目的要求给出了模型与解答过程，但是也存在一些缺点。1.在2009和2010年的水资源短缺风险预测中，对人口的预测应该使用logstic模型，本文使用灰色系统进行预测，会造成一些误差，但是考虑到预测的年限较少，所以误差可以忽略。在其他的指标预测中，用于预测的年份数选取没有较好的方法。2.本文的研究中忽略了突发事件的影响，是模型过于理想化，这在以后的研究中应该进一步的完善。6.2模型的优点本文存在缺点的同时，也有一些优点。1.本文定义了风险潜在函数，对研究问题有很大的帮助。2.第一问采用主成分分析定义的权重无主观性，有较强的说服力。3.第二问中将风险等级划分后，本文并未草草了事，而是运用聚类分析进行验证；4.第三问中运用灰色模型对水资源缺乏的风险进行预测后，还运用精度较高的神经网络模型进行检验。这样的检验使得本文的分析较为严谨。同时，本文也针对分析结果对相关部门提了一些控制水资源短缺的实用建议。6.3模型的推广本文运用了科学的方法对水资源的短缺进行了评价和预测，有较强的推广性。本文的模型不仅可以用于水资源短缺的风险研究，还可以用于其它的很多项目，例如工程建设中的风险研究、工程造价中的成本研究、供应链中的影响因素的研究等。七参考文献[1]韩中庚，数学建模方法及其应用，北京：高等教育出版社，2005年6月第1版[2]姜启源、谢金星，数学建模案例精选，北京：高等教育出版社，2006年7月第1版[3]薛长虹、于凯，大学数学实验，成都：西南交通大学出版社，2003年10月第1版[4]北京市2009年度统计年鉴，/tjnj/2009-tjnj/，访问日期（2011-5）八附录8.1主成分分析源代码（matlab）clear,clcX=[7.2 0.0518 0.0499 0.0014 0.0431 24.18 14.37 4.37 897.1 8.3333 49.2 0.0518 0.0499 0.0026 0.0431 31.83 13.77 4.94 904.3 8.3333 49.2 0.0518 0.0499 0.0025 0.0397 31.6 12.21 4.3 919.2 8.3333 49.2 0.0518 0.0499 0.0018 0.0397 28.81 13.89 4.52 935 8.3333 49.2 0.0518 0.0499 0.002 0.0397 31.6 11.24 4.72 950 8.3333 49.2 0.0518 0.0499 0.002 0.0397 21.84 14.376 4.017 965 8.3333 49.2 0.0518 0.0499 0.0014 0.0397 10.12 17.2 4.39 981 8.3333 410 0.0518 0.0935 0.0015 0.0379 19.46 9.91 7.18 1028 8.3333 410 0.0518 0.0935 0.0015 0.0379 9.68 14.01 7.26 1047 8.3333 410 0.0518 0.0619 0.0015 0.0379 21.99 14.04 6.4 1061 8.3333 411 0.0518 0.158 0.0023 0.0379 24.42 13.77 6.45 1075 8.3333 410.7 0.0769 0.0853 0.0014 0.0329 21.74 12.34 7.04 1086 8.3333 410.7 0.0719 0.0779 0.0013 0.0329 22.7 11.9 7.43 1094 8.3333 4.347810.5 0.0939 0.0974 0.0018 0.2222 19.94 15.51 10.98 1102 3.5714 2.631612.8 0.1639 0.1522 0.002 0.2222 20.35 15.28 9.59 1112 4.1667 3.030311.7 0.0526 0.0718 0.0012 0.0408 20.93 14.57 10.37 1125 5 3.571411.6 0.0767 0.0579 0.0017 0.0171 19.33 13.78 11.77 1251.1 5.5556 410.2 0.0511 0.0398 0.0014 0.0171 18.95 11.76 9.3 1259.4 3.5714 2.777812.1 0.1153 0.1176 0.0023 0.0171 18.12 11.1 11.1 1240 3.7037 2.564112 0.0684 0.0692 0.0014 0.0171 17.39 10.84 12.2 1245.6 2.6316 1.470614.2 0.2564 0.1727 0.0037 0.0171 18.45 10.56 12.7 1257.2 1.8182 1.08715.2 0.2907 0.1406 0.0027 0.0078 16.49 10.52 13.39 1363.6 1.1111 0.854716.2 0.216 0.189 0.003 0.007 17.4 9.2 12.3 1385.1 0.8929 0.833317.4 0.578 0.3846 0.0027 0.0055 15.5 7.5 11.6 1423.2 0.7299 0.740718.4 0.1071 0.2392 0.0022 0.0047 13.8 8.4 13.6 1456.4 0.6061 0.662319.04 0.093 0.1582 0.0021 0.0039 13.5 7.7 13.4 1492.7 0.5405 0.602420.21 0.0837 0.2179 0.0024 0.0031 13.2 6.8 14.5 1538 0.4785 0.549521.52 0.0814 0.2353 0.0031 0.003 12.8 6.2 15.3 1581 0.4292 0.507622.79 0.0904 0.2899 0.0021 0.0028 12.4 5.8 16.6 1633 0.3906 0.469522.92 0.0773 0.1869 0.0016 0.003 12 5.2 17.9 1695 0.3571 0.4386];Z=zscore(X);[pc,score,variance,t2]=princomp(Z)percent=variance/sum(variance)leiji=cumsum(percent)8.2聚类分析源代码（matlab）clc,cleara=[7.2 0.0518 0.0499 0.0014 0.0431 24.18 14.37 4.37 897.1 8.3333 49.2 0.0518 0.0499 0.0026 0.0431 31.83 13.77 4.94 904.3 8.3333 49.2 0.0518 0.0499 0.0025 0.0397 31.6 12.21 4.3 919.2 8.3333 49.2 0.0518 0.0499 0.0018 0.0397 28.81 13.89 4.52 935 8.3333 49.2 0.0518 0.0499 0.002 0.0397 31.6 11.24 4.72 950 8.3333 49.2 0.0518 0.0499 0.002 0.0397 21.84 14.376 4.017 965 8.3333 49.2 0.0518 0.0499 0.0014 0.0397 10.12 17.2 4.39 981 8.3333 410 0.0518 0.0935 0.0015 0.0379 19.46 9.91 7.18 1028 8.3333 410 0.0518 0.0935 0.0015 0.0379 9.68 14.01 7.26 1047 8.3333 410 0.0518 0.0619 0.0015 0.0379 21.99 14.04 6.4 1061 8.3333 411 0.0518 0.158 0.0023 0.0379 24.42 13.77 6.45 1075 8.3333 410.7 0.0769 0.0853 0.0014 0.0329 21.74 12.34 7.04 1086 8.3333 410.7 0.0719 0.0779 0.0013 0.0329 22.7 11.9 7.43 1094 8.3333 4.347810.5 0.0939 0.0974 0.0018 0.2222 19.94 15.51 10.98 1102 3.5714 2.631612.8 0.1639 0.1522 0.002 0.2222 20.35 15.28 9.59 1112 4.1667 3.030311.7 0.0526 0.0718 0.0012 0.0408 20.93 14.57 10.37 1125 5 3.571411.6 0.0767 0.0579 0.0017 0.0171 19.33 13.78 11.77 1251.1 5.5556 410.2 0.0511 0.0398 0.0014 0.0171 18.95 11.76 9.3 1259.4 3.5714 2.777812.1 0.1153 0.1176 0.0023 0.0171 18.12 11.1 11.1 1240 3.7037 2.564112 0.0684 0.0692 0.0014 0.0171 17.39 10.84 12.2 1245.6 2.6316 1.470614.2 0.2564 0.1727 0.0037 0.0171 18.45 10.56 12.7 1257.2 1.8182 1.08715.2 0.2907 0.1406 0.0027 0.0078 16.49 10.52 13.39 1363.6 1.1111 0.854716.2 0.216 0.189 0.003 0.007 17.4 9.2 12.3 1385.1 0.8929 0.833317.4 0.578 0.3846 0.0027 0.0055 15.5 7.5 11.6 1423.2 0.7299 0.740718.4 0.1071 0.2392 0.0022 0.0047 13.8 8.4 13.6 1456.4 0.6061 0.662319.04 0.093 0.1582 0.0021 0.0039 13.5 7.7 13.4 1492.7 0.5405 0.602420.21 0.0837 0.2179 0.0024 0.0031 13.2 6.8 14.5 1538 0.4785 0.549521.52 0.0814 0.2353 0.0031 0.003 12.8 6.2 15.3 1581 0.4292 0.507622.79 0.0904 0.2899 0.0021 0.0028 12.4 5.8 16.6 1633 0.3906 0.469522.92 0.0773 0.1869 0.0016 0.003 12 5.2 17.9 1695 0.3571 0.4386];[m,n]=size(a);d=zeros(m,m);fori=1:mforj=i+1:md(i,j)=mandist(a(i,:),a(j,:)');endenddnd=nonzeros(d);nd=union(nd,nd)fori=1:m-1nd_min=min(nd);[row,col]=find(d==nd_min);tm=union(row,col);tm=reshape(tm,1,length(tm));s(i)={char(['第',int2str(i),'次合成，平台高度为',num2str(nd_min),'时的分类结果为：',int2str(tm)])};%上面大括号{}代表建立数组nd(find(nd==nd_min))=[];iflength(nd)==0breakendends(:)y=pdist(a,'cityblock');yc=squareform(y)z=linkage(y)[h,t]=dendrogram(z)8.3灰色系统预测源代码（matlab）clc,clear

x0=[1110.710.710.512.811.711.610.212.11214.215.216.217.418.4

19.0420.2121.5222.7922.92];

n=length(x0);xx=1:20;

plot(xx,x0,'o-')

holdon

disp('级比')

lamda=x0(1:n-1)./x0(2:n)

range=minmax(lamda)

x1=cumsum(x0)

fori=2:n

z(i)=0.4*x1(i)+0.6*x1(i-1);

end

B=[-z(2:n)',ones(n-1,1)];

Y=x0(2:n)';

u=B\Y

x=dsolve('Dx+a*x=b','x(0)=x0');

x=subs(x,{'a','b','x0'},{u(1),u(2),x1(1)});

yuce1=subs(x,'t',[0:n-1]);

digits(6),y=vpa(x)%为提高预测精度，先计算预测值，再显示微分方程的解

yuce=[x0(1),diff(yuce1)]

disp('残差')

epsilon=x0-yuce%计算残差

disp('相对误差')

delta=abs(epsilon./x0)%计算相对误差

disp('级比偏差值')

rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda%计算级比偏差值

%若以上各项检验符合标准，则以下可预测将后的所需的数据

yuce2=subs(x,'t',[0:n-1+2]);%这里的n-1为验证模型时的数据量的上限，可以加上或减去一定的数值进行相应的预测

digits(2),y=vpa(x)%为提高预测精度，先计算预测值，再显示微分方程的解

yuce3=[x0(1),diff(yuce2)]

yy=1:22;

plot(yy,yuce3,'+-')

holdoff;8.4神经网络预测源代码（matlab）%通用BP神经网络P=[0.1449293520.0413286760.0294740210.5200189410.6581288080.0320072710.095266910.010.2718168330.103054640.6433037760.5465181090.5599028590.5141053390.4683078190.4115898720.5180594770.6500614960.6101280470.505081121]';t=[0.6794];H=[0.0413286760.0294740210.5200189410.6581288080.0320072710.095266910.010.2718168330.103054640.6433037760.5465181090.5599028590.5141053390.4683078190.4115898720.5180594770.6500614960.6101280470.5050811210.6794]';net=newff([01;01;01;01;01;01;01;01;01;01;01;01;01;01;01;01;01;01;01;01],[22,1],{'tansig','logsig'},'trainlm');%输入参数依次为:'样本P范围',[各层神经元数目],{各层传递函数},'训练函数'%训练函数traingd--梯度下降法，有7个训练参数.%训练函数traingdm--有动量的梯度下降法,附加1个训练参数mc(动量因子，缺省为0.9)%训练函数traingda--有自适应lr的梯度下降法,附加3个训练参数:lr_inc(学习率增长比，缺省为1.05;%lr_dec(学习率下降比，缺省为0.7);max_perf_inc(表现函数增加最大比，缺省为1.04)%训练函数traingdx--有动量的梯度下降法中赋以自适应lr的方法，附加traingdm和traingda的4个附加参数%训练函数trainrp--弹性梯度下降法,可以消除输入数值很大或很小时的误差，附加4个训练参数:%delt_inc(权值变化增加量，缺省为1.2);delt_dec(权值变化减小量，缺省为0.5);%delta0(初始权值变化，缺省为0.07);deltamax(权值变化最大值,缺省为50.0)%适合大型网络%训练函数traincgf--Fletcher-Reeves共轭梯度法;训练函数traincgp--Polak-Ribiere共轭梯度法;%训练函数traincgb--Powell-Beale共轭梯度法%共轭梯度法占用存储空间小,附加1训练参数searchFcn(一维线性搜索方法,缺省为srchcha);缺少1个训练参数lr%训练函数trainscg--量化共轭梯度法,与其他共轭梯度法相比，节约时间.适合大型网络%附加2个训练参数