基于快速傅里葉變換的降水空間變異函數(shù)有效性驗(yàn)證

陳紫陽(yáng)，胡慶芳，雍 斌，李伶杰，王銀堂

(1.河海大學(xué)地球科學(xué)與工程學(xué)院，江蘇 南京 211100；2.南京水利科學(xué)研究院水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210029)

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

在描述降水空間變異性的眾多方法或指標(biāo)中，空間變異函數(shù)(或變方差函數(shù))應(yīng)用較為廣泛[10-12]?？臻g變異函數(shù)的計(jì)算方法有兩種：傳統(tǒng)方法(基于某種理論變異函數(shù)的參數(shù)擬合方法)和基于快速傅里葉變換(fast Fourier transformation，F(xiàn)FT)的計(jì)算方法[13-15]。前者作為傳統(tǒng)的計(jì)算方法，需要選擇特定模型，針對(duì)實(shí)測(cè)數(shù)據(jù)得到的經(jīng)驗(yàn)半方差函數(shù)進(jìn)行參數(shù)擬合，且一般只能得到特定方向的變異函數(shù)值；而后者則不需要進(jìn)行參數(shù)擬合，能夠直接提供二維矩陣形式的變異函數(shù)計(jì)算結(jié)果，可以反映任意方向上的降水空間變異性，這是相對(duì)傳統(tǒng)方法的突出優(yōu)點(diǎn)。

基于FFT的降水空間變異函數(shù)計(jì)算方法(以下簡(jiǎn)稱“FFT方法”)在降水空間變異分析和定量估計(jì)中已得到應(yīng)用，如Velasco-Forero等[16]采用FFT方法計(jì)算了西班牙Catalonia地區(qū)6場(chǎng)降雨的變異函數(shù)，并應(yīng)用于降水空間插值中；Cecinati[17]則將該方法應(yīng)用于英國(guó)雷達(dá)-雨量計(jì)降水聯(lián)合估計(jì)中。但總的來(lái)說(shuō)，關(guān)于該方法的研究仍然很不充分，本文以淮河中上游某區(qū)域?yàn)槔_(kāi)展基于FFT方法的有效性研究，并與傳統(tǒng)的空間變異函數(shù)計(jì)算方法進(jìn)行比較，在不同時(shí)間尺度上對(duì)FFT方法的有效性加以驗(yàn)證。

1 空間變異函數(shù)計(jì)算方法

1.1 傳統(tǒng)方法

空間變異函數(shù)是地質(zhì)統(tǒng)計(jì)學(xué)中用于衡量區(qū)域化變量在不同位置差異性的統(tǒng)計(jì)指標(biāo)[18-19]。對(duì)于區(qū)域化變量Z(xi)和某一空間分離距離h，若Z(xi)均值存在且不取決于空間位置xi，同時(shí)對(duì)于任意xi和h，Z(xi+h)-Z(xi)的方差存在，且不取決于xi，則Z(xi)的空間變異函數(shù)可表達(dá)為

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

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

(1)

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

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

(2)

式中：C(h)為協(xié)方差函數(shù)；ρ(h)為空間自相關(guān)系數(shù)；σ2為方差。

傳統(tǒng)的空間變異函數(shù)計(jì)算主要包括以下兩個(gè)步驟：

a. 根據(jù)分離距離為h的樣本數(shù)據(jù)對(duì)，采用下式計(jì)算經(jīng)驗(yàn)變異函數(shù)值：

(3)

式中Ph為分離距離為h的樣本數(shù)據(jù)對(duì)的數(shù)量。

b. 選擇滿足正定性要求的有效變異函數(shù)模型對(duì)經(jīng)驗(yàn)變異函數(shù)值進(jìn)行參數(shù)擬合，得到理論變異函數(shù)。常用的理論變異函數(shù)模型主要有球狀模型、指數(shù)模型和高斯模型等。以高斯模型為例，其對(duì)應(yīng)的變異函數(shù)表達(dá)式為

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

(4)

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

由于區(qū)域化變量在不同方向可能會(huì)具有不同的結(jié)構(gòu)性和變異性(即空間各向異性)，故為全面了解區(qū)域化變量的空間分布特性，需要對(duì)不同方向分別計(jì)算空間變異函數(shù)，同時(shí)還需要根據(jù)各向異性的具體性質(zhì)，對(duì)不同方向上的變異函數(shù)進(jìn)行擬合，這正是采用傳統(tǒng)方法較為麻煩之處。

1.2 FFT方法

與傳統(tǒng)方法不同，F(xiàn)FT方法不需要針對(duì)各個(gè)方向采用特定模型對(duì)經(jīng)驗(yàn)變異函數(shù)值進(jìn)行擬合，它是一種非參數(shù)計(jì)算方法。這種方法首先計(jì)算出區(qū)域化變量的協(xié)方差矩陣，然后用FFT將協(xié)方差矩陣轉(zhuǎn)換到頻域進(jìn)行處理，最后再通過(guò)快速傅里葉逆變換(inverse fast Fourier transform,IFFT)將協(xié)方差函數(shù)變換回時(shí)域，最后利用變異函數(shù)和協(xié)方差函數(shù)之間的轉(zhuǎn)換關(guān)系得到二維矩陣形式的空間變異函數(shù)。

