2019年基于貴州省部分區(qū)域雷達監(jiān)測降雨量的分析.doc

1、第一屆“九章數學杯”數學建模挑戰(zhàn)賽承諾書我們仔細閱讀了第一屆“九章數學杯”數學建模挑戰(zhàn)賽的競賽規(guī)則。我們完全明白，在競賽開始后參賽隊員不能以任何方式（包括電話、電子郵件、網上咨詢等）與隊外的任何人（包括指導教師）研究、討論與賽題有關的問題。我們知道，抄襲別人的成果是違反競賽規(guī)則的, 如果引用別人的成果或其他公開的資料（包括網上查到的資料），必須按照規(guī)定的參考文獻的表述方式在正文引用處和參考文獻中明確列出。我們鄭重承諾，嚴格遵守競賽規(guī)則，以保證競賽的公正、公平性。如有違反競賽規(guī)則的行為，我們接受相應處理結果。我們允許海南大學數學建模協(xié)會公布論文，以供同學之間學習交流，海南大學數學建模協(xié)會網站以非

2、商業(yè)目的的論文交流不需要提前取得我們的同意。參賽隊員（簽名 ） ：隊員1：趙滎仁隊員2：李銳3：秦昊基于貴州省部分區(qū)域雷達監(jiān)測降雨量的分析摘要針對雷達檢測系統(tǒng)對貴州省2015 年 5 月 7 日 0 時至 11 時降雨天氣的記錄，本文將采用插值擬合、積分求和等方式求解總降雨量。模型修正改進分析中利用查閱資料所得結果總結出最佳有效模型。對于問題一，我們首先使用相關繪圖軟件對雷達回波圖中的貴陽市進行剪輯分割， 對貴陽市區(qū)域面積上各時間段的雷達回波圖的降水量進行樣點采集和預處理。 其次， 在不考慮雷達系統(tǒng)出現故障以及每半小時內出現降雨量不均等等問題下使用MATLABI預處理之后的雷達回波圖轉化為灰度

3、圖，將之轉化成矩陣以知道每個像素點的數值，再通過數值統(tǒng)計計算出每種降雨量程度的分布面積。最終根據面積與基本反射率乘積在時間上進行累加，從而計算出貴陽市逐小時降雨量及總降雨量。對于問題二，根據附錄表格及雷達監(jiān)測范圍125km以內確定包括貴陽市在內 的清鎮(zhèn)市、修文縣、平壩縣、花溪區(qū)、長順縣、開陽縣、惠水縣、龍里縣、貴定縣、都勻市、麻江縣、息烽縣等13 個區(qū)域進行降雨量總量的強度的比較。在問題一模型的基礎上運用插值擬合分析，并在整個降雨時間內積分求和以分別計算出各個地區(qū)的總降雨量，最終比較13 個區(qū)域的總降雨量，得出最大總降雨量的區(qū)域。對于判斷看圖識降雨量的方式精確問題，我們先分析引起產生誤差的因素

4、，再次，根據題意我們先使用卡爾曼濾波模型計算出的總降雨量V準為基準參量，再利用雷達回波圖計算出的總降雨量 V與V準的比值Q為估算的準確程度。關鍵詞：總降雨量基本反射率灰度圖 插值擬合卡爾曼濾波一、問題重述2015年5月7日，貴州省迎來的一場大范圍降雨天氣，此次降雨天氣過程 從00:00開始至11:00結束，整個過程被雷達記錄了下來（每半小時一次），基 本上反映了這次降雨天氣的整個過程。再此基礎上計算出貴陽市逐小時降雨量及 總降雨量。并計算此次降雨過程中降雨量最大的地方可能會出現在哪個區(qū)域； 最 終對這種看圖識降雨量的方式準確程度做出判斷估算。二、問題分析對于問題一，針對貴陽市逐小時降雨量及總降

