數(shù)學建模論文水資源短缺風險綜合評價

1、論文題目：水資源短缺風險綜合評價摘要水資源，是指可供人類直接利用，能夠不斷更新的天然水體。主要包括陸地上的地表水和地下水。風險，是指某一特定危險情況發(fā)生的可能性和后果的組合。本文提出了馬氏判別法、模糊聚類、bp神經(jīng)網(wǎng)絡等三種方法對北京市水資源短缺風險進行綜合評價，針對問題一基于附表1通過馬氏判別法篩選出影響水資源短缺的主要風險因子，針對問題二通過模糊聚類的方法，分了水資源短缺的四個等級，在問題三中通過構建神經(jīng)網(wǎng)絡，測出了2010年、2011年和2012年的水資源總量和用水總量，為解決水資源短缺風險，提出了南水北調、再生水的利用、污水處理等幾種措施，并分析了在進行這幾項措施后歷年風險等級的下降情

2、況，最后向水行政主管部門書寫了一份建議報告，基于建立的水資源短缺風險評價模型提出了建議。關鍵詞：馬氏判別法、模糊聚類、bp神經(jīng)網(wǎng)絡一、問題的重述 水資源短缺風險，泛指在特定的時空環(huán)境條件下，由于來水和用水兩方面存在不確定性，使區(qū)域水資源系統(tǒng)發(fā)生供水短缺的可能性以及由此產(chǎn)生的損失。北京是世界上水資源嚴重缺乏的大都市之一，其人均水資源占有量不足300m3，為全國人均的1/8，世界人均的1/30，屬重度缺水地區(qū)，附表中所列的數(shù)據(jù)給出了1979年至2000年北京市水資源短缺的狀況。北京市水資源短缺已經(jīng)成為影響和制約首都社會和經(jīng)濟發(fā)展的主要因素。如何對水資源風險的主要因子進行識別，對風險造成的危害等級進

3、行劃分，對不同風險因子采取相應的有效措施規(guī)避風險或減少其造成的危害，這對社會經(jīng)濟的穩(wěn)定、可持續(xù)發(fā)展戰(zhàn)略的實施具有重要的意義?？紤]以下問題：1 評價判定北京市水資源短缺風險的主要風險因子；2 建立一個數(shù)學模型對北京市水資源短缺風險進行綜合評價， 作出風險等級劃分并陳述理由。對主要風險因子,如何進行調控，使得風險降低？3 對北京市未來兩年水資源的短缺風險進行預測，并提出應對措施。4 以北京市水行政主管部門為報告對象，寫一份建議報告。二、模型假設及符號說明1、模型的假設（1）以一年為一個周期，考慮在一個周期內的水資源短缺情況；（2）相鄰的任意兩年內水資源來源和使用均互不影響；（3）不考慮影響水資源短

4、缺情況較小的因素，如：降雨量、地下水位埋深；（4）只考慮影響水資源短缺情況較大的因素（工業(yè)用水、農(nóng)業(yè)用水、第三產(chǎn)業(yè)及生活等其他用水和水資源總量等）；（5）影響水資源短缺較大的因素之間相互獨立；（6）未來幾年內，國家政策及相關因素對水資源影響不大；2、符號說明 缺水量 風險程度 模糊集 最少的缺小量 最大的缺水量 三、問題的分析問題一，要求評價判定北京市水資源短缺風險的主要風險因子。主要風險因子也就是從對水資源影響較大的因素中篩選出提供信息量最大且之間相關程度較低的因子，馬氏距離是一種有效的計算兩個未知樣本集的相似度的方法，我們采用馬氏距離判別法使得相鄰兩類之間的馬氏距離最大，從而篩選出主要風險

5、因子。問題二，這是一個模糊聚類問題。根據(jù)北京市1979年到2008年的水資源狀況計算出每年的缺水量，通過構造隸屬函數(shù)得出各年份的風險程度，再通過spss軟件對30年風險程度進行系統(tǒng)聚類，再根據(jù)結果進行模糊聚類，進而劃分出風險等級。再通過對比措施前后水資源短缺風險等級得出措施的有效性。問題三，通過bp神經(jīng)網(wǎng)絡預測模型預測數(shù)未來兩年的水資源短缺的風險程度，從而確定風險等級，再根據(jù)相應主要因子的調控提出相應解決方案。 我們綜合考慮問題一到題三，在此基礎上書寫建議報告，完成問題四。四、模型的建立及求解1、對于問題一： 影響北京水資源短缺風險的因素可歸納為以下兩個方面：（1）自然因素：人口數(shù)；入境水量；

