基于半方差函數(shù)和Anusplin方法華北地區(qū)降水量空間分布研究

王克磊，曲華斌，王 爽

(山東省海河淮河小清河流域水利管理服務(wù)中心，山東 濟(jì)南 250100)

降水是自然界水循環(huán)過程中最活躍的環(huán)節(jié)之一，對農(nóng)業(yè)生產(chǎn)、水利建設(shè)、生態(tài)健康等有重要影響[1-4]。降水量時(shí)空分布模式研究是水文氣象研究中的一項(xiàng)關(guān)鍵問題，而以半方差函數(shù)為核心地統(tǒng)計(jì)學(xué)是以空間變量為對象的統(tǒng)計(jì)學(xué)分支，可對降水要素進(jìn)行空間分析和模擬，并在析取空間要素半方差結(jié)構(gòu)方面得到廣泛應(yīng)用[2-5]。Anusplin插值法是專門針對氣象要素的一種空間分析工具，其將氣象數(shù)據(jù)點(diǎn)轉(zhuǎn)換為光滑連續(xù)曲面，以避免空間分析中的異常值風(fēng)險(xiǎn)。當(dāng)前華北地區(qū)城市化日益加深、農(nóng)業(yè)生產(chǎn)規(guī)模不斷擴(kuò)張，加之氣候變化與環(huán)境污染等深刻影響著該地自然生境[6-8]。而對該地自然資源可持續(xù)利用和低影響開發(fā)，需獲取準(zhǔn)確而精細(xì)的分布式降水量信息。本文擬采用半方差函數(shù)和Anusplin插值法定量區(qū)域降水量模式特征和空間規(guī)律。

1 數(shù)據(jù)和方法

1.1 華北地區(qū)簡述

文中華北地區(qū)特指北京、天津、河北、山西和內(nèi)蒙古自治區(qū)局部，經(jīng)緯度范圍介于35°35′N～42°40′N，113°20′E～119°50′E，總面積約60萬km2。該地位于中緯度太平洋西岸向亞洲內(nèi)陸延伸區(qū)，海拔范圍在0～3058m之間，總體地勢自西向東、自北向南傾斜，形成土石山、丘陵、平原河谷、草原沙漠等地貌。屬黃河、海河等流域，水文特征表現(xiàn)為流速平緩、冬季冰封、長期枯水。受海陸位置與季風(fēng)環(huán)流影響，該形成溫帶季風(fēng)氣向溫帶大陸性氣候過度區(qū)，物候分明、冬寒夏季熱、雨熱同季，年平均氣溫9～14℃，年降水量在400～800mm，無霜期140～220d，積溫1000～3500℃，干旱日數(shù)100～250d。區(qū)域典型氣象災(zāi)害包括干旱、冰雹、雷電、和霧霾等。

1.2 數(shù)據(jù)源和預(yù)處理

本研究選取均勻覆蓋華北地區(qū)的328個(gè)氣象站觀測資料，降水量觀測時(shí)間從2001年到2020年，如圖1所示。為保證數(shù)據(jù)質(zhì)量和完整性，首先對全部站點(diǎn)數(shù)據(jù)進(jìn)行識別缺漏、異常值刪除和插補(bǔ)等預(yù)處理，然后根據(jù)每一站點(diǎn)日降水量累積成年數(shù)據(jù)，并得到該站點(diǎn)多年平均降水量。

圖1 文中華北地區(qū)位置和雨量站分布

1.3 地統(tǒng)計(jì)學(xué)半方差函數(shù)方法

地統(tǒng)計(jì)學(xué)原本是針對地質(zhì)自然地理要素的一種空間統(tǒng)計(jì)學(xué)方法，其核心是半方差函數(shù)，主要作用是量化地理變量在不同空間位置上的相關(guān)性，即空間自相關(guān)。通常降水量在區(qū)域尺度上空間自相關(guān)程度可用半方差函數(shù)γ(x)表示，其計(jì)算公式如下：

