機器學習在人工增雨效果統計檢驗中的應用

李 丹 林 文* 劉 群 馮宏芳 胡淑萍 汪智海

1)(福建省氣象科學研究所, 福州 350001) 2)(福建省災害天氣重點實驗室, 福州 350001) 3)(中國氣象局海峽災害天氣重點開放實驗室, 福州 350001) 4)(閩南師范大學數學與統計學院, 漳州 363000) 5)(福建省寧德市古田縣氣象局, 寧德 352000)

引 言

近年人工影響天氣事業不斷發展,在農業抗旱、水庫蓄水、凈化空氣等工作中發揮了重要作用,作為人工影響天氣作業的重要環節,人工增雨作業效果的科學評價受到人們高度關注[1-3]。經過大量實踐研究,國內外學者將人工增雨效果檢驗方法總結為物理檢驗、數值模擬檢驗和統計檢驗[4-15]。物理檢驗主要針對人工影響天氣作業前后云系的宏微觀變化特征,結合云降水形成及其催化原理,對比目標云系被影響后產生的相關物理響應(微物理結構或宏觀動力結構變化等),對作業產生的直接效果進行定性或定量分析[16-18];數值模擬檢驗是利用能夠描述云降水微物理過程及人工增雨催化過程的數值模式,通過改變催化條件,定量模擬出人工增雨作業后的云系發展變化和產生降水情況,并與未催化的自然發展云系觀測結果進行對比,得到作業效果[19-21];統計檢驗是對間接效果即產生的地面降水進行定量分析,主要通過不同人工增雨試驗方案,利用數理統計理論對作業效果進行定量檢驗分析[22-25]。研究表明,序列分析、雙比分析、區域對比分析、區域歷史回歸分析和基于聚類統計的浮動對比區等方法是常用的統計檢驗方法。其中,區域歷史回歸分析檢驗功效高、準確度和靈敏度較好[26-27],是國內外人工影響天氣作業效果檢驗比較推薦的效果檢驗方法,該方法基于歷史降水長序列數據,通過分析目標區和對比區的降水相關性,估算目標區作業期間的自然雨量,并與目標區實際雨量進行比較,得到的差值則視為增雨作業效果。該方法假定作業期間目標區和對比區的雨量相關關系與歷史相同天氣下兩區域雨量相關關系相同,但實際上該假定很難滿足,這是因為實際降水自然變率較大,目標區與對比區的關系缺乏穩定性,從而導致評估結果穩定性較差,這也是區域歷史回歸方案最大的困難[28-29]。如果能夠選擇恰當的對比區,且樣本足夠多,該方案的評估功效仍較高。國內外許多抗旱增雨的業務性試驗效果評估均采用該方案。此外,房彬等[30-31]和翟羽等[32]等利用聚類分析進一步改進非隨機化人工增雨效果區域歷史回歸方法,將物理協變量作為控制因子和網格插值計算雨量,提出一種新的浮動對比區歷史回歸統計檢驗方法(簡稱CA-FCM方法),并用于河南人工影響天氣作業效果分析,結果表明效果評估功效顯著提高[32]。隨著機器學習在氣象領域的不斷發展和廣泛應用,發現其能夠更好地自適應數據變化并提取特征,具有強大的非線性建模能力[33-37]。

為了盡可能客觀、定量地檢驗人工增雨作業效果,本文基于2014年1月—2023年1月福建地區降水數據和作業信息,利用機器學習和多種數學統計方法[38],對比不同回歸統計檢驗方案,旨在進一步優化基于浮動對比區的最佳自然雨量估測模型,為開展區域內人工增雨效果統計檢驗提供參考。

本文插圖中所涉及的國界和行政區域界線基于審圖號為GS(2017)3320號標準地圖制作,底圖無修改。

1 研究區域和數據

本文研究區域為福建省古田人工增雨效果檢驗隨機試驗區,覆蓋閩北地區古田、屏南、周寧和建甌等地,其范圍為26.1°～27.3°N,118.0°～119.5°E,選取該區域內古田水庫流域人工增雨效果檢驗外場試驗作業影響區域為目標區(圖1)。效果檢驗必須準確預計效果出現的區域、時間、量值等[39],因此效果評估前首先需確定效果所在的區域即目標區。