5、雨量的求解，在雷達回波圖能正確反映降雨量的條件下以面雨量大小在時間上的累積進行計算。 所以必須求出貴陽市面積 以及降雨量情況分布面積。運用 MATLAB軟件將預處理之后的雷達回波圖轉化 為灰度圖，將之轉化成矩陣以知道每個像素點的數值，再通過數值統(tǒng)計計算出每 種降雨量程度的分布面積，再進行疊加求和。而對于求解貴陽市逐小時降雨量由 于給出的雷達回波圖是每間隔半小時監(jiān)測出的數據， 因此則還需進一步的進行預 測處理。而對問題二而言，貴陽市也處于監(jiān)測范圍內，它也可能是總降雨量最大的區(qū)域， 因此需在問題一的基礎上進行總降雨量最大的區(qū)域估計。其次運用插值擬合和積分方法對其余12個區(qū)域進行總降雨量求解，最終比

6、較13個區(qū)域的總降雨量來估 計總降雨量最大的地方。針對問題三，根據題意求解對看圖識 降雨量的方式準確程度判斷問題， 應先采用最佳 有效模型測算出貼合實際降雨情況的總降雨量 V準，再應用問題一二的模型求解出的總降雨量與V準進行比值推算，從而作為看圖識降雨量的準確程度估算值三、符號說明面雨量Ni特定區(qū)域面積有限元dA上雨量各基本反射率在雷達回波圖中的像素點個數hi基本反射率，單位6mm/m4Vj表示不同基本反射率的下標 各區(qū)域的總降雨量，單位m3表示區(qū)域名稱 表示逐小時間隔段 準確程度 準確程度的比較基準參量四、模型假設4.1 對問題一二的基本假設(1)雷達監(jiān)測系統(tǒng)運行正常。(2)在雷達監(jiān)測每間隔

7、半小時之內的大氣對流穩(wěn)定。(3)在125km數據范圍內各區(qū)域地形地貌相似。(4)在層狀云降雨時，不同的降雨空間分辨率相同。(5)地面物質對雷達監(jiān)測不產生反射影響。4.2 對問題三的基本假設(1)雷達監(jiān)測系統(tǒng)運行正常。(2)在125km數據范圍內各區(qū)域地形地貌相似。五、模型建立與模型求解5.1 貴陽市逐小時降雨量及總降雨量模型求解5.1.1 基本反射率下的面積整合使用相關繪圖軟件對雷達回波圖中的貴陽市進行剪輯分割， 對貴陽市區(qū)域面 積上各時間段的雷達回波圖的降水量進行樣點采集和預處理 (部分圖片見下圖1 所示)000:24:4503:03:1305:22:5009:02:30圖1貴陽部分回波圖預

8、處理然后運用MATLAB H=imread(文件名)語句將預處理之后的雷達回波圖轉 化為灰度圖，將之轉化成矩陣以知道每個像素點的數值，具體執(zhí)行程序見附錄。 再通過數值統(tǒng)計計算出每種降雨量程度的分布面積大小，再進行疊加求和，即像素點個數總和。5.1.2 面雨量的定義在分析預報水情變化時面雨量運用相當廣泛， 一般將其定義為相關流域（或 區(qū)域）各個雨量站測得的降水量（點雨量）計算出流域（或區(qū)域）平均降水量， 這種降水量稱為面雨量。在水文學中，將面雨量表示為：(1)-1R =- RdAA A式中R為面雨量，A為特定區(qū)域面積，R為有限元dA上雨量1。區(qū)域面雨量能 客觀的反映該區(qū)域降雨情況。5.1.3 逐

9、小時降雨量與總降雨量求解根據預處理之后的雷達回波圖和 MATLAB勺數值轉換分別統(tǒng)計出每種降雨量 程度的像素點個數總和設基本反射率為,降雨量為Vj立方米，假設像素點近似等于單位面積，則有關系式：Vj = Ni xhi （2）在計算貴陽市逐小時降雨量時，題目給出的附件中雷達回波圖是以每半個小 時為間隔進行監(jiān)測測量，而待解決的問題是要求以逐小時為單位時間問隔求出每 個小時內的降雨量，因此需要進行時間段修正，將每半小時之內檢測一次的降雨 量延長至每一個小時再進行降雨量累積。公式（2）變?yōu)椋簍2Vk = I N M hdt (3)此時對應各時間段進行的像素點數乘以基本反射率的積分結果就將為貴陽市逐 小