(1)

式中，N—雨量站總數(shù)；h—空間距離，km；Z(xi)—在位置xi處降水量值，mm。

通常規(guī)定當(dāng)h=0時(shí)，此時(shí)的γ(x)值命名為塊金值值(Nugget)。可預(yù)見γ(x)隨著空間距離增加而逐漸增大，當(dāng)γ(x)趨于穩(wěn)定時(shí)的值命名為偏基臺值(Partial Sill)，基臺值(Sill)為Nugget與Partial Sill的加和，而此時(shí)的h命名為變程(A)，降水量在A范圍內(nèi)存在相關(guān)性，超出該范圍則相關(guān)性為0。塊金系數(shù)Nugget/Sill為降水量空間自相關(guān)性的定量參數(shù)，當(dāng)Nugget/Sill<0.25時(shí)，說明其空間自相關(guān)性強(qiáng)烈；當(dāng)其介于[0.25～0.75]時(shí)，說明其空間自相關(guān)性為中等；當(dāng)Nugget/Sill>0.75時(shí)，說明其空間自相關(guān)性微弱[4-6]。

1.4 Anusplin插值法

Anusplin插值法是一種廣泛應(yīng)用于地理學(xué)和環(huán)境科學(xué)領(lǐng)域的插值方法，可處理不規(guī)則網(wǎng)格數(shù)據(jù)。本研究利用Anusplin插值法得到華北地區(qū)的降水量分布圖，并將其與Ordinary Kriging插值法(OK)進(jìn)行對比。為確保插值空間分辨率和精度，對原始數(shù)據(jù)進(jìn)行l(wèi)og對數(shù)變換處理，在生成插值結(jié)果后再進(jìn)行反歸一化處理。設(shè)光滑函數(shù)為f、系數(shù)為b；xi為T維變量；則空間位置i處降水量值zi表示如下：

zi=f(xi)+bTyi+ei，(i=1，2……，N)

(2)

式中，ei—隨機(jī)誤差[6-8]。

2 結(jié)果與分析

2.1 雨量站降水量數(shù)據(jù)統(tǒng)計(jì)特征

利用Excel 2016軟件對研究區(qū)328個(gè)氣象站點(diǎn)降水資料進(jìn)行描述性統(tǒng)計(jì)，其結(jié)果見表1。可見研究區(qū)降水量最高值為705.9mm，出現(xiàn)在冀東北興隆縣；最低值為二連浩特站的175.5mm；最大、最小值相差530.4mm。統(tǒng)計(jì)得到全區(qū)多年平均降水量為420.8mm，按照全國降水量氣候分區(qū)屬半干旱區(qū)。觀測到其降水量中值與平均值相差不大，為437.1mm，離差系數(shù)達(dá)22.24%，屬中等程度變異性。利用Kolmogorov-Smirnov檢驗(yàn)表明站點(diǎn)序列降水量的sig值為0.03，未通過5%閾值檢驗(yàn)水平，不符合正態(tài)假設(shè)。對此采用log變換使其符合正態(tài)分布形式。

表1 華東地區(qū)雨量站降水量和降雨侵蝕力數(shù)據(jù)統(tǒng)計(jì)特征(n=328)

2.2 華北地區(qū)降水量空間自相關(guān)特征

半方差函數(shù)可直觀描述華北地區(qū)降水要素在不同空間位置之間變化。為擬合降水量空間變異結(jié)構(gòu)，利用GS+9.0軟件中的高斯、線性、球狀和指數(shù)模型對歸一化后的降水量數(shù)據(jù)進(jìn)行空間擬合，其結(jié)果見表2。通過比較4種空間模型的決定系數(shù)R2和殘差RSS，可知當(dāng)采用高斯模型時(shí)，擬合得到的R2最高，達(dá)0.68，相應(yīng)地殘差最小，僅為1.020E-02。在此基礎(chǔ)上利用ArcGIS 10.8軟件繪制了華北地區(qū)降水量半方差結(jié)構(gòu)，如圖2所示?？梢娖渥兂碳s為50km，塊金系數(shù)達(dá)53.62%，說明研究區(qū)降水量在50km范圍內(nèi)存在中等程度空間自相關(guān)性。由于該塊金系數(shù)介于25%～75%之間，說明區(qū)域降水量空間變異性由結(jié)構(gòu)性因素(氣候循環(huán)系統(tǒng)、海陸位置等)控制，即而隨機(jī)因素(如土地利用、人為活動)的影響相對微弱。