圖1 福建省古田人工增雨效果檢驗基地

本文選取2014年1月—2023年1月試驗區自動氣象站小時降水數據,結合福建省人工影響天氣作業歷史數據,根據作業點經緯度、作業起止時間以及目標云系移動范圍建立歷史數據庫并進行標記分類,期間福建省境內共開展人工影響天氣作業1.06萬余次,試驗區內開展人工增雨作業隨機試驗約80次,其中隨機試驗主要在石坑和西溪2個作業點進行火箭冷云催化作業。為排除由于人工影響天氣作業對自然降水數據的影響,將人工影響天氣作業后4 h內目標云系移動過程中覆蓋到試驗區域的自動氣象站降水數據標記為人工增雨作業影響數據庫,未覆蓋到試驗區域的降水數據則標記自然降水數據庫,并將試驗區域內3個自動氣象站雨量不小于0.1 mm的小時降水數據視為1次有效樣本,從而保證樣本量和數據質量。人工增雨作業信息數據來源于福建省氣象局人工影響天氣指揮中心。

2 基于區域歷史回歸方法確定最佳對比區

區域歷史回歸方案為選擇1個或多個與目標區天氣地理條件相似、降水相關性較好的區域為對比區,然后根據兩區域歷史雨量建立區域歷史回歸方程,將對比區的雨量代入方程求得作業區自然降水估測值,并與作業區降水實測值對比以確定增雨效果。區域歷史回歸分析要求對比區與目標區相互獨立,但天氣系統、降水分布和地形等方面的相似度較高,因此本文以降水相似度和地形相似度為主要參數指標。此外,對比區應選擇在目標區上風向,確保不受人工影響天氣作業催化劑擴散作用的影響。

2.1 主要參數指標

2.1.1 降水相似度

降水相似度決定了對比區與目標區歷史自然降水序列關系模型的可解釋性,以及基于對比區雨量對目標區自然雨量的可估測能力。雖然用日尺度以上雨量作為統計變量相對簡單且方便,但因時間跨度較大,且包含自然降水時段,不利于對影響時間有限的單次人工增雨作業進行合理準確評估。因此選擇目標區內各自動氣象站小時雨量數據,通過插值獲得該區域內小時平均面雨量,表征區域降水強度。

對自動氣象站逐時、逐日、逐候雨量進行K-S檢驗(Kolmogorov-Smirnov test)分析,自然雨量數據集不滿足正態分布特征,通過對數變換或者六次方根變換有所改善,不適合直接使用線性關系模型進行兩區域數據擬合。劉晴[24]和程鵬等[25]提出為滿足統計變量服從正態分布的要求,統計變量相對最優的選擇是候雨量或旬雨量的六次方根值。本文對小時雨量進行六次方根變換,采用線性回歸分析方法分析目標區和對比區的小時平均面雨量,將二者相關系數作為兩區域降水相似度。

2.1.2 地形相似度

地形相似度是基于地圖影像數據,利用地形特征數據劃分每種地物的分布區域,以區域為基本單位提取影像特征并進行對比。本文綜合考慮兩個區域地圖要素的形狀相似、位置相似以及信息內容相似程度,即屬性特征相似性。將圖片信息轉化為數組,通過灰度化處理簡化圖像色彩,計算所有像素的灰度平均值,通過二值化得到圖像的哈希值,比較像素的哈希值差異,利用漢明距離法得到兩區域的相似度。

由兩區域地形相似度與降水相似度擬合關系可知二者呈較明顯正相關關系(圖2)。因此,基于自然降水數據庫和GIS空間數據,利用浮動對比區方法選取對比區時,將降水相似度和地形相似度較高(面積、形狀與目標區相同)的區域確定為最佳對比區。

圖2 目標區與對比區地形相似度與降水相似度擬合關系

2.2 最佳對比區確定