6、水資源總量；地下水位埋深；（2）社會經(jīng)濟環(huán)境因素：工業(yè)用水量；污水處理率；農(nóng)業(yè)用水量；第三產(chǎn)業(yè)及生活等其他用水量。對附表一及北京市2009年統(tǒng)計年鑒中的數(shù)據(jù)研究分析得出：影響水資源短缺的主要因素為工業(yè)用水量、農(nóng)業(yè)用水量、第三產(chǎn)業(yè)及生活等其他用水量和水資源總量。通過逐步判別法對四各主要因素進行篩選，在判斷過程的每步種引入判別能力較強的變量或者剔除判別能力較弱的變量，并進行統(tǒng)計檢驗，最終的判別函數(shù)中的變量即為主要風險因子。 根據(jù)缺水量=總用水量-水資源總量及附表一數(shù)據(jù)信息計算出每年的缺水量，用0和1（1表示缺水量為正值，即該年度缺水，0表示缺水量為負值，即該年度不缺水）對缺水量進行分組，將樣本集合

7、按照缺水量分組情況劃分成k=2個總體?，F(xiàn)共有n=30的樣本，每個樣本有p=4項觀測指標，利用=3個變量來對k個總體進行判別。（1）計算出樣本點的總協(xié)方差矩陣（2）計算樣本點的組內離差和而后進行變量選擇：定義為矩陣t和w的行列式比值，即=，逐一檢驗工業(yè)用水量、農(nóng)業(yè)用水量、第三產(chǎn)業(yè)及生活等其他用水量和水資源總量對判別效果的貢獻，設判別函數(shù)中已經(jīng)有了農(nóng)業(yè)用水量，記為x，這時考慮增加變量工業(yè)用水量x，計算偏統(tǒng)計量： =(x,x)/(x)有結論： ，如果，則可增加變量工業(yè)用水量，反之，刪除工業(yè)用水量，按此循環(huán)直至四個主要因素都進行全部通過判別函數(shù)。使用spss軟件，帶入具體30年各指標用水量的數(shù)值，分組

8、定義范圍為1,2，方法選擇mahalanobisdistance mahalanobis距離方法，的概率為0.05，的概率為0.10，得出如下結論：分析中的變量步驟容差要刪除的 f 的顯著水平最小 d 方組之間1水資源總量1.000.0042水資源總量.879.001.5841.00 和 2.00農(nóng)業(yè)用水.879.0322.7721.00 和 2.003水資源總量.847.0285.3911.00 和 2.00農(nóng)業(yè)用水.371.0012.7831.00 和 2.00第三產(chǎn)業(yè)及生活等其他用水.359.0154.8301.00 和 2.00表11 逐步判別法進入變量不在分析中的變量步驟容差最小容差

9、要輸入的 f 的顯著性最小 d 方組之間0農(nóng)業(yè)用水1.0001.000.166.5841.00 和 2.00工業(yè)用水1.0001.000.132.6961.00 和 2.00第三產(chǎn)業(yè)及生活等其他用水1.0001.000.183.5381.00 和 2.00水資源總量1.0001.000.0042.7721.00 和 2.001農(nóng)業(yè)用水.879.879.0324.8301.00 和 2.00工業(yè)用水.891.891.6382.8641.00 和 2.00第三產(chǎn)業(yè)及生活等其他用水.852.852.8752.7831.00 和 2.002工業(yè)用水.678.668.1146.1561.00 和 2.0

10、0第三產(chǎn)業(yè)及生活等其他用水.359.359.0158.1971.00 和 2.003工業(yè)用水.438.233.9268.2031.00 和 2.00表12 逐步判別法剔除變量由表11可知，水資源總量農(nóng)業(yè)用水及第三產(chǎn)業(yè)及生活等其他用水刪除f的顯著性較小，在表12中，該三項不在分析中的顯著性較小，而工業(yè)用水不在分析中的顯著性較大，故進入判別函數(shù)的變量為農(nóng)業(yè)用水、第三產(chǎn)業(yè)及生活等其他用水、水資源總量三項因素，即影響水資源短缺的主要風險因子為農(nóng)業(yè)用水、第三產(chǎn)業(yè)及生活等其他用水和水資源總量。2、對于問題二：2.1、風險等級的劃分建立一個數(shù)學模型對北京市水資源短缺風險進行綜合評價，針對風險的模糊性，定義模

