數(shù)學建模之長江水質(zhì)監(jiān)測問題

1、長江水質(zhì)監(jiān)測摘要本文解決的是長江水質(zhì)的評價與監(jiān)測問題，通過分析過去十年不同監(jiān)測站收集到的長江水質(zhì)數(shù)據(jù)，運用不同的理論建立不同的模型，對長江過去十年的水質(zhì)情況作出評價，然后再預(yù)測未來十年長江水質(zhì)的變化情況。針對問題一：考慮到問題一中需要對長江水質(zhì)情況作出定量的評價，并分析各地區(qū)水質(zhì)的污染狀況，為此，建立模糊綜合評價模型確定了其隸屬度函數(shù)，建立評判因子的權(quán)重矩陣，求得最終結(jié)果為：水質(zhì)最差的地方是江西南昌滁槎 （15號），其次水質(zhì)差的地方為四川樂山岷江大橋（8號）、湖南長沙新港（12號）以及四川瀘州沱江二橋（10號），此四處水質(zhì)污染嚴重；水質(zhì)最好的地方是湖北丹江口胡家?guī)X （11號

2、）。針對問題二：根據(jù)長江的降解系數(shù)，可得到污染物隨時間的變化量。由于污染源的污染物排放量等于本地區(qū)污染物的流量與上游流下的污染物流量之差。因此，建立污染物流量隨時間變化的微分方程模型。最后求得：高錳酸鉀指數(shù)和氨氮的污染源主要集中在宜昌至岳陽之間。針對問題三：根據(jù)已知的過去10年的主要統(tǒng)計數(shù)據(jù)，建立了灰色預(yù)測模型。 在相對誤差較小的情況下對未來10年的水質(zhì)情況作出了預(yù)測，分析得出結(jié)論：未來10年可飲用水所占的比例越來越低，排污量有明顯的上升趨勢。針對問題四：在問題四中建立多元線性回歸方程，利用最小二乘法求解系數(shù)，在滿足問題四要求的前提下，求出未來10年的允許最大相對排污量，繼而求得未來

3、10年每年的相應(yīng)排污量，后者與前者的差值與未來10年的長江水總流量的乘積，求得最終結(jié)果如下表： 未來10年預(yù)處理的排污量年代2005200620072008200920102011201220132014預(yù)處理排污量(億噸)71.2483.1194.98106.86118.73130.60142.48154.35166.22178.09針對問題五：分析總結(jié)前幾個問題的結(jié)果，找出水質(zhì)污染的根本原因。結(jié)合考察團的調(diào)查結(jié)果，給出合理的建議和意見。 最后，對模型中運用的方法進行了優(yōu)、缺點評價，在模型的推廣中提出了可以建立類似模型解決生活中的一類問題。關(guān)鍵詞：模糊評價 微分方程 灰色預(yù)測 線性

4、回歸 1問題重述1.1問題背景長江是我國第一、世界第三大河流，其水質(zhì)的污染程度日趨嚴重，已引起了相關(guān)政府部門和專家們的高度重視。2004年10月，由全國政協(xié)與中國發(fā)展研究院聯(lián)合組成“保護長江萬里行”考察團，從長江上游宜賓到下游上海，對沿線21個重點城市做了實地考察，揭示了一幅長江污染的真實畫面，其污染程度讓人觸目驚心。為此，專家們提出“若不及時拯救，長江生態(tài)10年內(nèi)將瀕臨崩潰”。 長江水是許多人賴以生存和發(fā)展的資源，保護水資源就是保護我們自己，就長江近年來的水質(zhì)情況，采取合理的保護和治理措施刻不容緩。1.2解決問題為了制定出合理的治理長江水質(zhì)污染的方案，根據(jù)長江地區(qū)近兩年多主要水質(zhì)指

