變差函數理論在濕季降水量推算全年降水量中的應用

趙楠楠 楊振京 楊慶華

摘要 基于石家莊市16個氣象站點1981—2010年逐日降水資料，根據變差函數理論分析并獲取了時間變差函數參數，利用該參數作為克里金插值的約束信息推算全年降水量。結果表明，濕季（5—9月）推算結果明顯優于其他時段，且具有較強的穩定性;該時段的推算誤差率普遍小于11%，誤差率的平均值為-1.24%，可作為推估全年降水量的有效方法。

關鍵詞 時間變差函數;降水量;克里金插值;降水特征;濕季

中圖分類號 P 331 ?文獻標識碼 A

文章編號 0517-6611（2019）12-0225-06

doi：10.3969/j.issn.0517-6611.2019.12.062

開放科學（資源服務）標識碼（OSID）：

Abstract Based on the daily rainfall data at 16 meteorological stations in Shijiazhuang from 1981 to 2010，the semivariogram was analyzed to calculate the variogram parameters （range;sill and nugget），and then the annual precipitation was estimated by Ordinary Kriging（OK） interpolation constrained by these parameters.The results showed that the estimated precipitation of wet season was closer to the actual precipitation compared to others periods and more stability.The calculated error rates of test data in Shijiazhuang were normally less than 11% and the mean error rate was 1.24%.It indicated that time semivariogram was an useful method to estimate annual precipitation under the condition of wet seasonal precipitation was known.

Key words Time semivariogram;Precipitation;Kriging interpolation;Precipitation characteristics;Wet season

連續可靠的降水數據是水文分析、水資源管理和氣候研究的重要前提。但是，受人為、經濟和氣候等因素影響難以獲取連續、完整的降水數據，因此通過已有數據推算補全缺失資料就顯得十分必要。氣候要素本身就是一個隨機變量[1]，統計學方法是氣候學家們用來研究氣候要素常用的方法。利用統計學方法研究降水的文獻很多，主要是針對其時空分布[2-4]、頻數[5]、降水周期[6]、強度特征[5]和變化趨勢[7-8]、預估[9]等。降水資料缺失主要分為2個方面：空間資料缺失，如王剛等[10]通過與該站臨近的、相關性高的站點補齊缺失數據;時序資料缺失，張強等[11]用三次樣條函數內插方式補齊缺失數據;晏利斌[7]利用多元回歸插值方法補齊缺失數據;張志萍等[12]則通過雨量站分區間利用相關站點汛期降水量占全年比率插補全年降水量。通常，氣象站點有完整的濕季降水資料而缺少全年完整數據，如張志萍等[12]曾指出大理河流域許多雨量站只記錄了5—10月份的降水資料，顯然其不能代表年降水量，所以需要通過濕季降水數據擴推全年降水量。然而，對于年內缺失濕季之外降水推算全年降水量的研究很少。變差函數是研究變量空間變化極為有利的工具，是地質統計學的基本內容之一，利用地質統計學研究降水[13-14]、水位[15]等水文參數空間變異特征取得了較好的研究效果。因此，筆者選取石家莊16個氣象站點的累年逐日降水資料，分析、獲取各站點時間變差函數參數，以這些參數為插值約束進而通過濕季降水量推算全年降水量。

1 資料與方法

1.1 數據來源

中國氣象局氣象數據中心開放全國國家級站點地面氣象資料，覆蓋了絕大多數的縣級區域。該研究整理和分析了經中國氣象局檢驗修正后（包括16個縣站點1981—2010年累年逐日降水量、石家莊站點1955—2015年逐月降水量）的地面氣象資料，這些資料具有可靠性高、時間連續性好等特點，得到廣泛認可[10]。

1.2時間變差函數 變差函數是用來研究區域變量的空間變化特征的有力工具，根據變差函數的研究理論，把變差函數中的位置和距離替換為時間和時間差，則可以用來分析變量的時間相關性，稱其為時間變差函數[16]，時間變差函數γ（t，f）可表達為：

γ（t，f）=12E[Z（t）-Z（t+f）]2（1）

式中，t代表時間，f代表隨機過程中的時間間隔，Z（t）為t時刻隨機變量的取值，E[]為隨機變量的期望。如果假設隨機過程是二階平穩的，時間變差函數表達為：

γ（f）=12E[Z（t）-Z（t+f）]2（2）

根據變差函數理論，當f=0時，γ（f）的取值最小，兩者之間的相關性最強;在f趨于變程（α）的過程中γ（f）值逐漸增大，兩者之間大的相關性隨之減小;當f>α時，γ（f）值不再改變，兩者之間沒有相關性。

由于受實際取樣有限等因素的影響，需先制作試驗變差函數，試驗變差函數計算公式為：