表2 降水要素的半方差結(jié)構(gòu)模型及參數(shù)

圖2 華北地區(qū)降水量半方差函數(shù)圖

2.3 華北地區(qū)降水量空間分布特征

運(yùn)用R語言中Anusplin程序包設(shè)計(jì)降水量空間插值的程序，得到華北地區(qū)1km空間分辨率的降水量柵格面，如圖3所示?？梢姡珹nusplin插值法生成的降水量范圍介于156～758mm，空間平均值為445.6mm，離差系數(shù)達(dá)26.78%，這與站點(diǎn)統(tǒng)計(jì)特征一致，同時(shí)反映了柵格面數(shù)據(jù)可靠性。為比較該插值方法與經(jīng)典OK插值法之間的差異性，兩種插值效果圖如圖3所示。可見，二者得到的區(qū)域降水量宏觀分布特征基本一致，即均呈現(xiàn)由西北向東南遞增格局，其中河北大部、晉南地區(qū)降水量明顯高于內(nèi)蒙地區(qū)，其中蒙東地區(qū)降水量略高于蒙西地區(qū)。需指出的是，Anusplin插值的降水量柵格面更具自然平滑特征(圖5a)，規(guī)避了OK方法產(chǎn)生的牛眼和機(jī)械性漸變(圖3b)，并且其清晰識別了特殊下墊面(例如城鎮(zhèn)等)條件下出現(xiàn)異常降水量分布信息，然而這一細(xì)節(jié)特征被OK算法概化。綜合而言，Anusplin插值生成的降水量柵格面更符合區(qū)域?qū)嶋H，且華北地區(qū)的降水量具有較強(qiáng)的空間變異，表現(xiàn)為大塊狀的區(qū)域差異[9]。

為進(jìn)一步比較不同插值法生成的降水量柵格精度差異，以站點(diǎn)觀測降水量為真值，計(jì)算了不同柵格的相對誤差。Anusplin方法的插值精度R2達(dá)0.62，比OK方法高出37.78%，而MAE和RMSE分別達(dá)58.9、89.4mm，相對于OK模型分別降低了15.49%、12.78%，見表3。

表3 不同插值方法的精度表現(xiàn)

3 結(jié)語

本文利用地統(tǒng)計(jì)學(xué)半方差函數(shù)和Anusplin插值法探究了華北地區(qū)的降水量空間分布特征。結(jié)果表明，華北地區(qū)降水量空間擬合模型以高斯模型最佳，局域降水量在50km范圍內(nèi)存在中等程度空間相關(guān)性；該地多年平均降水量具有強(qiáng)烈空間異質(zhì)性，總體呈現(xiàn)自西北向東南遞增格局，表現(xiàn)為塊狀區(qū)域性差異；相對于傳統(tǒng)OK模型，Anusplin插值法更適宜創(chuàng)建該地降水量柵格面，體現(xiàn)在空間平滑和識別微域降水量特征方面；生成的降水量柵格面可為該地水資源管理和地質(zhì)災(zāi)害預(yù)警提供科學(xué)數(shù)據(jù)源。后續(xù)研究應(yīng)關(guān)注半方差函數(shù)中相關(guān)參數(shù)對降水量空間結(jié)構(gòu)擬合精度的影響，例如二階高斯模型可容忍數(shù)據(jù)噪聲，并且在Anusplin插值法中融入空間自相關(guān)信息，進(jìn)而改善插值精度。