變差函數(shù)理論在濕季降水量推算全年降水量中的應(yīng)用

趙楠楠 楊振京 楊慶華

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

關(guān)鍵詞 時間變差函數(shù);降水量;克里金插值;降水特征;濕季

中圖分類號 P 331 ?文獻標(biāo)識碼 A

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

doi：10.3969/j.issn.0517-6611.2019.12.062

開放科學(xué)（資源服務(wù)）標(biāo)識碼（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

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

1 資料與方法

1.1 數(shù)據(jù)來源

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

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

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

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

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

根據(jù)變差函數(shù)理論，當(dāng)f=0時，γ（f）的取值最小，兩者之間的相關(guān)性最強;在f趨于變程（α）的過程中γ（f）值逐漸增大，兩者之間大的相關(guān)性隨之減小;當(dāng)f>α?xí)r，γ（f）值不再改變，兩者之間沒有相關(guān)性。

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

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

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

γ（f）=C0f=0

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

C0+Cf>α（4）

式中，C0為塊金效應(yīng)，C為基臺值，C0+C為總基臺值。

1.3 時間變差函數(shù)方法的步驟與原則

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

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

2 石家莊地區(qū)年降水特征

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

3 時間變差函數(shù)方法在石家莊地區(qū)的應(yīng)用

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

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

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

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

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

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

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

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

3.4 年降水量計算結(jié)果及分析

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

參考文獻

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

[2] 張秀娟，陳曉光，王堯，等.西北四省區(qū)降水的時空變化特征分析[J].安徽農(nóng)業(yè)科學(xué)，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] 徐利崗，杜歷，姚海嬌，等.中亞干旱區(qū)降水時空變化特征及趨勢分析[J].干旱區(qū)資源與環(huán)境，2015，29（11）：121-127.

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

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

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

[8] 吳慧，吳勝安.近48年海南省極端降水時空變化趨勢[J].安徽農(nóng)業(yè)科學(xué)，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] 王剛，嚴(yán)登華，張冬冬，等.海河流域1961年-2010年極端氣溫與降水變化趨勢分析[J].南水北調(diào)與水利科技，2014，12（1）：1-6，11.

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

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

[13] 李巍，范文義，毛學(xué)剛，等.降雨量空間插值方法比較研究[J].安徽農(nóng)業(yè)科學(xué)，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] 裴韜，周成虎，李全林，等.基于變差函數(shù)分析的地震時間相關(guān)性定量估算[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] 莊立偉，王石立.東北地區(qū)逐日氣象要素的空間插值方法應(yīng)用研究[J].應(yīng)用氣象學(xué)報，2003，14（5）：605-615.