γ*（f）=12N（f）N（f）i=1E[Z（ti）-Z（ti+f）]2（3）

式中，N（f）是f間隔下樣本點對的數目。試驗變差函數確定之后，通過一種合適的理論變差函數對其進行擬合，從而形成理論變差函數模型。常用的變差函數理論模型有球狀模型、高斯模型、指數模型、線性模型等，由于球狀模型可以較好地擬合各站點試驗變差函數計算結果，故該研究使用的理論模型為球狀模型，其表達式為：

γ（f）=C0f=0

C0+C（32fα-12f3α3）0

C0+Cf>α（4）

式中，C0為塊金效應，C為基臺值，C0+C為總基臺值。

1.3 時間變差函數方法的步驟與原則

①對研究區周邊地區站點的累年逐日降水資料進行整理與分析，獲得研究區相關站點的時間變差函數各項參數。如果研究區有氣象站點，可直接使用該站點通過分析后的變差函數參數;如果沒有站點，則可以通過周邊站點插值后獲取。

②獲取累年逐日最大降水量（Pday）以及累年逐月最大降水量（Pmonth）。分析研究區濕季各月降水量，超過Pmonth降水量的月份把超出降水量分配給鄰月，優先分配給鄰月降水量大的月份;月間降水量分配完畢之后，進行月內降水量分配，月內每天的降水量以Pday優先分配，不能以整數天分配Pday的，Pday分配完后，最后一天分配剩余降水量，把這些分配降水的天放置在該月中間。③利用第一步獲得的時間變差函數各項參數作為約束，對濕季的日降水量進行普通克里金插值，最后提取該年每天降水量，求和后即可得到該年的全年降水量。

2 石家莊地區年降水特征

石家莊氣象站點在20世紀50年代中期就有連續且完整的年降水資料，其中5—9月多年平均降水量占全年降水量的80%以上，因此該研究中把該時段定義為研究區的濕季。選取石家莊站點的1955—2015年的逐月降水資料，計算每年濕季降水量占該年降水總量的比率（圖1）。從圖1可以看出，近61年石家莊地區年降水量變化幅度較大，其中1996年降水量最大，高達1 097.1 mm，1972年降水量最小，僅226.0 mm，兩者數量懸殊相差甚大。從比率曲線可以看出，石家莊濕季降水量（Pwet）占該年降水量的比率變化范圍也較大，在1968年比率最小，僅55.04%，在1988年比率最大，達97.72%，幾乎全年降水量都在濕季完成，同時也反映了年內降水量的分配也是極不均勻的。此外，降水量最多（最少）的年份，其占全年降水比率并不是最大（最小），且每年比率各不相同，因此常規的方法難以推算該區年降水量。

3 時間變差函數方法在石家莊地區的應用

3.1 石家莊市氣象站點累年逐日降水資料相關性分析

通過對收集到的石家莊16個縣站點累年逐月、累年逐日降水數據統計與分析，發現累年逐月、累年逐日數據通過球狀模型都可以很好地擬合原始數據。但是，經進一步研究（以行唐為例），發現累年逐月降水數據規律性較差（圖2），雖然根據原始數據可以較好地擬合變差函數，但其可能是由于數據點太少而造成的“假象”，如果將其應用到后期的插值過程中，可能會產生較大的誤差，從而導致錯誤的結果，因此該研究僅對累年逐日降水數據的時間變差函數進行分析與利用。

石家莊16個縣站點累年逐日降水數據都符合正態分布規律，滿足高斯域變差函數基本要求。各站點試驗時間變差函數的擬合結果見圖3（一維隨機變量），通過球模型可以很好地擬合原始數據（R2均在0.94以上），擬合出的各項參數見表1。從圖3和表1可以看出，各站點之間的塊金效應不同，變化范圍較大，為0.54～1.40，相差近2倍;基臺值為5.124～7.976，變程為136.4～144.2，它們主要受控于累年日降水量多少及其時間分布。平山站點的塊金效應最大（1.40），但是其變程不是最大的;同樣，行唐的塊金效應最小（0.54），其變程并不是最小的。上述情況可以用變差函數性質解釋：變程越大，隨機變量的非均質性越小，但塊金效應增大會在一定程度上增加隨機變量的非均質性。

3.2 石家莊站點各項時間變差函數參數的確定

