基于kriging方法的無(wú)定河流域年降水插值及空間場(chǎng)變異分析

基于kriging方法的無(wú)定河流域年降水插值及空間場(chǎng)變異分析

地質(zhì)統(tǒng)計(jì)是基于區(qū)域變量理論和區(qū)域變量理論的基礎(chǔ)，以突變函數(shù)為基本工具，研究和展示空間中的各種自然現(xiàn)象，并表現(xiàn)出一定的結(jié)構(gòu)性和隨機(jī)性。發(fā)展到現(xiàn)今已不單單應(yīng)用于地質(zhì)找礦和礦產(chǎn)估量中,它的使用范圍由逐步涉及環(huán)境科學(xué)、水科學(xué)、海洋學(xué)、農(nóng)田水利、土壤學(xué)、森林及漁業(yè)等領(lǐng)域,發(fā)展成為用來(lái)研究自然界具有隨機(jī)性和規(guī)律性雙重特征變量的具有普遍意義的科學(xué)方法。在水科學(xué)中,克立格方法已應(yīng)用于土壤水、徑流模數(shù)的估值,并且在布設(shè)水文與氣象站點(diǎn)時(shí),克立格也提供了強(qiáng)有力的理論依據(jù)。降雨量是具有隨機(jī)性與規(guī)律性雙重特征的區(qū)域化變量,本文以無(wú)定河流域1990年至1997年的年降雨量為對(duì)象,運(yùn)用克立格方法分析了該流域降雨量的空間變異特征。1kriger方法關(guān)于空間和空間的設(shè)計(jì)。據(jù)經(jīng)驗(yàn)克立格方法(Kriging)是由南非的地質(zhì)學(xué)家、采礦工程師克立格(D.G.Krige)于1951年提出,后來(lái)法國(guó)學(xué)者馬特隆(G.Matheon)于1962年提出“區(qū)域化變量”的概念后,得到不斷發(fā)展和完善。它是一種特定的滑動(dòng)加權(quán)平均法。這種方法假定采樣點(diǎn)間的距離和方向可反映一定的空間關(guān)聯(lián),并利用它們來(lái)解釋空間的變異性。它主要是利用一定的數(shù)學(xué)函數(shù)對(duì)特定點(diǎn)或是給定搜索半徑內(nèi)的所有點(diǎn)進(jìn)行擬合來(lái)估計(jì)每個(gè)點(diǎn)的值,對(duì)處理含有距離和方向上的偏差的數(shù)據(jù)尤為適用。要對(duì)Kriging方法中的區(qū)域化變量進(jìn)行估值,首先必須確定影響因素的范圍,即變程。而變程的確定,必須通過(guò)半變異函數(shù)模型的計(jì)算得出,這樣才能排除隨機(jī)因素的影響。同時(shí),對(duì)于Kriging方法中的協(xié)方差矩陣的計(jì)算,也需要預(yù)先知道半變異函數(shù)值,所以半變異函數(shù)是Kriging方法的基礎(chǔ)之基礎(chǔ),重中之重?？梢哉f(shuō)所有的Kriging計(jì)算都是基于半變異函數(shù)的分析而得出的。1.1空間協(xié)方差示克立格法的分析工具就是半變異函數(shù)。半變異函數(shù)r(h)來(lái)自于英文單詞(Semivariogram),如果用P表示空間內(nèi)一點(diǎn),用Z(P)表示空間內(nèi)這點(diǎn)的函數(shù),那么在所研究范圍內(nèi)的空間協(xié)方差用r(h)表示,其數(shù)學(xué)表達(dá)式為:r(h)=12Ν(h)Σ(Ζ(Xi)-Ζ(Xi+h))2(1)r(h)=12N(h)Σ(Z(Xi)?Z(Xi+h))2(1)式中,N(h)為相距h的數(shù)據(jù)對(duì)數(shù),h是距離矢量。從上面的定義可知,r(h)僅僅是任意兩點(diǎn)之間距離h的函數(shù)。根據(jù)r(h)和h點(diǎn)繪成圖(稱為試驗(yàn)半變異函數(shù)圖),用數(shù)學(xué)公式擬合,就可得到理想的半變異函數(shù)模型。1.2點(diǎn)krige估值的計(jì)算對(duì)任一變量在點(diǎn)x處的估計(jì)值Z*x,可以通過(guò)該點(diǎn)影響范圍內(nèi)n個(gè)有效觀測(cè)值Z(Xi)的線性組合得到,即:Ζ*x=nΣi=1λiΖ(Xi)i=1,2??,n(2)Z?x=Σi=1nλiZ(Xi)i=1,2??,n(2)式中,λi是賦予觀測(cè)值Z(Xi)的權(quán)重,表示各個(gè)觀測(cè)值Z(Xi)對(duì)估計(jì)值Z*x的貢獻(xiàn)。在進(jìn)行點(diǎn)Kriging估值計(jì)算時(shí),信息樣品就是點(diǎn)樣品值。這時(shí)所有權(quán)系數(shù)的計(jì)算十分簡(jiǎn)單,只需計(jì)算距離,克立格矩陣[K]只依賴于樣品X,第二元矩陣[Ma]為信息樣X(jué)i與待估點(diǎn)Xo之間的協(xié)方差。本文中所處理的其實(shí)就是點(diǎn)Kriging估值。在水科學(xué)中,克立格方法已經(jīng)涉及到地下水、地表水、土壤水等很多的領(lǐng)域。它不但有自己一套比較系統(tǒng)的數(shù)學(xué)、統(tǒng)計(jì)分析工具,而且對(duì)于結(jié)論也有比較好的驗(yàn)證方式。對(duì)于降雨量而言,一方面它受到經(jīng)緯度、高程等規(guī)律性因素的影響,同時(shí)又受到許多復(fù)雜的隨機(jī)性因素的影響,剛好滿足克立格的估值條件和要求。因此,可以運(yùn)用克立格方法對(duì)其進(jìn)行插值分析和計(jì)算。2降雨插值的確定克立格插值是基于正態(tài)分布進(jìn)行的。而在自然界中,很多不能取負(fù)的變量與對(duì)數(shù)正態(tài)分布有密切的關(guān)系。比如風(fēng)速、人口密度、生物個(gè)體大小、大氣氣溶膠粒徑、土壤和巖石中大多數(shù)微量元素等的頻率分布,都屬于或者接近對(duì)數(shù)正態(tài)分布,降雨量也是其中的一例。因此,在進(jìn)行降雨量的插值時(shí),首先要對(duì)數(shù)值進(jìn)行對(duì)數(shù)正態(tài)變換。本文以黃河中游無(wú)定河流域?yàn)檠芯繉?duì)象,選取其中67個(gè)雨量站1990～1997年的降雨量為數(shù)據(jù)源。同時(shí)克立格方程的矩陣形式或矩陣中的具體符號(hào)含義也有所變化:對(duì)于點(diǎn)的插值變?yōu)?lnΖ*v=C+nΣi=1λilnΖ(Xi)(5)lnZ?v=C+Σi=1nλilnZ(Xi)(5)式中,C、λi為待定系數(shù),Xi為定義于信息支撐Vi(i=1,2,…,n)的n個(gè)信息樣品觀測(cè)值。3降水的單次測(cè)量與半變異函數(shù)降雨、徑流、侵蝕產(chǎn)沙之間存在著緊密的關(guān)系,若能找出降雨與徑流之間、徑流與產(chǎn)沙之間的數(shù)量關(guān)系,就可以根據(jù)降雨的變化趨勢(shì)預(yù)測(cè)徑流、侵蝕產(chǎn)沙的變化趨勢(shì)。這對(duì)研究區(qū)域的水土流失、環(huán)境變遷,尤其在黃土高原地區(qū)是很有現(xiàn)實(shí)意義的。無(wú)定河地處黃土高原區(qū),是黃河中游最大的多沙支流。本文應(yīng)用地質(zhì)統(tǒng)計(jì)方法對(duì)無(wú)定河流域的降水進(jìn)行了時(shí)空變異性分析。無(wú)定河發(fā)源于陜西省靖邊、定邊、吳旗三縣交界的白于山,是黃河中游河口-龍門(mén)區(qū)間最大的支流,由西向東經(jīng)內(nèi)蒙古伊克昭盟和陜西省榆林、延安地區(qū),于清澗縣河口村匯入黃河。無(wú)定河干流全長(zhǎng)491.1km,全河比降為1.97‰,流域內(nèi)總面積為30261km2,其中水土流失面積達(dá)23137km2,占總面積的76.5%,侵蝕模數(shù)達(dá)8000t/km2.a。流域?qū)贉貛Т箨懶愿珊蛋敫珊禋夂蝾愋?。按地形地貌及水土流失特點(diǎn),可分為風(fēng)沙區(qū)、河源梁澗區(qū)及丘陵溝壑區(qū)。該區(qū)降雨量具有量少、集中、強(qiáng)度大的特點(diǎn)。據(jù)1957～1989年白家川站的實(shí)測(cè)資料統(tǒng)計(jì)可知,無(wú)定河流域多年平均年降雨量409mm,降水年內(nèi)分配的特點(diǎn)是:降雨以汛期(6～9月)為主,主要集中在7、8月份,而且這些雨量多由若干次高強(qiáng)度的暴雨、大暴雨構(gòu)成。汛期降雨量約占全年雨量的75%,極端最大月降雨量可達(dá)到年降雨量的70%以上。此外,年降雨量的變率也很大,最高年變率達(dá)6～71。流域內(nèi)有雨量站90多個(gè),其中近一半是1975年以后設(shè)立的。主要雨量站的位置分布如圖1所示。表1中列舉了韓家峁、靖邊、橫山、殿市、丁家溝、李家河等6個(gè)代表性雨量站1990～1997的降雨量數(shù)據(jù)。該流域年降雨量經(jīng)對(duì)數(shù)化后的基本統(tǒng)計(jì)資料如表2。變異系數(shù)最小0.28,最大為0.42,說(shuō)明了降雨量在時(shí)間上分布的不均勻性。對(duì)樣點(diǎn)按照公式(1)作出實(shí)驗(yàn)半變異函數(shù)圖。根據(jù)圖分析發(fā)現(xiàn),試驗(yàn)半變異函數(shù)符合球狀模型。對(duì)于較易識(shí)別的實(shí)驗(yàn)半變異函數(shù)圖可用直接法對(duì)球狀模型擬合,求出待定參數(shù)。對(duì)于用直接法難以得出模型參數(shù)值的情況,可以運(yùn)用最小二乘法的原則,采取加權(quán)多項(xiàng)式回歸法對(duì)變異曲線進(jìn)行球狀模型的最優(yōu)擬合(圖2)。球狀模型的一般形式:r(h)={C0+C[32ha-12h3a3]h≤aC0+Ch>a(6)確定理論變異函數(shù)模型的最終目的在于,定量地揭示降雨量區(qū)域化變量的空間變異規(guī)律,并將其結(jié)果直接用于降雨量的空間最優(yōu)插值。因此,如果理論變異函數(shù)模型確定得好且較為符合實(shí)際,就應(yīng)該最終體現(xiàn)在根據(jù)周圍信息點(diǎn)對(duì)實(shí)測(cè)點(diǎn)的空間估計(jì)值與實(shí)測(cè)值之差的平方平均最小,我們把它作為擬合理論變異函數(shù)模型“最優(yōu)性”檢驗(yàn)的標(biāo)準(zhǔn)。計(jì)算半變異函數(shù)的幾點(diǎn)認(rèn)識(shí)和說(shuō)明:在實(shí)際計(jì)算半變異函數(shù)時(shí),由于數(shù)據(jù)的間距不規(guī)則,通常的做法是按區(qū)間[k-ε(h),k+ε(h)]將所有的距離組合成距離組,用距離在[k±ε(h)]內(nèi)的所有數(shù)據(jù)對(duì)來(lái)計(jì)算r(kh)(其中ε(h)稱為距離誤差限),一般在計(jì)算時(shí)選擇距離誤差限ε(h)≤2/h。本文中定為2/h。對(duì)于一個(gè)特定的數(shù)據(jù)樣點(diǎn)集,間距小,則意味著實(shí)驗(yàn)半變異函數(shù)曲線不平滑,落在每個(gè)距離段內(nèi)的點(diǎn)對(duì)數(shù)相應(yīng)也少,不利于分析。而間距的增加,意味著塊金效應(yīng)的增大,這說(shuō)明在某種程度上半變異函數(shù)受到人為因素的影響,所以在計(jì)算和擬合時(shí),分隔距離的選取就顯得尤為重要。本文通過(guò)設(shè)置滑動(dòng)桿的方法,在4km到20km這個(gè)距離范圍搜索最佳分隔距離。同時(shí),因?yàn)樽儺惡瘮?shù)只有在最大間距的1/2內(nèi)才有意義,所以只采用間距小于或等于最大距離一半的點(diǎn)對(duì)參與實(shí)際計(jì)算。在所分析的系列中,所有年份都存在塊金效應(yīng)。塊金效應(yīng)反映了隨機(jī)變化,塊金值通常表示由實(shí)驗(yàn)誤差和小于實(shí)驗(yàn)取樣尺度引起的變異,較大的塊金方差值表明較小尺度上的某種過(guò)程不容忽視。若塊金效應(yīng)為零,也并不表示數(shù)據(jù)無(wú)測(cè)量誤差或隨機(jī)特性,只是說(shuō)明隨機(jī)效應(yīng)的尺度小于相關(guān)尺度——變程,設(shè)置的距離段太大而掩飾了或不足以表現(xiàn)出這種隨機(jī)效應(yīng)或測(cè)量誤差。因此,將零塊金效應(yīng)理解為塊金值很小更為合理。在所研究的時(shí)間序列中,變程從20.17km到68.65km不等。變程反映的是空間點(diǎn)的自相關(guān)性,變程大則自相關(guān)的尺度廣。在我們所研究的時(shí)段中,1994年的變程最大,1991年的變程最小,說(shuō)明前者的降雨量變化比較緩和,而后者的降雨量空間變化相對(duì)較烈。除了用σ*e來(lái)檢驗(yàn)擬合模型的好壞外,通常我們也可運(yùn)用對(duì)比檢驗(yàn)的方法來(lái)分析所選方法的適用性。4降雨空間變化的空間異質(zhì)性分析通過(guò)參數(shù)的求得,運(yùn)用選定的模型可以對(duì)區(qū)域內(nèi)的點(diǎn)或塊段進(jìn)行插值或是成圖。對(duì)于降雨量而言,我們根據(jù)選定的克立格方法和擬合的模型,做出各個(gè)年份的降雨量等值線圖,并基于此作出降雨量三維空間分布圖(圖3,圖版4)。這樣,對(duì)于降雨量的空間分布和時(shí)間變異特征,我們就能比較容易地識(shí)別出來(lái)。同時(shí)若能結(jié)合地形和高程圖,運(yùn)用數(shù)學(xué)統(tǒng)計(jì)方法,我們還能找出影響降雨量的主導(dǎo)因素。特別是三維空間分布圖,我們能夠直觀地獲得降雨量的空間變異情況,同時(shí)對(duì)比不同年份的圖幅,我們還能獲得時(shí)間上的變化情況。選擇分布較為集中的50個(gè)站點(diǎn)作為空間插值成圖區(qū),即圖1中的虛線區(qū)域。對(duì)于選定的區(qū)域而言,其東北、西北部降雨量變化幅度較小,同時(shí)降雨量也比較少,在圖上表現(xiàn)出比較平緩的趨勢(shì)(1992年除外),流域的中南部變化較大,在圖上表現(xiàn)為比較劇烈的起伏。而1997年的降雨量空間場(chǎng)變化和以往各年相比又有明顯的不同,首先該年的南部,尤其是中南部降雨量較少(與以往年份相比)