10、時降雨量。由以上式子運用 MATLAB算出貴陽市逐小時降雨量具體數值結果見表格 （2）。表格（1）基本反射率對應降水量分布表dBZ顏色降雨量45-50紅色20mm/h40-45土黃色15mm/h30-40亮黃色10mm/h25-30綠色5mm/h20-25淺綠色3mm/h10-20淺藍色1mm/h0-10藍色0mm/h貴陽市逐小時降雨量記錄時刻有效面 積 像素值紅色 像素 值土黃色 像素值亮黃 色 像素 值綠色 像素 值淺綠 色 像素 值淺藍 色 像素 值藍色像素值降雨量（mn）5.6 ;23:45:228940742723780.05525.7 ；00:24:458940129558350.

11、01465.7 ；01:09:47894017954717990.16135.7 ；01:37:548940721792660.06035.7 ；02:00:2489405831342720.21065.7 ;02:28:298940240511724599359917437223.94185.7 ;03:02:1389408112001368139518399723565.18585.7 ;03:30:18894097810351933115217378915897.41355.7 ;04:49:0389401841470781204911841233.31275.7 ;05:50:5989

12、4020291771912456852.37415.7 ;06:52:528940103112961861260972.50255.7 ;07:26:408940631180711236823761.66625.7 ;08:23:03894078260110170.18675.7 ;09:02:30894046105117500.09985.7 ;10:27:01894023210.00265.7 ;11:00:5089400總降雨里：21.63注：降雨量=三（相應顏色彳a素值分別/有效面積像素值）*相應顏色的降雨量每個時刻（逐時刻）降雨量已用紅體字標出，總降雨量即逐時刻降雨量之和。對表格（2

13、）貴陽市逐小時降雨量分布表進行繪圖整理如圖 2所示：圖2貴陽市逐小時降雨量分布圖5.2 最大降雨量區(qū)域評估依據題目附錄表格中給出的數據，可以假定雷達監(jiān)測的125km范圍內的清鎮(zhèn) 市、修文縣、平壩縣、花溪區(qū)、長順縣、開陽縣、惠水縣、龍里縣、貴定縣、都 勻市、麻江縣、息烽縣12個區(qū)域做為總降雨量最大的區(qū)域監(jiān)測評估。由于貴陽 市完全處在雷達監(jiān)測的范圍之類，而且在求解問題一的基礎上已經求出貴陽市的 總降雨量強度較大，由此可推斷貴陽市也有可能成為總降雨量最大的區(qū)域。因此為減小模型的求解誤差，我們將貴陽市也列如總降雨量最大的區(qū)域監(jiān)測評估之 中。至此將有13個地區(qū)參與評估。在求解貴陽市逐小時降雨量及總降雨量

14、的過程中，已經假設雷達回波圖能準確反映降雨實際降雨量，因此可以在第一問貴陽市逐小時降雨量及總降雨量的解 答基礎上采用插值擬合方法擬合計算出各個區(qū)域降雨量的擬值函數方程，再根據各個區(qū)域的預處理之后的雷達回波圖和MATLAB勺數值轉換分別統(tǒng)計出每種降雨量程度的像素點個數總和、基本反射率和降雨量三者的關系，再在整個監(jiān)測過程 中進行時間上的積分求和，列出下列關系式：Vj =Njk。出 （4）由公式（4）可計算出各個地區(qū)的總降雨量。根據公式（4）和附錄表格中的數據用 MATLAB1值擬合分別得出各個區(qū)域的 函數表達式、擬合圖形以及進行積分求和之后的總降雨量： 清鎮(zhèn)市降雨量：y1 = -0.4503t +

15、4.2067 ;當 t=10 ,得 y=-0.2963當 y=0,得 t=9.3420s所以降雨量 s=0.5* （9.3420*4.2067-0.6580*0.2963 ） =19.5520mm6修文縣降雨量：y2 =-0.2036t 2.0600當 t=10, y=0.024故 s=0.5* (2.0600+0.024 ) *10=10.4200mm平壩縣降雨量：y3 = -0.6158t 5.4867當 t=10s, y=-0.6713當 y=0, t=8.91s故 s=0.5*(5.4867*8.91-1.09*0.6713)=24.0774mm16M平壩B時間花溪區(qū)降雨量：y4 u-