古田試驗區的主要天氣系統為低渦切變、暖區輻合和高空槽,降水云系多為向東北方向移動的積層混合云[40]。由圖3可知,紅色方框為古田隨機試驗人工影響天氣目標區,結合歷史天氣類型及云系主要移動路徑,在目標區上游和側方設計連續多個形狀大小一樣的區域為浮動對比區(藍色方框)?；跉v史自然降水數據,將20個浮動對比區與目標區進行面雨量相關性和地形相似度分析?？紤]到地形變化對天氣系統的影響因素,本文優先選擇相關系數最高的區域,若相關系數相同,則優先選擇地形相似度最高的區域作為最佳對比區(圖3黑色方框),兩區雨量相關系數為0.63,地形相似度為53.08%。

圖3 古田隨機試驗目標區、浮動對比區與最佳對比區設置

3 不同統計方法的自然雨量關系模型對比

區域歷史回歸建立兩區域雨量關系模型的目的是利用對比區的自然雨量合理預期目標區作業期的自然雨量,并將其與作業影響后目標區的降水實測值對比,確定增雨效果。通過引入機器學習,對比不同統計方法的自然雨量關系模型,基于最優雨量關系模型得到目標區最接近自然降水的降水估測值,進而得到人工增雨作業后的合理增雨量。

下文采用線性擬合、多項式回歸、樣條回歸、機器學習進行對比分析[41-42]。選取均方根誤差和確定系數比較不同擬合方法的優劣。均方根誤差是回歸模型的擬合標準差,越接近于0,模型預測結果越精準;確定系數反映因變量的變化能由自變量解釋的比例,表征回歸模型的可靠程度,其正常取值范圍為[0,1],該數值越大,代表模型解釋能力越強,對數據預測效果更好。

3.1 基于平均面雨量的多種統計關系模型對比評估

圖4為不同雨強(I)的樣本量。由圖4可知,小時雨量主要集中在雨強小值區域。將降水分為4類:弱降水(0.1 mm·h-1≤I<5 mm·h-1)、一般降水(5 mm·h-1≤I<10 mm·h-1)、中等降水(10 mm·h-1≤I<25 mm·h-1)和強降水(I≥25 mm·h-1),其中弱降水樣本量占比為95.98%,強降水占比僅為0.03%(表1)。

表1 2014—2023年不同等級降水樣本量

圖4 不同雨強樣本量統計

分別選取對比區小時平均面雨量和平均面雨量六次方根變換值作為預測變量,選取目標區小時平均面雨量和平均面雨量六次方根變換值作為響應變量,利用線性回歸、多項式回歸(二項、三項、四項、五項)和樣條回歸(回歸次數分別為1、2、3、4)多種模型,對比不同降水等級樣本和總樣本數據的擬合結果,建立的回歸模型均方根誤差如圖5所示。由圖5可知,選取平均面雨量作為預測變量,一般降水和中等降水樣本均方根誤差為3.3～4.5 mm,總樣本的均方根誤差相對較小,約為1.1 mm,其次為弱降水數據。平均面雨量經六次方根變換后,總樣本的統計結果同樣表現相對較好,均方根誤差相對較小。在后續分析中可以基于總樣本進行統計分析,其中針對總樣本構建相關模型時,四項式回歸的均方根誤差最小,其次為二次樣條回歸。

圖5 不同雨強下各模型的降水估測均方根誤差對比

圖6 基于總樣本的不同模型降水估測的均方根誤差和確定系數

3.2 基于面雨量空間格點的卷積神經網絡機器學習關系模型

考慮到兩區域面積較大,降水的空間分布明顯不均,單純以區域平均面雨量代表區域降水存在局限性,且兩區域的上下游關系導致降水存在明顯時間序列效應,同一降水過程目標區和對比區的歐式距離最近。因此嘗試由點到點轉化為面到面,利用CNN機器學習方法研究兩區域面雨量關系,該模型能夠涵蓋云系的移向和移速,即降低歐氏距離的影響因素。如圖7所示,目標區有33個自動氣象站,對比區有18個自動氣象站。對目標區和對比區進行空間均勻格點化,并利用反距離權重、克里金、線性和三次樣條4種插值方法計算區域內面雨量空間分布。以2016年4月9日19:00(北京時,下同)為例,4種方法插值結果如圖8所示。由圖8可見,克里金插值法導致分析區域內18.9 mm的強降水中心丟失,整體面雨量分布偏小;線性和樣條插值均存在邊界效應,邊界處的數據極易不穩定、失真;反距離權重方法能相對準確地反映降水強度空間分布特征,插值效果相對較好。

