基于快速傅里葉變換的降水空間變異函數有效性驗證

2020-05-26 03:56:30陳紫陽胡慶芳李伶杰王銀堂

水資源保護 2020年3期

陳紫陽，胡慶芳，雍斌，李伶杰，王銀堂

(1.河海大學地球科學與工程學院，江蘇南京 211100；2.南京水利科學研究院水文水資源與水利工程科學國家重點實驗室，江蘇南京 210029)

降水是最基本的氣象水文要素之一，受大氣運動狀態和地理環境因素的共同影響，降水具有復雜的多尺度空間變異性。深入解析降水空間變異性是科學認識降水時空演化規律并定量估計降水空間分布的重要基礎[1-3]。降水空間變異性除通過等值線等方式直接展現外，更主要的是采用不均勻指數、空間相關系數、信息熵、經驗正交函數分解(empirical or thogonal function, EOF)等統計指標或方法加以分析[4]。如李邦東等[5]采用不均勻指數(變差系數)分析了我國東北地區1961—2010年降水事件，發現在1993年后該地區降水空間不均勻性增加，在一定程度上增加了區域旱澇風險；胡慶芳等[6]采用Moran指數等6種全局或局部空間相關性指標分析了贛江流域降水空間變異性，得出年、月、日降水量整體上呈顯著時變空間正相關性，但在局部上具有非平穩性和奇異性的結論；張繼國等[7]基于信息熵傳輸模型，解析了淮河流域年降水在東西方向相關性強、南北方向弱的空間結構特征；徐利崗等[8]基于EOF發現我國北方荒漠區降水總體表現為“相間復雜”和“東西相反”兩種空間結構形態；方國華等[9]采用降水質心時間為評價指標，研究了淮河流域降水結構的時空演變規律。

在描述降水空間變異性的眾多方法或指標中，空間變異函數(或變方差函數)應用較為廣泛[10-12]。空間變異函數的計算方法有兩種：傳統方法(基于某種理論變異函數的參數擬合方法)和基于快速傅里葉變換(fast Fourier transformation，FFT)的計算方法[13-15]。前者作為傳統的計算方法，需要選擇特定模型，針對實測數據得到的經驗半方差函數進行參數擬合，且一般只能得到特定方向的變異函數值；而后者則不需要進行參數擬合，能夠直接提供二維矩陣形式的變異函數計算結果，可以反映任意方向上的降水空間變異性，這是相對傳統方法的突出優點。

基于FFT的降水空間變異函數計算方法(以下簡稱“FFT方法”)在降水空間變異分析和定量估計中已得到應用，如Velasco-Forero等[16]采用FFT方法計算了西班牙Catalonia地區6場降雨的變異函數，并應用于降水空間插值中；Cecinati[17]則將該方法應用于英國雷達-雨量計降水聯合估計中。但總的來說，關于該方法的研究仍然很不充分，本文以淮河中上游某區域為例，開展基于FFT方法的有效性研究，并與傳統的空間變異函數計算方法進行比較，在不同時間尺度上對FFT方法的有效性加以驗證。

1 空間變異函數計算方法

1.1 傳統方法

空間變異函數是地質統計學中用于衡量區域化變量在不同位置差異性的統計指標[18-19]。對于區域化變量Z(xi)和某一空間分離距離h，若Z(xi)均值存在且不取決于空間位置xi，同時對于任意xi和h，Z(xi+h)-Z(xi)的方差存在，且不取決于xi，則Z(xi)的空間變異函數可表達為

γ(h)=V[Z(xi+h)-Z(xi)]/2=

E{[Z(xi+h)-Z(xi)]2}/2

(1)

式中：V為方差；E為均值。顯然，若Z(xi)呈空間正相關，則|h|越大，變異函數值γ(h)越大；反之，γ(h)越小。在Z(xi)滿足二階平穩條件的情況下，變異函數與協方差函數以及空間相關系數存在以下轉換關系：

γ(h)=σ2-C(h)=σ2[1-ρ(h)]

(2)

式中：C(h)為協方差函數；ρ(h)為空間自相關系數；σ2為方差。

傳統的空間變異函數計算主要包括以下兩個步驟：

a. 根據分離距離為h的樣本數據對，采用下式計算經驗變異函數值：

(3)

式中Ph為分離距離為h的樣本數據對的數量。

b. 選擇滿足正定性要求的有效變異函數模型對經驗變異函數值進行參數擬合，得到理論變異函數。常用的理論變異函數模型主要有球狀模型、指數模型和高斯模型等。以高斯模型為例，其對應的變異函數表達式為

γ(h)=C0+C1[1-e-(h/a)2]

(4)

式中：a為變程，反映區域化變量的相關尺度；C0為塊金值，反映極短分離距離之間變異函數的不連續性；C1為變異系數，反映區域化變量變化速率的大小。a、C0、C1為高斯模型的參數。

由于區域化變量在不同方向可能會具有不同的結構性和變異性(即空間各向異性)，故為全面了解區域化變量的空間分布特性，需要對不同方向分別計算空間變異函數，同時還需要根據各向異性的具體性質，對不同方向上的變異函數進行擬合，這正是采用傳統方法較為麻煩之處。

1.2 FFT方法

與傳統方法不同，FFT方法不需要針對各個方向采用特定模型對經驗變異函數值進行擬合，它是一種非參數計算方法。這種方法首先計算出區域化變量的協方差矩陣，然后用FFT將協方差矩陣轉換到頻域進行處理，最后再通過快速傅里葉逆變換(inverse fast Fourier transform,IFFT)將協方差函數變換回時域，最后利用變異函數和協方差函數之間的轉換關系得到二維矩陣形式的空間變異函數。