11、糊集合，其中表示風險程度，表示缺水量，構造隸屬函數(shù)，將1979-2008年的缺水量換算成對應的風險程度用0到1區(qū)間的數(shù)字表示；構造的隸屬函數(shù)如下： ，其中表示缺水量的最小值，表示缺水量的最大值，p為大于等于1的正整數(shù)，不妨取p為2，分別計算北京30年的風險程度如下表1：年份缺水量（x億立方米）風險程度（u（x）19794.690.1241198024.540.8394198124.110.8172198210.620.2712198312.860.341519840.740.05761985-6.290.001619869.520.23961987-7.71019883.250.0969198

12、923.090.765619905.260.13581991-0.260.0448199223.990.811199325.550.892819940.450.0537199514.540.39961996-5.860.0028199718.070.536419982.730.088199927.491200023.540.7882200119.70.6064200218.50.5544200317.40.5089200413.20.3529200511.30.291720069.80.24742007110.282520080.90.0598對于1979-2008年的缺水量與風險程度的比較：

13、通過以上兩圖的比較可以說明隸屬度真正反映了30年缺水的狀況，很好的反應了風險的模糊性。通過對北京30年間風險程度的不同進行分類采用系統(tǒng)聚類的方法，開始時每個對象自成一類，然后每次將最相似的兩類合并，合并后重新計算新類與其他類的距離或相似度，這一過程一直繼續(xù)，一直到所有對象歸為一類為止，并類的過程用一張普類系譜圖表示。由spass軟件可求得對北京30年風險的分類樹狀圖：由該圖可知，對于這30年的數(shù)據(jù)共可分為30組、5組、4組、2組、1組。根據(jù)系統(tǒng)聚類的結果，然后采用模糊聚類的方法分別分出5組、4組、2組，鑒于30組和1組情況太過極端對風險的分類太過分散和統(tǒng)一，因此舍去這兩種情況，分別采用模糊聚類

14、的方法對風險分成5組、4組和2組。模糊聚類的步驟：第一步：建立數(shù)據(jù)觀測矩陣sample；sample=4.69 24.54 24.11 10.62 12.86 0.74 -6.29 9.52 -7.71 3.25 23.09 5.26 -0.26 23.99 25.55 0.45 14.54 -5.86 18.07 2.73 27.49 23.54 19.7 18.5 17.4 13.2 11.3 9.8 11 0.9;第二步：標定模糊相似矩陣data=(),矩陣元素表示向量和的相似程度，由隸屬函數(shù)w可求得data；data= 0.1241 0.8394 0.8172 0.2712 0.341

15、5 0.0576 0.0016 0.2396 0 0.0969 0.7656 0.1358 0.0448 0.8110 0.8928 0.0537 0.3996 0.0028 0.5364 0.0880 1.0000 0.7882 0.6064 0.5544 0.5089 0.3529 0.2917 0.2474 0.2825 0.0598第三步：在matlab中編程實現(xiàn)聚類分析，找出中心點center；（1）、分成兩組的情況：風險等級類中心風險描述級低度風險0.1529可忽略的風險級高度風險0.7571無法承受的風險對于三十年的劃分圖：可以看出對30年的等級劃分并不明顯，不能起到警示和對水資

16、源短缺風險程度的正確描述，因此這種劃分標準方法并不可行。（2）分成5組的情況風險等級類中心風險描述風險等級0.0538可忽略的風險風險等級0.2621低度風險風險等級0.3624中度風險風險等級0.5531較高風險風險等級0.8399高度風險圖中可知，五級的劃分在較低風險時過于密集，在風險程度為0.5以下時，在實際情況下，造成的危害并不是很大，因并沒有必要對0.5以下的風險劃分三個級別。(3)、分成4組的情況風險等級類中心風險描述風險等級0.8410高風險風險等級0.5473中度風險風險等級0.2933低風險風險等級0.0571可忽略的風險有圖中可知，劃分四級對于三十年的數(shù)據(jù)具有較好的劃分，且

