基于kriging方法的無定河流域年降水插值及空間場變異分析

基于kriging方法的無定河流域年降水插值及空間場變異分析

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