時間變差函數的各項參數與日降水量及其時間分布有密切的聯系，而降水量的空間分布通常可以由插值方法確定，在氣象數據稀少等特殊區域，時間變差函數的各項參數也可以用降水量空間插值方法確定。克里金插值、距離反比插值是降水空間插值中最常用的方法。許多學者研究發現，克里金插值比距離反比插值精度更高，結果更可靠[17]。但是在對石家莊16個氣象站點的時間變差函數的塊金效應、基臺值和變程進行變差函數擬合時，發現這些參數空間變異性規律極差（圖4），如果使用克里金插值可能會產生較大的誤差。另外，莊立偉等[18]在研究東北地區逐日降水空間插值方法時發現，距離反比插值精度比克里金插值精度高，但平滑程度較小，更適合日降水量的空間插值，因此石家莊站點的各項參數由距離反比方法確定，其站點各項參數經插值后結果見表1。

3.3 石家莊站點Pday和Pmonth的確定和分配

石家莊站點時間變差函數各項參數獲取后，需要確定該站點的Pday和Pmonth。該站點距離藁城和正定2個站點較近，其大致位于2站的中間位置，因此石家莊站點的Pday和Pmonth由2處的平均值確定（表1）。Pday和Pmonth確定完以后開始進行月間與月內分配。下面以石家莊2016年的Pwet為例（表2）進行分配，具體步驟如下：7月降水量445.6 mm遠大于Pmonth（127.0 mm），故7月降水量分配127.0 mm，剩余部分優先分配鄰月降水量大的月份，依次分配給5、6、8、9月。分配完成后的各月降水量見表2。

月間降水量分配完畢后，進行月內降水量分配，其中，石家莊站點Pday為9.45 mm。以5月為例優先分配9.45 mm，其中88.4/9.45≈9.35，則連續分配9天Pday，第10天分配剩余降水量3.35 mm，即5月的11—19日日降水量均為9.45 mm，20日為3.35 mm。Pday分配完畢后，把分到降水量的10 d放在該月中間位置，其余日期降水量分配為0mm，另外4個月份做法相同。濕季之外的月份降水量設為未知，通過后期的插值獲得，到此，降水量分配完畢。

3.4 年降水量計算結果及分析

收集到了石家莊站點2011—2016年和2015年正定、井陘等6個站點的逐月降水數據進行驗證上述方法的可行性。首先，把各個站點按照上述的步驟處理完畢之后導入到ArcGIS 10.2.2軟件中，然后，在表1中對應站點時間變差函數參數約束下進行普通克里金插值，處理范圍為（0～366，-1～1），最后，提取該年內相應天數的插值數值求和后即可得到該年的年降水值。測驗數據的原始相關信息見表3，為證明濕季推算結果可靠性，提供了相鄰不同時段的計算結果（表4）。

從表3可以看出，同一站點不同年份以及同一年份不同站點之間Pwet、Pa和月降水量差別明顯，這符合石家莊地區比率、年降水量變化幅度較大以及年內降水分布不均的降水特征。

誤差率是推算精度的重要指標，從圖5可以看出，3—7、6—10和7—11月的計算降水量明顯的偏離Pa，它們的推算結果普遍難以接受（表4），不能獲得理想的推算效果。雖然4—8月推算結果在特殊情況下產生比5—9月份較優的推算結果，但是4—8月份的推算結果具有很大的不穩定性，最大誤差達37.68%，因此，如果采用4—8月降水資料推算Pa，其推算結果可能很大程度偏離Pa，無法獲取具有參考價值的推算結果。5—9月份推算結果則相對穩定，雖然可能產生高達16.91%的誤差，但是其推算結果一直游離在Pa附近，具有較高的參考價值。

利用5—9月降水資料推算實際降水量（Pa）在石家莊地區獲得了較好的應用，其計算誤差率普遍在11%以下，其中石家莊2013年計算誤差率僅為2.11%，說明該方法可以作為一種全年降水量的推算方法;但石家莊2011年和無極2015年推算誤差卻分別達16.91%和-12.35%（表4）。誤差較大的原因可能是由于Pmonth和Pday不能代表站點自身降水特征以及克里金的圓滑效應引起的。

由于無法收集充足的歷史數據確定去掉奇異值后的Pmonth和Pday，對無極和平山借用相鄰站點的Pmonth和Pday進行重新計算（表4），這2個站點的推算誤差得到了有效的縮小，故在推算前需要對Pmonth和Pday的代表性進行檢驗，盡可能地體現氣象站點自身的降水特征。

克里金插值的圓滑效應也是誤差來源之一。在占全年比率大致不變的情況下，Pwet過高或過小都會導致誤差率較大;在Pwet大致相近的情況下，占全年比率的過高或過小同樣也會導致誤差率較大（圖6）。前者在比率較大的情況下，由于克里金法的圓滑效應，降水量多或少會導致其周圍未分配降水量日期的插值結果多或少，并對誤差進行傳遞，從而造

成較大的誤差率。在Pwet大致相同的情況下也是這種原因。