17、落于每個等級的年份與實際情況也相符?；谝陨系姆治?，風險等級應分為4級即風險等級類中心風險描述級風險等級0.8410高風險級風險等級0.5473中度風險級風險等級0.2933低風險級風險等級0.0571可忽略的風險通過計算觀測量與各中心的聚類，根據(jù)距離最小的原則，把各觀測量分派到各類中心所在的類中去。2.2、對主要風險因子的調控由問題1可知水資源總量、農(nóng)業(yè)用水、第三產(chǎn)業(yè)及生活等其他用水為主要風險因子，因此調控的對象主要是針對這三個因素。在第三產(chǎn)業(yè)及生活等其他用水中，其中生活用水是不可壓縮的，隨著北京都市化進程的不斷加快，人口增長與人民生活水平的不斷提高，生活用水量會進一步加大，因此，只能通過減

18、少北京的人口數(shù)量來間接的減少生活用水量，同時也會相應的減少第三產(chǎn)業(yè)的用水量，但北京作為首都，這樣的舉措會實施起來相對困難。北京農(nóng)業(yè)用水量占北京總用水量的40左右，但由于受到基本農(nóng)田保護制度的政策的制約，進一步大幅度壓縮農(nóng)業(yè)用水的可能性不大。在問題1中可以看出水資源總量對水資源短缺有比較大的影響，因此增加水資源總量，將會大大緩解北京的缺水狀況，通過比較從2008年開始的南水北調工程可以看出2008、2009年的風險情況。如下圖;由圖中明顯可以感覺出2009年的風險等級下降了一個等級，對北京缺水狀況的改善相當?shù)拿黠@，在2008年中由于奧運會等特殊原因且加上南水北調首次實施，因此風險下降還不十分明顯

19、。 在用水方面可以提高再生水的利用，比較2005年到2009年再生水使用情況，如下圖：由圖中可以看出再生水的使用對緩解北京水資源的短缺效果顯著，特別是在2009年時，極大的緩解了水資源短缺情況。 科學技術的逐漸發(fā)展，污水處理技術水平也越來越來高，污水處理效率逐漸升高，因此提高對處理過的污水的利用將極大的改善北京現(xiàn)在的用水狀況，若被處理的污水能夠得到充分的利用。模擬的2001年到2009年的風險的等級圖 如下：上圖可知，對于污水的利用將極大改善北京的缺水狀況，在途中我們可以看出：年份風險等級處理前處理后2001高度風險中度風險2002高度風險中度風險2003高度風險中度風險2004高度風險低度風

20、險2005高度風險低度風險2006高度風險低度風險2007高度風險低度風險2008低度風險可忽略的風險2009高度風險低度風險3、對于問題三：未來水資源短缺風險預測模型bp神經(jīng)網(wǎng)絡的數(shù)學模型的建立，首先是基于下面最基本的數(shù)學公式。沿信息傳播的方向，給出網(wǎng)絡的狀態(tài)方程，用和 表示第層第 個神經(jīng)元的輸入和輸出，網(wǎng)絡的各層輸入輸出關系可描述為：第一層（輸入層）：將輸入引入神經(jīng)網(wǎng)絡。 第二層（隱層）： 其中是激發(fā)函數(shù)，可以取不同形式，如：s 函數(shù)：高斯基函數(shù)：以及徑向基函數(shù)、樣條基函數(shù)、小波函數(shù)等，本文采用的是高斯基函數(shù)。第三層（輸出層）： 基于本文模型的條件：i=10，k=1,l=9,已上就完成了神

21、經(jīng)網(wǎng)絡模型的基本網(wǎng)絡構建。下一步是要確定網(wǎng)絡的學習算法。習的基本思想是：通過一定的算法調整網(wǎng)絡的權值，使網(wǎng)絡的實際輸出盡可能接近期望的輸出。在本網(wǎng)絡中，采用誤差反傳（bp）算法來調整權值（即 bp 神經(jīng)網(wǎng)絡）。其基本原理是：當網(wǎng)絡的輸入（即相應影響因素）為 時，網(wǎng)絡的實際輸出為，網(wǎng)絡的期望輸出（即實際的學費值）為 。定義學習的目標函數(shù)為：（均方差法）同時 bp 算法通過下列公式來調整權值，使目標函數(shù)最?。海ㄗ钏傧陆捣ǎ┢渲袨閷W習率。其中具體分析每一層神經(jīng)網(wǎng)絡，可以得到：其中基于以下公式，偏差逐步反傳：.得出各因素的處理后權值： 根據(jù)得到的權值便可進行對未來數(shù)據(jù)的預測。這里用bp人工神經(jīng)網(wǎng)絡模型