5、標的檢測數(shù)據(jù)，現(xiàn)討論以下幾個問題： (1) 對長江近兩年多的水質(zhì)情況做出定量的綜合評價，并分析各地區(qū)水質(zhì)的污染狀況。  (2) 研究、分析長江干流近一年多主要污染物高錳酸鹽指數(shù)和氨氮的污染源的主要分布地區(qū)。  (3) 假如不采取更有效的治理措施，依照過去10年的主要統(tǒng)計數(shù)據(jù)，對長江未來水質(zhì)污染的發(fā)展趨勢做出預(yù)測分析。  (4) 根據(jù)預(yù)測分析，如果未來10年內(nèi)每年都要求長江干流的類和類水的比例控制在20%以內(nèi)，且沒有劣類水,那么每年需要處理多少污水？  (5) 發(fā)表對解決長江水質(zhì)污染問題切實可行的

6、建議和意見。2問題分析問題需要對長江水質(zhì)作出評價和預(yù)測，根據(jù)附件中已知的數(shù)據(jù)，建立相應(yīng)的評價預(yù)測模型，分析得出長江水質(zhì)過去十年以及未來十年的污染狀況。2.1問題一的分析問題提供的數(shù)據(jù)主要反映了三方面的內(nèi)容即：水質(zhì)劃分等級標準、17個測站點在28個月份中主要污染物的濃度和水質(zhì)等級。而觀測點水質(zhì)等級的確定是由主要污染物的濃度決定的，因此為了對長江及各觀測點作出定量的綜合評價我們需要將主要污染物的濃度歸一到一個比較量中。這個比較量要能夠快速而準確的反映出各測站點的水質(zhì)污染狀況，因此我們評價整個長江流域的水質(zhì)狀況采用模糊綜合評價模型。 雖然PH值是影響水質(zhì)等級的因素之一，但是給出的所有PH值

7、都在69這個范圍內(nèi)，即PH值對于水質(zhì)等級幾乎沒有影響，所以在以后的問題討論中不再考慮PH值的影響。2.2問題二的分析因為污染物在時間和空間上是動態(tài)變化的，為簡便起見，不考慮支流的因素，只從縱向的角度，結(jié)合本地與上游污水進行水質(zhì)的污染分析。根據(jù)常識可知，污染物的污染源就是新增污染物較多的地方。因此，考慮到降解系數(shù)，求得污染物隨時間的變化量，通過相鄰兩主要觀測點之間每千米的相對排污量判定污染源，即本地新增污染物的量等于本地區(qū)現(xiàn)有的污染物減去上游通過降解后到達該處的污染物的總量，再比上兩地之間的距離，得到相對排污量。建立污染物的量隨時間變化的微分方程模型，最后通過新增污染物的相對量大小關(guān)系得出主要的

8、污染源地區(qū)。2.3問題三的分析測第三問根據(jù)過去10年長江的總體水污染狀況的監(jiān)測數(shù)據(jù)，可以看出長江總體水流量變化不大，但年排污的總量在增加，這使得污染河段比例增加，污染的嚴重程度呈現(xiàn)快速增長的趨勢，即每年污染情況主要與當年的排污量和總水流量等因素有關(guān)。為此，首先利用回歸分析方法確定出可飲用水的比例與總排污量和總水流量的關(guān)系式，然后根據(jù)過去10 年排污量，利用灰色預(yù)測方法對未來的年排污量做出預(yù)測，最后根據(jù)總排污量的增長趨勢來推斷出可飲用水比例的變化趨勢，從而可以預(yù)測出未來10 年長江水質(zhì)的變化情況。2.4問題四的分析 在問題四中建立多元線性回歸方程，利用最小二乘法求解系數(shù)，在滿

9、足未來10年內(nèi)每年都要求長江干流的類和類水的比例控制在20%以內(nèi)，且沒有劣類水的前提下，求出未來10年的允許最大相對排污量，繼而求得未來10年每年的相應(yīng)排污量，后者與前者的差值與未來10年的長江水總流量的乘積，即可求得未來10年要處理的排污量。2.5問題五的分析分析總結(jié)前幾個問題的結(jié)果，找出水質(zhì)污染的原因。結(jié)合考察團的調(diào)查結(jié)果，給出合理的建議和意見。3模型假設(shè)與符號說明3.1模型假設(shè)假設(shè)一： 降解系數(shù)在一定時間段固定不變；假設(shè)二：長江干流的自然凈化能了近似均勻； 假設(shè)三：干流的相鄰兩觀測站點之間的排污口主要集中在下游處；假設(shè)四：干流污染物的富集主要受上游影響，支流的影響忽略