如在2011年石家莊站點，Pwet占全年比率高達89.05%，降水量為600.4 mm，克里金插值的圓滑效應使得該年產生較大的誤差率;而2013年石家莊站點比率高達89.26%，但是其Pwet僅為453.7 mm，雖然克里金法仍有圓滑效應存在，但是其引起的絕對誤差較小，因而其誤差率很小。2015年藁城站點則相反，其Pwet及其比率都較小，從而產生較大的誤差。

4 結論

利用石家莊16個氣象站點的逐日降水資料，借助變差函數相關理論，通過對濕季降水量推算全年降水量方法的研究，得出以下結論：

（1）不同時段的推算結果有很大差異，其中濕季的推算結果最為可靠。3—7、6—10和7—11月推算結果較大程度地偏離實際降水量（Pa）;4—8月推算結果有較大的不穩定性;而濕季推算結果更可靠，其推算結果和誤差的穩定性均較好。

（2）時間變差函數方法在石家莊測驗數據中的推算誤差普遍小于11%，其中石家莊2011年推算誤差最小，僅為2.11%;其誤差率的平均值為-1.24%，方差為8.66，具有較高的參考價值，可以作為一種推估全年降水量的有效方法。

（3）時間變差函數方法在石家莊2011年和無極2015年測試數據中誤差較大是由于Pmonth和Pday不能代表站點自身降水特征及克里金插值的圓滑作用引起的。其中無極和平山站點借助各自臨近站點的Pmonth和Pday使得誤差得到一定程度的縮小。

參考文獻

[1]江志紅，丁裕國，陳威霖.21世紀中國極端降水事件預估[J].氣候變化研究進展，2007，3（4）：202-207.

[2] 張秀娟，陳曉光，王堯，等.西北四省區降水的時空變化特征分析[J].安徽農業科學，2012，40（18）：9809-9812.

[3] CUI L F，WANG L C，LAI Z P，et al.Innovative trend analysis of annual and seasonal air temperature and rainfall in the Yangtze River Basin;China during 1960-2015[J].Journal of atmospheric and solarterrestrial physics，2017，164：48-59.

[4] 徐利崗，杜歷，姚海嬌，等.中亞干旱區降水時空變化特征及趨勢分析[J].干旱區資源與環境，2015，29（11）：121-127.

[5] 王志福，錢永甫.中國極端降水事件的頻數和強度特征[J].水科學進展，2009，20（1）：1-9.

[6] 劉健，夏軍，王明森，等.1961—2015年山東省降水周期變化特征[J].人民黃河，2017，39（4）：6-10.

[7] 晏利斌.1961—2014 年黃土高原氣溫和降水變化趨勢[J].地球環境學報，2015，6（5）：276-282.

[8] 吳慧，吳勝安.近48年海南省極端降水時空變化趨勢[J].安徽農業科學，2010，38（19）：10101-10103.

[9] ALIZADEH M J，KAVIANPOUR M R，KISI O，et al.A new approach for simulating and forecasting the rainfallrunoff process within the next two months[J].Journal of hydrology，2017，548：588-597.

[10] 王剛，嚴登華，張冬冬，等.海河流域1961年-2010年極端氣溫與降水變化趨勢分析[J].南水北調與水利科技，2014，12（1）：1-6，11.

[11] 張強，孫鵬，陳喜，等.1956～2000 年中國地表水資源狀況：變化特征、成因及影響[J].地理科學，2011，31（12）：1430-1436.

[12] 張志萍，冉大川，慕志龍.大理河流域降水資料插補方法探討[J].人民黃河，2006（12）：26-27，78.

[13] 李巍，范文義，毛學剛，等.降雨量空間插值方法比較研究[J].安徽農業科學，2014，42（12）：3667-3669.

[14] DIRKS K N，HAY J E，STOW C D，et al.Highresolution studies of rainfall on Norfolk Island Part II： Interpolation of rainfall data[J].Journal of hydrology，1998，208（3）：187-193.

[15] OHMER M，LIESCH T，GOEPPERT N，et al.On the optimal selection of interpolation methods for groundwater contouring： An example of propagation of uncertainty regarding interaquifer exchange[J].Advances in water resources，2017，109：121-132.

[16] 裴韜，周成虎，李全林，等.基于變差函數分析的地震時間相關性定量估算[J].地震，2002，22（2）：17-21.

[17] PIAZZA A D，CONTI F L，NOTO L V，et al.Comparative analysis of different techniques for spatial interpolation of rainfall data to create a serially complete monthly time series of precipitation for Sicily，Italy[J].International journal of applied earth observations & geoinformation，2011，13（3）：396-408.

[18] 莊立偉，王石立.東北地區逐日氣象要素的空間插值方法應用研究[J].應用氣象學報，2003，14（5）：605-615.