圖7 目標區和對比區格點化及區域自動氣象站分布

圖8 2016年4月9日19:00不同插值方法得到的面雨量分布

基于面雨量空間格點數據,利用CNN建立對比區-目標區雨量最佳關系模型。CNN主要適用于解決圖像處理、識別問題,無需事先提取特征,將圖像的原始像素直接作為輸入,大幅減少了傳統方法所需的大量重復、繁瑣的數據預處理工作。此外,卷積的局部感知、時空亞采樣和權值共享結構能夠大量減少需要訓練的參數數量,避免過擬合且可以降低模型復雜度。各輸入要素插值到格點上后用CNN提取特征,能更好地保留區域內云系生消發展過程中的移動、變化等特征,提高空間預測能力。

采用自適應矩估計、均方根傳遞和梯度隨機下降3種不同優化器訓練模型,通過調整過濾器和卷積核尺寸建立兩區域小時面雨量格點數據關系模型。其中,自適應矩估計(簡稱ADAM)優化器是一種能夠自適應地計算并調節每個參數學習率的方法,在均方根傳遞(簡稱RMSP)的基礎上結合動量梯度下降的方法,實現快速收斂;梯度隨機下降(簡稱SGD)優化器每次迭代隨機選擇1個樣本進行訓練,在樣本量較大的情況下,能夠極大加速每輪參數的更新速度。此處定義模型的均方根誤差為整個區域面雨量格點值均方根誤差的平均值。經過比較,ADAM優化器適合該模型,其均方根誤差為0.61 mm,其次為SGD優化器,其均方根誤差為0.67 mm。隨機選取2019年1月20日20:00,2019年5月5日22:00和2022年3月2日03:00共3個個例,比較3種優化器的性能(圖9)。由圖9可知,個例1(2019年1月20日20:00) RMSP優化器估測的降水分布和強度偏差最大,雨量大于0.7 mm的區域面積明顯偏大;個例2(2019年5月5日22:00) RMSP和SGD兩種優化器得到的目標區降水估測值與實測值偏差亦相對較大,且RMSP優化器估測的雨強明顯偏弱;個例3(2022年3月2日03:00) RMSP和SGD 兩種優化器估測的降水估測值均偏大,且雨量大于0.6 mm的區域面積偏大。綜上,基于ADAM優化器的目標區降水估測值與實測值較一致。

圖9 3種優化器降水估測對比

4 結論與討論

為進一步優化人工增雨作業效果統計檢驗方法,本文基于2014年1月—2023年1月小時降水數據,利用機器學習,結合線性擬合、多項式回歸和樣條回歸等多種數學統計方法,建立目標區和對比區間不同雨量關系模型并進行對比分析,得到如下結論:

1) 古田隨機試驗目標區與其周圍浮動對比區地形相似度和降水相似度呈較明顯的正相關關系,結合降水云系主要移動路徑,在目標區上游位置得到最佳對比區,雨量相關系數為0.63,地形相似度為53.08%。

2) 歷史降水數據樣本主要集中在弱降水等級,利用多種回歸方法(線性回歸、多項式回歸和樣條回歸),采用分段雨強建立兩區域間的線性回歸模型改進不明顯,其均方根誤差普遍高于總樣本的統計結果。

4) 基于面雨量空間格點數據,采用CNN 3種不同優化器(RMSP、ADAM和SGD)建立對比區-目標區雨量關系模型,其中ADAM優化器模型最優,利用對比區估算目標區自然雨量能力最強,其均方根誤差為0.61 mm。

分析表明:CNN機器學習方法更適合于建立兩區域間的雨量關系模型,且基于面雨量空間格點的機器學習雨量關系模型可在一定程度上消除平均面雨量進行六次方根變換時產生的量綱影響,并減少強降水中心的擾動。作為人工增雨作業效果定量評估的關鍵手段,可根據典型天氣形勢或主要移動路徑分類細化模型,進一步提高目標區與目標雨量關系模型的可信度和準確度。