基于半方差函數和Anusplin方法華北地區降水量空間分布研究

王克磊，曲華斌，王 爽

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

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

1 數據和方法

1.1 華北地區簡述

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

1.2 數據源和預處理

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

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

1.3 地統計學半方差函數方法

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

(1)

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

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

1.4 Anusplin插值法

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

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

(2)

式中，ei—隨機誤差[6-8]。

2 結果與分析

2.1 雨量站降水量數據統計特征

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

表1 華東地區雨量站降水量和降雨侵蝕力數據統計特征(n=328)

2.2 華北地區降水量空間自相關特征

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

表2 降水要素的半方差結構模型及參數

圖2 華北地區降水量半方差函數圖

2.3 華北地區降水量空間分布特征

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

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

表3 不同插值方法的精度表現

3 結語

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