10、不計。3.2符號說明符號符號意義氧氣屬于第i等級的隸屬度高錳酸鹽屬于第i等級的隸屬度氨氮屬于第i等級的隸屬度第i種指標的權(quán)重各觀測點模糊關(guān)系矩陣第i種指標的實時均值 平均排污量 相對排污量 斷面水流平均流速 某組分子在處的濃度 降解常數(shù) 每年長江總水流量 每年總排污量 V類水所占的比例 劣V類水所占比例4模型的建立與求解4.1問題一的解答4.1.1模型一的建立與求解（1）確立評判指標：根據(jù)地表水環(huán)境質(zhì)量標準（GB38382002）中4個主要項目標準限值設(shè)置表，結(jié)合28個月的各觀測點的測量結(jié)果，選取溶氧量（DO）、高錳酸鹽指數(shù)(CODMn)和氨氮（NH3-N）為評價因子，設(shè)置評價因子集：，求出1

11、7個觀測點28個月的個評價因子的均值，得到17個點的評價因子集。（2）建立評價集 序 號項目標 準值分類類 類 類 類 類 劣類 1 溶解氧(DO)    7.5（或飽和率90%）65 3 2 0 2 高錳酸鹽（CODMn) 246 10 15 3氨氮（NH3-N） 0.15 0.51.0 1.5 2.0 4 PH值（無量綱） 6-94個指標將水質(zhì)分為6個等級。 （3）建立隸屬函數(shù)，進行單因素評價：由于水質(zhì)污染程度和水質(zhì)分級標準都是模糊的，所以用隸屬度來描述分級界限較為合理，先根據(jù)各指標的6級標準，做出6個級別的隸屬函數(shù)，其中DO的評價指標以數(shù)值大為優(yōu)，其余兩個指

12、標以數(shù)值小為優(yōu)，函數(shù)如下：DO（A）以升半梯形分布建立隸屬函數(shù) COOMn(B)和NH3-N（C）以升半梯形分布建立隸屬函數(shù)：有了各指標的隸屬函數(shù)，就可以進行單因素評價，將各觀測點的3個指標帶入相應(yīng)的隸屬函數(shù)，計算出指標的隸屬度，得到各觀測點模糊關(guān)系矩陣（j=1，2，3.，17），具體數(shù)據(jù)見附錄一。 （4）建立評價因素的權(quán)重集：由于DO,COOMn和NH3-N等污染指標對水質(zhì)的影響不同，因此對各指標應(yīng)賦予不同的權(quán)重，根據(jù)污染物對水質(zhì)的污染大、權(quán)重的的原則來決定權(quán)重的大小。定義：對于DO越大越優(yōu)型：，對于COOMn和NH3-N越小越優(yōu)型：，其中，為第i種指標的權(quán)重，是第i種指標的實時均值，分別為

13、多級濃度標準的最大值和最小值。再將權(quán)重歸一化處理，評價因素權(quán)重的集合為A=0.145,0.226,0.629。（5） 模糊綜合評價： 評價因素的權(quán)重集A乘以單因子矩陣得到模糊綜合評價結(jié)果，。將6個等級看作一種相對位置使其連續(xù)化，設(shè)各等級分別用1，2，3，4，5，6表示，并稱為各等級的秩，然后用B中對應(yīng)的量將各級的秩加權(quán)求積，得到被評等級的相對位置，即： 式中為隸屬于第k級的隸屬度。具體計算結(jié)果見下表：表1：各觀測站水質(zhì)等級序號123456789水質(zhì)級1.29911.41161.40381.45691.24041.50001.16482.50001.4450序號1011121314151617水

14、質(zhì)級2.12541.00002.38881.91561.50004.15721.43241.41654.1.2問題一的結(jié)果分析由表1可知，各觀測站中水質(zhì)等級最大的為江西南昌滁槎 （15號），可知此處水質(zhì)情況最差，其次水質(zhì)差的地方為四川樂山岷江大橋（8號）、湖南長沙新港（12號）以及四川瀘州沱江二橋（10號），此四處水質(zhì)污染嚴重；水質(zhì)最好的地方是湖北丹江口胡家?guī)X （11號）。4.2問題二的解答 因為污染物在時間和空間上是動態(tài)變化的，為簡便起見，不考慮支流的因素，只從縱向的角度，結(jié)合本地與上游污水進行水質(zhì)的污染分析。建立污染物流量隨時間變化的微分方程模型，最后通過新增污染物的流

