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

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

引 言

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

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

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

1 研究區(qū)域和數(shù)據(jù)

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

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

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

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

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

2.1 主要參數(shù)指標

2.1.1 降水相似度

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

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

2.1.2 地形相似度

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

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

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

2.2 最佳對比區(qū)確定

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

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

3 不同統(tǒng)計方法的自然雨量關系模型對比

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

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

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

圖4為不同雨強(I)的樣本量。由圖4可知,小時雨量主要集中在雨強小值區(qū)域。將降水分為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 不同雨強樣本量統(tǒng)計

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

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

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

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

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

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

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

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

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

圖9 3種優(yōu)化器降水估測對比

4 結論與討論

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

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

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

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

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