16、0.3285t 4.0067當 t=10 , y=0.7217,故 s=0.5* (0.7217+4.0067 ) *10=23.6420mm長順縣降雨量：y5 = -0.4612t 7.1867當 t=10, y=2.5747,故 s=0.5* (7.1867+2.5747 ) *10=48.807mm開陽縣降雨量：y6 =-0.0285t 0.3267當 t=10, y=0.0417故 s=0.5* (0.3267+0.0417 ) *10=1.8420mm惠水縣降雨量：y7 = -0.6503t 8.1867當 t=10, y=1.6837所以 s=0.5* (1.6837+8.1867

17、 ) *10=49.352mm35 11+130 -龍里縣降雨量:y8 =0.1400t 2.8000當 t=10 , y=1.4所以 s=0.5* (1.4+2.8 ) *10=21mm龍里縣貴定縣降雨量：y9 = -0.0485t 0.7333當 t=10, y=0.2483所以 s=0.5* (0.2483+0.7333 ) *10=4.908mm貴定縣0.5時間uM1.5都均市降雨量：y10 -0.3606t 0.2867當 t=10, y=-3.3193當 y=0, t=0.7950所以 s=-0.5* (0.2867*0.7950-9,205*3.3193 ) =15.1631mm

18、麻江縣降雨量：y11 =-0.3703t 0.0333當 t=10, y=-3.6697當 y=0, t=0.09所以 s=-0.5 (0.09*0.0333-9.91*3.6697 ) =18.1819mm息烽縣降雨量:y12 =-0.0200t 0.1800當 t=10, y=-0.02當 y=0, t=9所以 s=0.5* (9*0.18-1*0.02 ) =0.8mm貴陽市降雨量：y13 = -0.1308t 2.5794當 t=10 , y=1.2714所以 s=0.5* (2.5794+1.2714 ) *10=19.254mm綜上所述，進行最終總降雨量比較得降雨量最大的地區(qū)是惠水

19、縣。在雷達監(jiān)測圖中的位置用黑灰色標示出，如圖（3）所示：圖（3）總降雨量最大區(qū)域圖5.3 看圖識降雨量的方式準確程度判斷根據Pessoa等用不同空間分辨率的雷達降雨研究結論，不同的降雨空間 分辨率對于模擬的精度有很大的影響，降雨積累的時間越短，結果的誤差越大， 相對長時間的降雨累計可以消除由于短時間降雨累計的誤差，然而降雨的空間分辨率并不是越小越好。采用多大的時間和空間分辨率更合理， 與所選的模型、流 域特征、降雨的類型、降雨的強度以及流域的尺度等因素有關，并且大氣對流程 度以及水分蒸發(fā)對雷達監(jiān)測的基本反射率都存在一定影響， 其次，看圖識降雨量 主要是依據雷達監(jiān)測的基本反射率來判斷， 方式簡單

20、直觀，但基本反射率對降水 粒子直徑有所要求，且影響雷達監(jiān)測的因素較多，因此看圖識降雨量的方式必定 存在一定誤差。查閱資料可知采用卡爾曼濾波算法得出的監(jiān)測結果更加準確，更加契合實際降雨量情況，因此我們近似將卡爾曼濾波算法得出的總降雨量V準作為準確精度的比較基準參量，設看圖識降雨量的方式計算出的總降雨量為 V，準確程度為Q,則:-V 0,、Q= 1000%(5)V準由公式（5）可算出看圖識降雨量的方式的準確程度。六、模型結果分析與模型優(yōu)缺點在本次模型建立中，對于問題一貴陽市逐小時降雨量和總降雨量的求解， 模 型相對簡單，計算比較容易也容易理解，利用降雨量在時間和空間上的分段累積 求出的結果大體接近

21、真實值，但仍然存在一定誤差，不能真正準確的計算出實際 降雨量值。而對于問題二的求解過程中，本模型是基于問題一的模型上進行進一 步的優(yōu)化改進，即插值擬合與積分累加的聯(lián)合求解， 使總降水量的計算過程更加 方便簡潔。不足之處就是在于模型的建立必須依賴假設條件的多樣性，模型成立的要求 條件較多，例如要求地形相似，忽略大氣對流及水汽蒸發(fā)等，然而這些因素在實 際情況中是仍然存在的，所以模型求解出來的誤差相對較大。 而針對問題三，由 于雷達監(jiān)測站網不夠多，數據量小，采用卡爾曼濾波算法不能完整求出精確的總 降雨量，因此看圖識降雨量的方式的準確程度只能是一個估量值。七、改進方向以上模型結果分析與模型優(yōu)缺點可清楚

