趙楠楠 楊振京 楊慶華

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

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

中圖分類號(hào) P 331 ?文獻(xiàn)標(biāo)識(shí)碼 A

文章編號(hào) 0517-6611（2019）12-0225-06

doi：10.3969/j.issn.0517-6611.2019.12.062

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

1 資料與方法

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

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

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

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

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

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

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

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

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

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

γ（f）=C0f=0

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

C0+Cf>α（4）

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

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

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

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

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

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

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

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

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

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

3.2 石家莊站點(diǎn)各項(xiàng)時(shí)間變差函數(shù)參數(shù)的確定

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

3.3 石家莊站點(diǎn)Pday和Pmonth的確定和分配

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

月間降水量分配完畢后，進(jìn)行月內(nèi)降水量分配，其中，石家莊站點(diǎn)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個(gè)月份做法相同。濕季之外的月份降水量設(shè)為未知，通過后期的插值獲得，到此，降水量分配完畢。

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

參考文獻(xiàn)

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

[2] 張秀娟，陳曉光，王堯，等.西北四省區(qū)降水的時(shí)空變化特征分析[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ū)降水時(shí)空變化特征及趨勢(shì)分析[J].干旱區(qū)資源與環(huán)境，2015，29（11）：121-127.

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

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

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

[8] 吳慧，吳勝安.近48年海南省極端降水時(shí)空變化趨勢(shì)[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年極端氣溫與降水變化趨勢(shì)分析[J].南水北調(diào)與水利科技，2014，12（1）：1-6，11.

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

[12] 張志萍，冉大川，慕志龍.大理河流域降水資料插補(bǔ)方法探討[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ù)分析的地震時(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é)報(bào)，2003，14（5）：605-615.