22、對未來年份的水資源總量及總用水量作出預測，進而得到未來年份的缺水量，通過前面得到的隸屬函數(shù)求得缺水程度，最后根據(jù)模糊聚類得到的風險等級劃分標準確定未來年份的水資源風險等級。具體實現(xiàn)過程如下：首先，將已知年份的水資源總量、總用水量整理成將要用matlab進行bp網(wǎng)絡訓練的輸入數(shù)據(jù)與目標數(shù)據(jù)。其次，將已得數(shù)據(jù)進行歸一化處理，建立bp神經(jīng)網(wǎng)絡，設定必要參數(shù)并進行訓練，得到合適的網(wǎng)絡。建立bp神經(jīng)網(wǎng)絡模型，分別對2010、2011、2012年的水資源總量和用水總量進行預測； 建立的神經(jīng)網(wǎng)絡示意圖通過用建立的神經(jīng)網(wǎng)絡擬合以往的歷史數(shù)據(jù)如下圖所示：通過上圖發(fā)現(xiàn)建立的bp神經(jīng)網(wǎng)絡與以往的歷史數(shù)據(jù)擬合較好，可

23、以用來預測；通過matlab的求解，得出2010、2011、2012年的數(shù)據(jù)如下：年份水資源總量（億立方米）用水總量（億立方米）缺水總量（億立方米）201025.716135.1689.4519201113.32434.12120.797201221.641634.496112.8545通過問題二建立的風險綜合評價，對2010、2011、2012年的數(shù)據(jù)進行風險等級的劃分如下圖所示：通過上圖可以看出2010年、2012年屬于低度風險（級）、2011年屬于中度風險（級），基于以上結論，我們采取以下措施：在2010年、2012年時缺水的危險等級較低，特別是在2010年時風險等級與級極為接近，201

24、2年時缺水等級雖然屬于低度危險，但風險等級與中度風險等級也極為接近，因此在2012年時應該做好防范工作。在這三年中2011年的危險等級最高，基本接近級，因此在2011年需要采取相應的解決方案應對缺水危機。在2011年和2012年應加大對于北京的入境量，提高南水北調的送水量，另一方面應該應提高再生水的利用率和污水處理能力，使得污水能盡快得到循環(huán)利用。4、對于問題四：通過對歷年北京市缺水風險的觀察可以看出：隨著年份的增加，風險不趨于穩(wěn)定處于較大的波動狀態(tài)，同時，預測未來兩年的缺水等級也不宜樂觀。為了使北京市的順利發(fā)展和社會的和諧進步，本著“節(jié)流優(yōu)先，治污為本，多渠道開源”的城市水資源可持續(xù)利用的新

25、戰(zhàn)略，提出如下建議：（1）節(jié)約用水節(jié)約用水是當務之急也是長遠發(fā)展的戰(zhàn)略方針，在優(yōu)先保證城市生活用水和重點工業(yè)供水的前提下，要適度的壓縮農(nóng)業(yè)用水、，加強工業(yè)、農(nóng)業(yè)節(jié)水力度，調整產(chǎn)業(yè)結構，大力發(fā)節(jié)水型工業(yè)、農(nóng)業(yè)。實施滴灌代替大水澆灌，讓每一滴水都盡量達到百分百的利用；收集并有效利用雨洪，既可防治雨洪災害，緩解城市雨洪壓力，同時又增加了可用水資源等。最后，還要提高人們的節(jié)水意識，從節(jié)約生活中的每一滴水做起降低，盡量降低人們的生活用水量。最好價格杠桿來控制人們的用水量，這樣也可以有效的節(jié)約水資源的利用（2）提供用水為保證城市供水安全，應科學適當?shù)脑黾拥叵滤_采量，合理開發(fā)利用，對已確定的應急供水水源地