22、的知道，此次模型的建立條件假設約 束太多與實際上的總降雨量存在較大的出入， 所以要想減小模型結果誤差需要對 模型進行優(yōu)化改進。在模型求解面雨量及總降雨量過程中，不僅僅只用MATLABI圖片轉化為灰度圖，再根據像素點總數乘以基本反射率， 也應結合附錄表格中給出的對貴州省 12個區(qū)域的監(jiān)測數據運用 Voronoi圖網進行區(qū)域降水面積求解，再根據相應的權 重比比較兩種方法哪一種測出的總降雨量與實際降雨相近，誤差更小，更合理， 由此確定最佳有效模型。而看圖識降雨量方式的準確程度判斷應在最佳有效模型求解的基礎上進行 影響檢測因素的權重分析，再按影響檢測結果的因素所占權重比對看圖識降雨量 方式及最佳有效模

23、型進行評估對比，進一步估算出準確程度。八、 參考文獻1 葛徽衍，張永紅 .流域面雨量的一種估算方法-網絡插值法J. 陜西氣象，2006， （ 1） ： 34-35.2Pessoa M L ， Bras R L ， Williams E R Use of weather radar for flood forecasting in the sieve river basin： A sensitivity analysisJ J Appl Meteor， 1993,32： 462-4753司守奎，孫璽菁，張德存，周剛，韓慶龍，數學建模算法與應用.北京：國防工業(yè)出版社， 2013.1九、 附錄1.1

24、 雷達回波圖轉化為灰度圖H=imread（ 文件名 ）1.2 貴陽市逐小時降雨量x=1:1:12;7.4135 3.3127 2.3741 2.5025 1y=0.05520.0146 0.0603 3.94186662 0.1867 0.0998 0.0026;plot（t,y）title（ 貴陽市逐小時降雨量）xlabel（ 時間 ） ylabel（ 降雨量 ）2.1 第二問代碼 x1=1:1:10;y1=0.1 6.1 5.4 1.1 2.5 2.1 0 0 0 0;plot（x1,y1,*） xlabel（ 時間 ） ylabel（ 降雨量 ） title（ 清鎮(zhèn)市逐小時降雨量） a=

25、polyfit（x1,y1,1） a =-0.45034.2067 x2=1:1:10; y2=0 0.9 5.6 0.5 0.1 2.3 0 0 0 0; plot（x2,y2,*） xlabel（ 時間） ylabel（ 降雨量 ） title（ 修文縣） a=polyfit（x2,y2,1） x3=1:1:10; y3=0 6.2 10.8 0.4 3.3 0.3 0 0 0 0; plot（x3,y3,*） xlabel（ 時間 ） ylabel（ 降雨量 ） title（ 平壩縣 ） a=polyfit（x3,y3,1）a =-0.61585.4867 x4=1:1:10;y4=0

26、0 4.4 12.3 0.8 4 0.5 0 0 0;plot（x4,y4,*）xlabel（ 時間）ylabel（ 降雨量 ）title（ 花溪區(qū)）a=polyfit（x4,y4,1）a =-0.32854.0067 x5=1:1:10;y5=0 0 1 28.8 0.7 16 0 0 0 0;plot（x5,y5,*）xlabel（ 時間）ylabel（ 降雨量 ）title（ 長順縣）a=polyfit（x5,y5,1）a =-0.46127.1867 x6=1:1:10;y6=0 0 0.9 0.4 0.1 0 0.2 0.1 0 0;plot（x6,y6,*）xlabel（ 時間）ylabel（ 降雨量 ）title（ 開陽縣）a=polyfit（x6,y6,1）-0.02850.3267 x7=1:1:10;y7=0 0 0 34.8 7.8 2.8 0.7 0 0 0;plot（x7,y7,*）xlabel（ 時間）ylabel（ 降雨量 ）