FFT方法的具體步驟如下：

a. 首先根據樣本數據的空間分布特征，將研究區域劃分成一系列矩形格網，格網橫向和縱向的分辨率分別為x、y。各格網中區域化變量的數值為分布在格網中的采樣點的均值，分離距離h=(mx,ny)(m、n為格網相對平移行列數，m、n可為負數)。為方便計算，格網行列數m、n取2k+1(k>2)。

b. 采用下式計算分離距離h對應的協方差函數值C(h)：

(5)

由于采樣點數量有限，因此并非研究區域內的所有格網均是有數據的，這導致可能沒有樣本數據對用于計算某些h對應的協方差函數值，因此在用FFT進行變換之前要將協方差函數進行初步插值，本文采用下式進行計算：

(6)

其中

式中：Cmn為待插值點的協方差值；m、n為所在行列數；P為參與插值的點的數量；λi為插值權重；Ci為插值點的協方差值；pi為計算相應格網協方差函數值所采用的樣本數據對數量；di為插值格網與待插值格網間的距離。由于區域化變量常具有各向異性，因此在對協方差函數插值時需要考慮方向。首先連接中心點與待插值格網，然后在其連線左右一定容差角度范圍內，篩選出與待插值格網相距較近的插值點，最后根據這些距離待插值點較近點的協方差函數值即可擬合出待插值格網的協方差函數值。

c. 采用下式所示的FFT將協方差矩陣從時域轉換到頻域：

(7)

式中：M、N為協方差矩陣的行數和列數；Smn為頻域中第m行、第n列的頻譜值。為保證正定性條件，采用式(7)獲得頻譜圖后需要用式(8)所示的移動窗口平滑算法對各個頻譜值進行優化，逐步加大窗口寬度直到每個頻譜值為非負。如果窗口太大而頻譜值仍然為負，則可以直接將所有頻譜值加上所有值中最小值的絕對值(最小值太小會使整個頻譜基數變大從而使最終的變異函數變化平緩)或者直接將此結果設置為0來進行校正。

(8)

d. 采用式(9)的IFFT將Smn轉換回時域，在協方差函數矩陣的基礎上，根據式(2)得到優化后的變異函數矩陣。

(9)

2 研究區域與數據

選擇淮河流域中上游一長寬均為170 km的方形研究區域(圖1)來驗證FFT方法計算結果的合理性。研究區域位于安徽省阜陽附近(東經115°10′～116°51′，北緯32°8′～33°41′)，地跨淮河干流兩岸，屬于我國南北過渡帶，高程介于-11～409 m，多年平均降水量為900～1 000 mm。

研究區域內共有雨量站127個，面積約28 900 km2，站點密度約為230 km2一個。研究數據為從《淮河流域水文年鑒》中收集的127個站點2010—2015年連續的年、月、日降水數據,降水數據經過了嚴格的質量控制。

圖1 研究區域與雨量站點分布

3 結果分析與討論

根據研究區域內2010—2015年的降水數據，采用FFT方法計算了年、月、日3種時間尺度上的降水空間變異函數，從東-西(E-W)、東南-西北(SE-NW)、南-北(S-N)、西南-東北(SW-NE)4個方向比較了FFT方法與傳統方法的差異。

3.1 年降水量空間變異函數

圖2為FFT方法計算得到的研究區域年降水量(2012年)的空間變異函數(為二維矩陣，圖中每個像元的大小為10 km×10 km，圖像大小為 330 km×330 km，以向E、N為正，向W、S為負，圖像中心處表示分離距離為0的變異函數值，其他各位置表示不同分離距離和方向上的變異函數值)。顯然，空間變異函數矩陣具有中心對稱性，研究區域年降水量在空間上呈現出不同程度、不同特點的帶狀各向異性。對于2012年年降水量，其在SE-NW方向空間變異較快，而在SW-NE方向的變異相對較慢。從圖2還可發現，變異函數的塊金值約為 1 800 mm2，基臺值大致為18 000 mm2，變程為80 km。在變程范圍內，年降水量是結構性和隨機性的統一體。隨著分離距離的不斷增加，年降水量的空間結構性逐步降低、隨機性不斷增強。

圖2 FFT方法計算的年降水量空間變異函數

圖3為兩種方法計算的年降水量(2012年)空間變異函數對比(為對比兩種方法的差異性，圖中也給出了采用經驗變異函數直接計算得到的散點值)，其中傳統方法采用了多種模型進行擬合，最終選取了擬合效果最好的球形模型。總體而言，傳統方法計算結果與經驗變異函數計算值更為接近，且傳統方法和FFT方法所得結果在多數情況下具有相近變化特征，如E-W方向兩種方法所得年降水量空間變異函數值就很接近。

(a) E-W

(b) SE-NW

(d) SW-NE

3.2 月降水量空間變異函數

圖4為FFT方法計算的月降水量(2012年9月)空間變異函數。由圖4可知，研究區域內月降水量的空間帶狀各向異性特征也很明顯。月降水量在E-W和S-N方向變化慢，在SE-NW和SW-NE方向變化慢。從圖4還可以看出月降水量空間變異函數的塊金值約為100 mm2，基臺值大致為1 800 mm2，S-N方向變程約為100 km，E-W方向變程約為 80 km。與年降水量相似，隨著分離距離的增加，月降水量的半方差函數值不斷增大，隨機性不斷增強；直至分離距離達到和超過變程后，月降水量的半方差函數值趨于方差值。