26、應盡快投入開發(fā)和勘探工作，對其他地區(qū)繼續(xù)開展調查工作，尋找新的后備應急水源。南水北調進京后，提高了北京城市供水安全保證程度，支持了城市的可持續(xù)發(fā)展，所以在條件允許的前提下，適當提高南水北調的水量也可增加供水量。重要的一點，城市要大力的著手于污水的加工處理，從模型中可以得出，如果加大污水處理量或是提高污水處理率，可以大大緩解供水壓力，從而降低水資源短缺的風險。五、結果分析與檢驗該模型只考慮到影響水資源短缺的幾個主要因素，并未考慮其它次要因素，并且影響水資源的因素較多，也并不能找全所有的影響因素，降雨量的影響對于水資源的影響還是比較大的，一般來說，降雨量充沛的年份，缺水的風險就會較低，而干旱少雨的

27、年份缺水的風險就較大，基于這個原理，為檢驗模型，考慮1979到2009年各年的降雨量大小與各年風險的等級進行比較。如下圖所示：從圖中明顯可以看出降雨量大的年份，水資源短缺的風險的等級明顯較低，而在降雨量較小的年份，水資源短缺的風險等級明顯較高，因此可確定，在只考慮這幾個主要風險因子的情況下，我們這個模型仍然可以較為準確的評價水資源短缺的風險。在問題三中采用了bp神經(jīng)網(wǎng)絡對水資源總量和總水量進行了估計，通過訓練網(wǎng)絡，可以達到用十年的歷史數(shù)據(jù)估計下一年的用水總量和水資源總量，從而估計出該年的缺水量，為檢驗模型的正確性對以往2005、2006、2007的歷史數(shù)據(jù)進行估計，結果：年份總用水量水資源總量

28、實際預測實際預測200534.535.132723.224.8321200634.333.381524.523.875200734.835.264323.824.6547可以看出預測的數(shù)值與實際的情況誤差極小，因此該網(wǎng)絡對未來兩年的數(shù)據(jù)的預測的可信度較高，得到的結果較有說服力。六、模型的優(yōu)缺點1.模型的優(yōu)點： （1）問題一的模型忽略了次要因素的干擾，通過馬氏距離判別法，篩選出影響水資源短缺風險的主要因子，馬氏距離可以檢驗出兩個樣本之間的相似度，馬氏距離越大，兩個樣本的相關程度越小，因此采用這種方法對主要因子的篩選還是比較準確的； （2）問題二根據(jù)風險的模糊性，引入模糊集，建立隸屬函數(shù)，以每年的

29、缺水量作為自變量，定義出風險的程度，從而正確表示1979-2008年北京30年的缺水情況，并通過模糊聚類的方法對風險等級進行了劃分； （3）問題三建立bp神經(jīng)網(wǎng)絡，通過對水資源總量和用水總量數(shù)據(jù)處理訓練出網(wǎng)絡，然后對歷史數(shù)據(jù)估計，誤差較小情況下，對2010年、2011年、2012年進行預測，并對2011年可能會出現(xiàn)中度危險盡心了預警；2.模型的缺點： 在問題一中并未考慮其它次要因素，在對北京市水資源短缺的主要風險因子的評價過程中并未考慮其他因素的影響，這可能會造成對主要風險因子的劃分有偏差；問題二中隸屬函數(shù)的確定，并沒有系統(tǒng)的科學規(guī)律，在對風險程度的計算時，未必符合實際情況；問題三中神經(jīng)網(wǎng)絡的

30、預測中該模型還不是很穩(wěn)定，在對以往的歷史數(shù)據(jù)的預測時還會有較大的偏差。七、參考文獻1焦淑華等.2009.bp神經(jīng)網(wǎng)絡預測的matlab實現(xiàn).哈爾濱金融高等專科學校學學報2左其亭等.2003.水資源系統(tǒng)中的不確定性及風險分析方法.干旱區(qū)地理3北京市發(fā)改委課題組.2004解決北京市水資源問題的思路和主要措施.宏觀經(jīng)濟研究4王紅瑞.2009.基于模糊概率的水資源短缺風險評價模型及其應用.水利學報5范金城等.2002.數(shù)據(jù)分析.科學出版社6北京市水務局.2005北京市水資源公報7天津市人民政府.2009.關于下達2009年度南水北調工程基金地下水征收水量指標和基金上繳額度的通知.8張士鋒等.2010.