FFT方法的具體步驟如下：

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

b. 采用下式計(jì)算分離距離h對(duì)應(yīng)的協(xié)方差函數(shù)值C(h)：

(5)

由于采樣點(diǎn)數(shù)量有限，因此并非研究區(qū)域內(nèi)的所有格網(wǎng)均是有數(shù)據(jù)的，這導(dǎo)致可能沒(méi)有樣本數(shù)據(jù)對(duì)用于計(jì)算某些h對(duì)應(yīng)的協(xié)方差函數(shù)值，因此在用FFT進(jìn)行變換之前要將協(xié)方差函數(shù)進(jìn)行初步插值，本文采用下式進(jìn)行計(jì)算：

(6)

其中

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

c. 采用下式所示的FFT將協(xié)方差矩陣從時(shí)域轉(zhuǎn)換到頻域：

(7)

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

(8)

d. 采用式(9)的IFFT將Smn轉(zhuǎn)換回時(shí)域，在協(xié)方差函數(shù)矩陣的基礎(chǔ)上，根據(jù)式(2)得到優(yōu)化后的變異函數(shù)矩陣。

(9)

2 研究區(qū)域與數(shù)據(jù)

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

研究區(qū)域內(nèi)共有雨量站127個(gè)，面積約28 900 km2，站點(diǎn)密度約為230 km2一個(gè)。研究數(shù)據(jù)為從《淮河流域水文年鑒》中收集的127個(gè)站點(diǎn)2010—2015年連續(xù)的年、月、日降水?dāng)?shù)據(jù),降水?dāng)?shù)據(jù)經(jīng)過(guò)了嚴(yán)格的質(zhì)量控制。

圖1 研究區(qū)域與雨量站點(diǎn)分布

3 結(jié)果分析與討論

根據(jù)研究區(qū)域內(nèi)2010—2015年的降水?dāng)?shù)據(jù)，采用FFT方法計(jì)算了年、月、日3種時(shí)間尺度上的降水空間變異函數(shù)，從東-西(E-W)、東南-西北(SE-NW)、南-北(S-N)、西南-東北(SW-NE)4個(gè)方向比較了FFT方法與傳統(tǒng)方法的差異。

3.1 年降水量空間變異函數(shù)

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

圖2 FFT方法計(jì)算的年降水量空間變異函數(shù)

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

(a) E-W

(b) SE-NW

(c) S-N

(d) SW-NE

3.2 月降水量空間變異函數(shù)

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

圖4 FFT方法計(jì)算的月降水量空間變異函數(shù)

圖5為兩種方法計(jì)算的月降水量(2012年9月)空間變異函數(shù)對(duì)比，其中傳統(tǒng)方法選取了擬合效果最好的高斯模型。總體而言，其規(guī)律與年降水量類似，傳統(tǒng)擬合方法計(jì)算結(jié)果與經(jīng)驗(yàn)變異函數(shù)計(jì)算值更為接近，但傳統(tǒng)方法和FFT方法所得結(jié)果在多數(shù)情況下具有相近的變化特征。

3.3 日降水量空間變異函數(shù)

圖6為FFT方法計(jì)算的日降水量(2015年6月25日)空間變異函數(shù)。由圖6可知，日降水量在S-N方向變異較慢，E-W方向變異快，但即便如此仍有一定差異。日降水量空間變異函數(shù)的塊金值約為 2 mm2，基臺(tái)值大致為50 mm2，S-N方向變程約為110 km，E-W方向變程約為80 km。

圖7給出了兩種方法計(jì)算的日降水量(2015年6月25日)空間變異函數(shù)值對(duì)比，其中傳統(tǒng)方法選取了擬合效果最好的指數(shù)模型。FFT方法所得日降水量空間變異函數(shù)在長(zhǎng)距離上的表現(xiàn)沒(méi)有傳統(tǒng)方法效果好，但在短距離上兩者所得結(jié)果較為接近。

4 結(jié) 語(yǔ)

本文在介紹空間變異函數(shù)計(jì)算方法的基礎(chǔ)上，以淮河中上游某區(qū)域?yàn)槔?，通過(guò)與傳統(tǒng)方法的比較，在年、月、日時(shí)間尺度上對(duì)FFT方法的計(jì)算結(jié)果進(jìn)行了驗(yàn)證。在年、月時(shí)間尺度上，F(xiàn)FT方法和傳統(tǒng)方法所得降水空間變異函數(shù)具有相近變化特征，而在日時(shí)間尺度上，兩者在短距離上所得結(jié)果也較為接近。

(a) E-W

(b) SE-NW

(c) S-N

(d) SW-NE

圖6 FFT方法計(jì)算的日降水量空間變異函數(shù)

(a) E-W

(b) SE-NW

(c) S-N

(d) SW-NE

本文的研究證實(shí)了FFT方法在細(xì)致解析降水空間結(jié)構(gòu)特征上的價(jià)值，下一步將運(yùn)用該方法開(kāi)展降水空間定量估計(jì)的研究，并比較不同降水空間變異函數(shù)計(jì)算方法對(duì)降水定量估計(jì)結(jié)果的影響。