15、量大小關(guān)系得出主要的污染源地區(qū)。4.2.1模型二的建立與求解（1）某一污染物擴散所滿足的微分方程是一個拋物線方程，結(jié)合實際問題的假設(shè)，常假設(shè)其水流近似的處于穩(wěn)定狀態(tài)，斷面均勻，普通對流擴散方程為 上式是解決污染物的一般方程，結(jié)合問題二的情況，假設(shè)觀測站干流的污染物濃度在一個月內(nèi)保持不變，則，又因為,所以上式進一步簡化為： 其中，為斷面水流平均流速，為某組分子在處的濃度，為降解常數(shù)取0.2。由此解出長江干流污染物濃度與據(jù)觀測站米處的函數(shù)關(guān)系： （2）對于長江上任一段干流AB，A為第個觀測點，B為第個觀測點，兩觀測點的距離為，假設(shè)該段干流之間有個排污口（包括支流入口和直排口），為計算的方便，假設(shè)排

16、污口聚集在一起，其總的排污口的流量為，平均流量為，污染物的濃度為，即為AB段的排污量,則 對于任意一個江段AB，污染物的濃度，水流量，流速，距離和降解系數(shù)k都已知，則可以算出AB段污染物的排放量。（3） 根據(jù)所給數(shù)據(jù)可以算出，對于每個月每個江段的排污區(qū)間（j=1,2,.,13個月），取13個月的平均值表示平均排污量，則該江段每千米的排污量為相對排污量，每個江段都有一個相對排污量，以此為指標可以判斷長江干流的主要污染源。4.2.2模型二的求解按步驟運用MATLAB編求解（程序見附錄3），結(jié)果如下表： 表2：COOMn的排污量 江段 32510.5435266.77 4600229855.6211

17、104.1531404.31相對排污量 34.22 45.33 116.46 59.71 105.75 60.05排序651423 表3：NH3-N的排污量 江段 2822.62 3191.46 4451.77 2559.31 1191.38 158.46相對排污量 2.97 4.10 11.27 5.12 7.26 0.34排序 5413264.2.3問題二的結(jié)果分析問題二用長江干流7個觀測站點將長江分為6個江段，通過污染物流量隨時間的變化關(guān)系建立了微分方程模型，先計算出每月每段的高錳酸鹽和氨氮的量，再對六個污染源近一年所排放的污染物的數(shù)量求期望，分析得出污染物高錳酸鹽指數(shù)和氨氮的污染源主要

18、集中在宜昌至岳陽之間。 4.3問題三的解答4.3.1模型三的建立與求解當系統(tǒng)內(nèi)部信息和特性是部分已知的，另一部分是未知時，人們往往難以建立客觀的物理原形，內(nèi)部因素難以辨識，以及相互之間的關(guān)系較為隱蔽，難已準確了解這類系統(tǒng)的行為特征，因此，對于這類問題進行定量描述，即建立模型難度較大，便選擇建立灰色系統(tǒng)，對問題求解。1. 灰色動態(tài)模型 附件4中給出了1995-2004年長江各等級在豐水期，枯水期，水文年的百分比，為了預(yù)測此趨勢之下未來十年長江水質(zhì)的變化，通過對水文年全流域水質(zhì)的變化研究來對長江未來水質(zhì)污染的發(fā)展趨勢作出預(yù)測分析。建立灰色預(yù)測模型，首先舍去異常數(shù)據(jù)：1995年I