31、年水資源系統(tǒng)風險構成及其評價以北京市為例.自然資源學報9衷平等.2005.石羊河流域水資源短缺風險敏感因子的確定.干旱區(qū)資源與環(huán)境10李玉榕等.2006.一種基于馬氏距離的線性判別分析分類算法.計算機仿真11劉舒野等.2007.基于歸并聚類中心思想的模糊聚類分析方法研究.航空計算技術12朱起煌等.2000.概率模擬方法有助于預測項目的時限與成本.維普資訊13陳慶祥.1990.非等可能性隨即時間的概率模擬.電腦應用八、附錄：附表一：1979年2008年水資源情況年份總用水量(億立方米)農(nóng)業(yè)用水(億立方米)工業(yè)用水(億立方米)第三產(chǎn)業(yè)及生活等其它用水(億立方米)水資源總量（億方）197942.92

32、24.1814.374.3738.23198050.5431.8313.774.9426198148.1131.612.214.324198247.2228.8113.894.5236.6198347.5631.611.244.7234.7198440.0521.8414.3764.01739.31198531.7110.1217.24.3938198636.5519.469.917.1827.03198730.959.6814.017.2638.66198842.4321.9914.046.439.18198944.6424.4213.776.4521.55199041.1221.7412.

33、347.0435.86199142.0322.711.97.4342.29199246.4319.9415.5110.9822.44199345.2220.3515.289.5919.67199445.8720.9314.5710.3745.42199544.8819.3313.7811.7730.34199640.0118.9511.769.345.87199740.3218.1211.111.122.25199840.4317.3910.8412.237.7199941.7118.4510.5612.714.22200040.416.4910.5213.3916.86200138.917.

34、49.212.319.2200234.615.57.511.616.1200335.813.88.413.618.4200434.613.57.713.421.4200534.513.26.814.523.2200634.312.86.215.324.5200734.812.45.816.623.8200835.1125.217.934.2附表二：1979年2008年年降雨量統(tǒng)計年份降雨量（mm）年份降雨量（mm）1979718.41994813.21980380.71995572.51981393.21996700.91982544.41997430.91983489.91998731.71

35、984488.81999266.919857212000371.11986665.32001338.91987683.92002370.41988673.32003444.91989442.22004483.51990697.32005410.71991747.920063181992541.52007483.91993506.72008626.3附表三：2005年2009年污水總量、可再生水量、南水北調量年份污水總量可再生水量南水北調200511.32.6200612.93.6200712.984.9520085.9670.70220096.52.6matlab代碼：example.mcle

36、arclcdata=zeros(30,1);year=1979:2008;sample=4.69 24.54 24.11 10.62 12.86 0.74 -6.29 9.52 -7.71 3.25 23.09 5.26 -0.26 23.99 25.55 0.45 14.54 -5.86 18.07 2.73 27.49 23.54 19.7 18.5 17.4 13.2 11.3 9.8 11 0.9;sample=reshape(sample,30,1);for i=1:1:30 for j=1:1:1 data(i,j)=subordimproved(sample(i,j); ende

37、ndcenter,u,obj_fcn=fcm(data,4);plot(year,data,-ko,linewidth,2, markeredgecolor,k, markerfacecolor,r,markersize,10);axis(1979 2008 0 1);legend();title(1979-2008),xlabel(),ylabel();figure;plot(year,sample,o-);hold on;plot(year,sample,linestyle,none,marker,*,color,r);axis(1979 2008 -10 30);legend();tit

38、le(1979-2008),xlabel(),ylabel();center danger.mclearclcdata=zeros(30,1);year=1979:2008;sample=4.69 24.54 24.11 10.62 12.86 0.74 -6.29 9.52 -7.71 3.25 23.09 5.26 -0.26 23.99 25.55 0.45 14.54 -5.86 18.07 2.73 27.49 23.54 19.7 18.5 17.4 13.2 11.3 9.8 11 0.9;sample=reshape(sample,30,1);for i=1:1:30 for

39、j=1:1:1 data(i,j)=subordimproved(sample(i,j); endendplot(year,data,o-);axis(1979 2008 0 1);title(1979-2008),xlabel(),ylabel();y1=0.1753.*ones(30,1);y2=0.4206.*ones(30,1);y3=0.6943.*ones(30,1);hold onplot(year,y1,-,year,y2,-,year,y3,-);legend(, , , );danger1.mclearclcdata=zeros(30,1);year=1979:2008;s

40、ample=4.69 24.54 24.11 10.62 12.86 0.74 -6.29 9.52 -7.71 3.25 23.09 5.26 -0.26 23.99 25.55 0.45 14.54 -5.86 18.07 2.73 27.49 23.54 19.7 18.5 17.4 13.2 11.3 9.8 11 0.9;sample=reshape(sample,30,1);for i=1:1:30 for j=1:1:1 data(i,j)=subordimproved(sample(i,j); endendplot(year,data,o-);axis(1979 2008 0

41、1);title(1979-2008),xlabel(),ylabel();y1=0.455.*ones(30,1);hold onplot(year,y1,r-);legend(, );danger2.mclearclcdata=zeros(30,1);year=1979:2008;sample=4.69 24.54 24.11 10.62 12.86 0.74 -6.29 9.52 -7.71 3.25 23.09 5.26 -0.26 23.99 25.55 0.45 14.54 -5.86 18.07 2.73 27.49 23.54 19.7 18.5 17.4 13.2 11.3

42、9.8 11 0.9;sample=reshape(sample,30,1);for i=1:1:30 for j=1:1:1 data(i,j)=subordimproved(sample(i,j); endendplot(year,data,o-);axis(1979 2008 0 1);title(1979-2008),xlabel(),ylabel();y1=0.15795.*ones(30,1);y2=0.31225.*ones(30,1);y3=0.45775.*ones(30,1);y4=0.6965.*ones(30,1)hold onplot(year,y1,r-,year,

43、y2,r-,year,y3,r-,year,y4,r-);legend(, , , , );north.mclearclcsample1=0.9 13.7;data1=zeros(1,2);for i=1:1:2 data1(i)=subordimproved(sample1(i);endyear1=2008:2009;year2=2007:2009;subplot(1,2,1)plot(year1,data1,o,linewidth,2, markeredgecolor,k, markerfacecolor,r,markersize,10);y1=0.1753.*ones(3,1);y2=0

44、.4206.*ones(3,1);y3=0.6943.*ones(3,1);hold onplot(year2,y1,-,year2,y2,-,year2,y3,-);title();legend(, , , );subplot(1,2,2)sample2=1.602 16.3;data2=zeros(1,2);for i=1:1:2 data2(i)=subordimproved(sample2(i);endplot(year1,data2,o,linewidth,2, markeredgecolor,k, markerfacecolor,r,markersize,10);hold onpl

45、ot(year2,y1,-,year2,y2,-,year2,y3,-);title();predict.mclearclcsample1=9.45 20.797 12.855;data1=zeros(1,3);for i=1:1:3 data1(i)=subordimproved(sample1(i);endyear1=2010:2012;year2=2009:2013;plot(year1,data1,o-,linewidth,2, markeredgecolor,k, markerfacecolor,r,markersize,10);hold ony1=0.1753.*ones(5,1)

46、;y2=0.4206.*ones(5,1);y3=0.6943.*ones(5,1);hold onplot(year2,y1,-,year2,y2,-,year2,y3,-);title(2010-2012);legend(, , , );wasterecyle.mclearclcsample1=19.7 18.5 17.4 13.2 11.3 9.8 11 0.9 13.7;data1=zeros(1,9);for i=1:1:9 data1(i)=subordimproved(sample1(i);endyear1=2001:2009;year2=2000:2010;subplot(1,

47、2,1)plot(year1,data1,-o,linewidth,2, markeredgecolor,k, markerfacecolor,r,markersize,10);y1=0.1753.*ones(11,1);y2=0.4206.*ones(11,1);y3=0.6943.*ones(11,1);hold onplot(year2,y1,-,year2,y2,-,year2,y3,-);title();axis(2000 2010 0 1);legend(, , , );subplot(1,2,2)sample2=24.9 25.1 25.2 22.5 23.1 21.9 23.9

48、 12.96 26.7;data2=zeros(1,9);for i=1:1:9 data2(i)=subordimproved(sample2(i);endplot(year1,data2,-o,linewidth,2, markeredgecolor,k, markerfacecolor,r,markersize,10);hold onplot(year2,y1,-,year2,y2,-,year2,y3,-);title();axis(2000 2010 0 1);subordimproved.mfunction u=subordimproved(x)if(x-7.71&x27.49) u=(x+7.71)/35.2).2; else u=1; endendrecyle.mclearclcsample1=11.3 9.8 11 0.9 13.7;data1=zeros(1,5);for i=1:1:5 data1(i)=subordi