19、V類、劣V類數(shù)據(jù)，然后從附件4中提取出水文年全流域隨時間變化的原始數(shù)據(jù)。  對于原始數(shù)據(jù)進行累加處理，即   加生成數(shù)列為 ()()()(= 采用一階單變量微分方程進行擬合，得到白化方程的GM（1,1）模型 + =  其中式中的為待定系數(shù)。 綜上所述，灰色動態(tài)模型為：  =+-í+=î 2. 模型求解 取水文年全流域的飲用水預(yù)測為例，可通過對其變化研究對長江未來水質(zhì)污染的發(fā)展趨勢作出預(yù)測分析。 (1) 求解步驟： 第一步：從附

20、表中提取出數(shù)據(jù)  第二步：利用公式對做累加生成得   第三步：構(gòu)造矩陣B和數(shù)據(jù)向量，與滿足關(guān)系其中  第四步：求解常系數(shù)é ùéù=êúêúëûëû 第五步：得出表達式：  ´´ 令 +=-  所以表示還原后的值。 (2) 結(jié)果求解 按步驟運用MATLAB編求解（程序見附錄4），運用此程序根據(jù)過去十年的數(shù)據(jù)可以預(yù)測出未來10年的排污量和可

21、飲用水如下表：年代2005200620072008200920102011201220132014排污量303.01322.52343.29365.39388.92413.96440.61468.98499.18531.32 表4：未來10年排污量年份2005200620072008200920102011201220132014可飲水69.6168.0266.4764.9563.4762.0260.658.2256.8754.55 表5：未來10年可飲用水所占比例4.3.2問題三的結(jié)果分析從表5可以看出，可飲用水所占的比例越來越低，所以，可以看出水質(zhì)越來越差，從2012年起，可飲用水的百分比

22、甚至低于60%,如此嚴重的水污染，應(yīng)該引起足夠重視，采用合理的措施來減少對長江水的污染。4.4問題四的解答4.4.1模型的建立與求解 （1）長江干流的類和類水以及劣V類水的比例與每年長江總水流量和總排污量有關(guān)，因此分別以類和類水所占的比例和劣V類水所占比例為響應(yīng)變量,以歷年來的 = =為解釋變量，做出多元線性回歸，得到類和類水的總和及劣V類水與長江總水流量及總排污量的關(guān)系，做出一般的多元線性回歸模型：運用matlab用最小二乘法求解得到回歸系數(shù) 就求得了類和類水的總和及劣V類水與長江總水流量及總排污量的關(guān)系。問題四要解決問題是：未來10年內(nèi)每年都要求長江干流的類和類水的比例控制在20%以內(nèi)，且

23、沒有劣類水。則約束條件為： 綜上所述，所建立的模型為： 對上述回歸方程進行回歸檢驗時，發(fā)現(xiàn)回歸效果并不是太好，考慮到長江干流的類和類水的比例與江總水流量及總排污量不是直接的線性關(guān)系，因此對模型進行改進。（2） 優(yōu)化模型 長江干流的類和類水以及劣V類水的比例與每年總排污量和長江總水流量的比值相對排污量有關(guān)，因此分別以類和類水所占的比例和劣V類水所占比例為解釋變量，為響應(yīng)變量，做多元線性回歸。對表四的數(shù)據(jù)進行觀測，得到與時間成線性關(guān)系，因此列出如下線性回歸方程： 運用matlab用最小二乘法進行求解，得到回歸系數(shù)分別為 計算時剔除了2003年這個突變點，且檢驗到回歸性較好。得到的回歸方程如下： 將

24、帶入的表達式中，得到最大允許相對排污量，將帶入的表達式中得到未來10年每年預(yù)測的相對排污量，則污水處理量，因為每年的總水流量（除特殊年份如1998年）變化不大，故取前10年的總水流量均值為每年的總水流量，得到每年的污水處理量如下表：表6： 未來10年預(yù)處理的排污量年代2005200620072008200920102011201220132014預(yù)處理排污量(億噸)71.2483.1194.98106.86118.73130.60142.48154.35166.22178.094.4.2問題四的結(jié)果分析從表6可以看出，未來10年所需處理的污水越來越多，所以，可以看出水質(zhì)越來越差，如此嚴重的水污

25、染，應(yīng)該引起足夠重視，采用合理的措施來減少對長江水的污染。4.5解決長江水質(zhì)問題的建議通過問題一的計算可以看出，長江水質(zhì)的枯水期、豐水期與平水期都由水流量的不同導致水質(zhì)情況也不同。從問題三、四的求解結(jié)果可以看出，在未來的十年里，長江水質(zhì)不斷惡化，廢水排放量逐年遞增，需要處理的污水也越來越多。根據(jù)模型假設(shè)，在不發(fā)生特大洪澇災(zāi)害時長江的水流量是一定的，總結(jié)得出解決長江水質(zhì)問題的關(guān)鍵，在于減少廢水的排放量的同時增強廢水的處理能力。鑒于此，我們認為可以從以下幾個方面著手解決長江水質(zhì)問題：  （1）加大處罰力度。 （2） 對污染嚴重地區(qū)進行重點治理。 （3） 加強宣傳，增強

26、民眾保護長江的意識。 （4） 嚴格進行檢測，處理違規(guī)企業(yè)，嚴懲腐敗。 （5） 保護長江立法，引起政府和人民的高度重視。5模型的評價、改進與推廣5.1模型的評價5.1.1模型的優(yōu)點（1） 問題一中建立了模糊綜合評價模型，通過精確的數(shù)字手段處理模糊的評價對象，對蘊藏信息呈現(xiàn)模糊性的資料作出了比較科學、合理、貼近實際的量化評價；（2） 問題二的模型，簡化了許多因素，忽略了支流的影響，但這種理想化還是合理的，較符合實際情況;（3） 問題三采用灰色預(yù)測模型，具有少數(shù)據(jù)性、良好的時效性、較強的系統(tǒng)性和關(guān)聯(lián)性等特征，可以合理的對數(shù)據(jù)做出預(yù)測。5.1.2模型的缺點 在求解模型時，為了使問題

27、得到方便的解決，往往采用簡化的手段進行求解，因此求出的結(jié)果與真實值會存在一定的偏差。5.2模型的改進在問題三中運用了灰色動態(tài)模型對未來10年長江的水質(zhì)作出了預(yù)測，但并未對預(yù)測的結(jié)果進行檢驗，不知預(yù)測結(jié)果的準確度有多高，在此，最好能將模型三的基礎(chǔ)上對灰色模型的預(yù)測結(jié)果進行檢驗，以保證預(yù)測結(jié)果有較高的可信度。5.3模型的推廣模型一中，除了用模糊算法評價外，還可以用BP神經(jīng)網(wǎng)絡(luò)來求長江水質(zhì)的在綜合評價而；而在求解問題三和問題四時采用的灰色預(yù)測，還可以用回歸分析預(yù)測、時間序列預(yù)測和神經(jīng)網(wǎng)絡(luò)預(yù)測等代替。6參考文獻1 宋曉秋編著，模糊數(shù)學原理與方法，北京，中國礦業(yè)大學出版社，1999 

28、2 張俊福等編著，應(yīng)用模糊數(shù)學，北京，地質(zhì)出版社，20013 岳超源，決策理論與方法，北京：科學出版社，20034 鄧聚龍，灰理論基礎(chǔ),武漢：華中科技大學出版社，20025 韓中庚，數(shù)學建模方法及其應(yīng)用，北京：高等教育出版社，2005 附錄一各觀測點模糊關(guān)系矩陣：1號 2號3號 4號5號 6號7號 8號9號10號11號 12號13號 14號15號16號17號附錄二問題一的程序（以17號觀測站為例）：clc;cleara=1 0 0 0 0 0;0.5107 0.5107 0&

29、#160;0 0 0;0.3917 0.3917 0 0 0 0;b=0.145 0.226 0.629;c=b*a;d=0;e=0;for i=1:6    d=d+c(1,i)*i;endfor i=1:6    e=e+c(1,i);endB=d/e附錄三問題二的程序：clc;clearnum11=xlsread('11.xls');num29=xlsread('29.xls&

30、#39;);for i=1:6    Ca=num11(i,6);    Cb=num11(i+1,6);    Qa=num29(2,i+2);    Qb=num29(2,i+3);    d1=num29(1,i+1);    d2=num29(1,i+2);    d=(d2-d1)*1000;    Ua=num29(3,i+2);    Ub=num29(3,i+3);    U=(86400*(Ua+Ub)/2;    m=exp(-0.2*d/U);    W=Cb*Qb-Ca*Qa*m;   &#