土壤屬性空間預測精度與不確定性分析——有限最大似然法和高程輔助變量的應用

1、摘 要土壤屬性空間變化的準確預測是土壤學及環(huán)境科學等相關(guān)學科的一個共同主題。在全球變化、資源與環(huán)境、生態(tài)多樣性、食品安全以及人口與耕地等一系列問題的相關(guān)研究過程中，對土壤屬性空間分布信息，不論在信息的數(shù)量方面還是在準確性方面，都提出了越來越高的要求。本文從準確性和不確定性兩個方面，對有限最大似然法（REML）和高程輔助變量在土壤屬性空間預測中的應用進行了探討和研究，主要內(nèi)容與結(jié)論如下：（1）以北京市大興研究樣區(qū)作為研究樣區(qū)，以土壤有機質(zhì)、土壤含水量、土壤速效鉀和土壤有效錳四個土壤屬性作為目標變量，通過選擇不同的樣本點數(shù)，對比分析REML法和傳統(tǒng)的矩量法（MoM）在計算變異函數(shù)準確性方面的表現(xiàn)。

2、結(jié)果表明：在樣本點數(shù)從150個逐漸減少到50個過程中，當樣本點數(shù)小于70個后，所有目標變量的預測精度都開始明顯的下降。在樣本點數(shù)比較少的情況下，REML法估計的變異函數(shù)比MoM法估計的變異函數(shù)更準確一些。 （2）以土壤有機質(zhì)和土壤有效錳作為目標變量，通過選擇不同的空間尺度，分析空間尺度對土壤屬性變異函數(shù)準確性的影響。結(jié)果表明：空間尺度對于認識土壤目標變量的空間變異特征具有明顯的影響；REML法在提高空間預測精度方面效果是否明顯與研究區(qū)的空間尺度也緊緊相關(guān)，本文中，REML法估計的變異函數(shù)在20km和30km兩個尺度范圍內(nèi)比MoM法變異函數(shù)的預測精度高。（3）以平谷區(qū)為研究樣區(qū)，以高程作為輔助變

3、量，以土壤有機質(zhì)、土壤速效鉀、土壤有效磷、土壤有效鐵等為目標變量，通過利用普通克里格法（OK）、協(xié)克里格法（CK）和回歸克里格法（RK）三種方法，對比分析高程數(shù)據(jù)是否可以用來提高上述土壤屬性的空間預測精度。結(jié)果表明：高程作為輔助變量，可以用來提高土壤屬性的空間預測精度，但這種方法并不適合所有土壤屬性，本文中土壤有效銅、土壤有效鐵和土壤有效錳三種微量元素的預測精度沒能夠被提高；利用高程對土壤屬性進行空間預測時，在選用最適宜的方法之前，應該對土壤變量的空間結(jié)構(gòu)、土壤屬性的全局趨勢、土壤屬性與高程之間的線性相關(guān)關(guān)系、結(jié)構(gòu)相關(guān)關(guān)系等進行仔細分析。（4）利用高程、坡度和地形濕度指數(shù)等地形指標組合成不同的

4、外部趨勢模型，通過RK法和有限最大似然法為基礎(chǔ)的線性無偏最優(yōu)估計法（BLUP-REML）對平谷區(qū)土壤有機質(zhì)的空間分布進行預測。結(jié)果表明：地形指標高程、坡度和地形濕度指數(shù)等可以作為輔助變量用來提高土壤有機質(zhì)的空間預測精度，尤其是通過不同地形指標之間的組合，能進一步的提高預測精度；在樣本點數(shù)較少的前提下，在外部趨勢模型對土壤有機質(zhì)空間變化解釋程度偏低的情況下， BLUP-REML法在預測的穩(wěn)定性方面比RK法表現(xiàn)得更好。（5）將高程輔助數(shù)據(jù)與序貫高斯協(xié)模擬方法相結(jié)合，對密云樣區(qū)和平谷樣區(qū)的土壤有機質(zhì)的空間分布進行模擬，通過對模擬結(jié)果的準確性、局部不確定性和空間不確定性進行評價。結(jié)果表明，在條件高斯模

5、擬的計算環(huán)境下，可以利用高程作為輔助變量來提高土壤有機質(zhì)空間預測精度和空間不確定性模擬的準確性，但沒能提高局部不確定性模擬的準確性。關(guān)鍵詞：土壤屬性，空間預測精度，不確定性，地統(tǒng)計方法，高程，有限最大殘差法AbstractSpatial prediction of soil properties has become a common topic in soil science research. There are increasing demands for the information of spatial distribution of soil characteristics, b

6、oth in the amount of information and in terms of accuracy, for the study of the issues related to changes in the global resources and the environment, biodiversity, food security and population and arable land. In this paper, we discussed the issues of enhancing the accuracy of spatial prediction of

7、 soil from the use of geostatistics methods and auxiliary variables. Main contents and conclusions are as follows:Soil organic matter, soil water content, soil available potassium and manganese have been examined here from Beijing Daxing district with different sample sizes to compare the performanc

8、e of the MoM and REML. The results shows the size of data set have a marked impact on the accuracy of variogram. When the sample is less than 70 points, the prediction accuracy of target variables are beginning to decline significantly. The REML variograms are more precise than the MoM variograms in

9、 the cases of small size of sample.To assess the impact of spatial scale on variogram, the performance of MoM variogram and REML variogram were examined with the data of soil organic matter and soil available manganese from different scale plots. The results showed that the spatial scale had a obvio

10、us effect on the accuracy of variogram computed by both MoM and REML. In this paper, with soil organic matter and soil manganese as the target variable, REML variograms are more accurate than MoM variograms with the 20km scale and 30km scale plots.This paper compares the accuracy of predictions of t

11、he soil organic matter, soil available potassium, phosphorus, iron and other soil variables in Pinggu District from ordinary kriging, cokriging and regression kriging. For the last two, intensive elevation data were used as auxiliary variable. The results showed that elevation can be used to improve

12、 the accuracy of spatial prediction of soil properties, but this method is not suitable for all soil variables. The results of factor analysis suggested that linear relationship, structure relationship between elevation and soil properties and the global trend of soil properties should be carefully

13、checked before the selection of the most appropriate method.The performance of the empirical best linear unbiased predictor with residual maximum likelihood (BLUP-REML) with that of regression kriging (RK) were compared for predicting soil organic matter (SOM) with the presence of different external

14、 drifts. Terrain attributes (elevation, slope and topographic wetness index) calculated from DEM was used as external drift variable. The results suggest that topographical data can further improve the accuracy of the spatial predictions of SOM by using RK and REML. Although the improvement with REM

15、L-EBLUP is not noticeable compared with RK, the prediction capability of REML-EBLUP is more stable than RK especially when a smaller proportion of the variation in target variable is accounted for by a trend model.The sequential Gaussian simulation (SGS) and sequential Gaussian co-simulation (SGCS)

16、algorithms were compared with respect to the accuracy of predictions as well as to the uncertainty inherent in the spatial prediction of soil organic matter. The SGS algorithm accounted for only the SOM data. The SGCS accounted for both SOM data and intensive elevation data. The root mean square err

17、ors revealed that the more accurate simulations were those accounting for intensive elevation information by the SGCS method for the two areas compared with SGS. As regards modeling local uncertainty, SGCS performed better at modeling prediction uncertainty than SGS. In addition, the results of asse

18、ssing the standard deviation confirmed that the exhaustive elevation data could be used to reduce the spatial uncertainty of SOM by SGCS compared with the SGS algorithm.Keywords: soil properties, accuracy, uncertainty, geostatistics, residual maximum likelihood, elevation 93目 錄第一章 緒論11.1選題背景與研究意義11.

19、2 研究現(xiàn)狀綜述11.3主要研究內(nèi)容61.4 技術(shù)路線7第二章 樣本點數(shù)對變異函數(shù)準確性的影響92.1 引言92.2 研究目的92.3 理論與方法92.4 研究區(qū)概況與數(shù)據(jù)獲取112.5 結(jié)果分析132.6本章小結(jié)29第三章 空間尺度對變異函數(shù)的影響313.1 前言313.2 研究目的313.3 研究區(qū)介紹與數(shù)據(jù)獲取313.4 結(jié)果分析333.5 本章小結(jié)40第四章 利用高程輔助提高土壤屬性空間預測精度424.1 引言424.2 研究目的424.3 理論與方法434.4 研究區(qū)概況與數(shù)據(jù)獲取454.5 結(jié)果分析474.6 本章小結(jié)56第五章 利用不同外部趨勢變量預測土壤有機質(zhì)空間分布575.1

20、 引言575.2 研究目的575.3 理論與方法575.4 研究區(qū)概況與數(shù)據(jù)獲取595.5 結(jié)果分析615.6 本章小結(jié)70第六章 高程輔助變量在不確定性分析中的應用726.1 引言726.2 研究目的726.3 理論與方法726.4 研究區(qū)概況與數(shù)據(jù)獲取756.5結(jié)果分析766.6 本章小結(jié)82第七章 主要結(jié)論與研究展望837.1 主要結(jié)論837.2 研究展望84參考文獻85致 謝90個人簡歷91中國農(nóng)業(yè)大學博士學位論文 第一章 緒論第一章 緒論1.1選題背景與研究意義土壤作為地球表層系統(tǒng)自然地理環(huán)境的重要組成部分，在全球氣候變化、資源與環(huán)境、生態(tài)多樣性、食品安全以及人口與耕地等一系列問題中

21、起著十分重要的作用。例如，土壤產(chǎn)生與排放CO2、CH4和N2O等溫室氣體的過程，是陸地生態(tài)系統(tǒng)碳氮循環(huán)的一個重要過程，是土壤碳氮庫的主要輸出途徑，對全球氣候變化，尤其是全球變暖起著重要的影響作用。隨著上述熱點問題的逐步深入研究，對相關(guān)的土壤物理、化學和生物過程的理解和認識顯得越來越重要，因此對土壤屬性空間變化的信息，不論在信息量方面還是在信息的準確性方面，都提出了越來越高的要求。例如，在區(qū)域土壤水分和溶質(zhì)運移動力學模擬研究中，需要區(qū)域內(nèi)每個具體點位的土壤飽和與非飽和導水率和土壤水分特征曲線等土壤水力學參數(shù)作為模型運算的基礎(chǔ)；在現(xiàn)代精準農(nóng)業(yè)中，為了實現(xiàn)田間的點位管理，就必須對田間具體點位的土壤養(yǎng)

22、分信息有一個清楚的認識。然而，土壤作為母質(zhì)、地貌、氣候、生物和時間等成土因素相互作用的產(chǎn)物，土壤屬性在空間和時間尺度上表現(xiàn)出復雜的差異性。由于土壤的空間變異特性，在目前相關(guān)土壤研究工作中經(jīng)常面臨的一個主要問題是缺少必要的、準確的和足夠的土壤屬性空間分布信息。田間調(diào)查采樣及隨后的室內(nèi)分析是目前唯一能夠獲取土壤準確信息的方法。土壤的田間調(diào)查一方面需要大量的人力、物力和時間，另一方面獲取的土壤信息卻十分有限，不能滿足各方面對土壤信息的需求。為了有效的解決土壤信息的供求矛盾，同時也為了客觀真實的反映區(qū)域內(nèi)總的土壤變化特征，就需要通過一定的方法利用田間調(diào)查的準確土壤信息對未采樣點的土壤屬性值進行預測。在

23、提高土壤屬性空間預測精度方面，目前地統(tǒng)計學中備受關(guān)注的兩條主要途徑是：（1）利用有限最大似然法估計變異函數(shù)；（2）利用在空間上與土壤相關(guān)的變量來輔助預測。本文將從樣本點數(shù)、空間尺度以及外部趨勢度等方面對上述兩條途徑的應用進行探討。研究結(jié)果將有有助于對土壤屬性空間預測精度的預測方法體系有更全面的認識，對正確認識土壤空間變化提供更完善的理論基礎(chǔ)。這對于當前我國開展的測土施肥和耕地質(zhì)量動態(tài)監(jiān)測工作具有重要的實踐意義。1.2 研究現(xiàn)狀綜述為了獲取更準確的土壤屬性預測值，土壤學家在土壤空間預測方面一直努力尋求更好的預測方法，一些先進的預測方法不斷地被引進到土壤學中。例如，Burgess and Webs

24、ter （1980）將普通克里格方法應用到土壤屬性的空間預測；Lark and Webster（1999）將小波分析方法應用到土壤學中進行尺度相關(guān)的分析；Christakos（1998）首先用貝葉斯最大熵方法解決土壤屬性時空變化問題。在這些方法中，地統(tǒng)計法方法得到了廣泛地關(guān)注和研究。地統(tǒng)計學起源于20 世紀60 年代，早期主要應用于采礦業(yè)和石油勘探中。地統(tǒng)計學的研究對象是空間變量（也叫區(qū)域化變量）。空間變量的特點是變量在空間上不一定是完全隨機的或完全獨立的，即在空間上可能是相互聯(lián)系的。地統(tǒng)計學將空間變量在空間某一位置的測量值看作是滿足固有假設(shè)條件的隨機函數(shù)在這個位置上的一次實現(xiàn)值。由于不同位置

25、的樣本值之間可能是相互聯(lián)系的，就需要分析樣本的空間位置是否含有必要的信息，來揭示變量性質(zhì)隨距離的相關(guān)關(guān)系，即變量的空間變異結(jié)構(gòu)(張仁鐸, 2005)。有關(guān)地統(tǒng)計學的基本理論在此不再累述，具體內(nèi)容參考相關(guān)文獻（Goovaerts，1997；Stein，1999；王政權(quán)，1999；Webseter and Oliver，2001；Wachernagel，2003；張仁鐸，2005；）。如今，越來越多涉及到空間分析的學科求助于地統(tǒng)計學的研究工具，地統(tǒng)計學已經(jīng)被廣泛用于土壤學、生態(tài)學、環(huán)境科學、地理學等諸多領(lǐng)域的研究中。自從Burgess and Webster（1980）將地統(tǒng)計學引入土壤學后，地統(tǒng)

26、計學已經(jīng)在土壤學的各個領(lǐng)域得到廣泛地研究與應用。例如，土壤養(yǎng)分方面（胡克林等，1999； 張世熔等，2003；許文強等，2006；王淑英等，2008）、土壤有機質(zhì)方面（張世熔等，2002；黃元仿等，2004）、土壤質(zhì)地方面（高峻等，2003）、土壤污染方面（楊勁松等，2006； 張紅艷等，2006；鐘曉蘭等，2007）、土壤鹽分方面（胡克林，2000；楊勁松等, 2006）、土壤水分方面（李笑吟等, 2006; 楊勁松等, 2006）等。上述文獻的側(cè)重點都是利用克里格方法對各土壤屬性的空間分布進行插值預測，然后對土壤屬性的預測結(jié)果進行評價。在應用克里格法對土壤屬性進行預測的過程中，還需要對克里

27、格方法的適用范圍、克里格法預測能力表現(xiàn)以及克里格法預測結(jié)果存在的誤差等問題有所了解和認識。相關(guān)文獻對這些問題進行了探討和研究。下面從變異函數(shù)準確性、地統(tǒng)計混合模型（geostatistics hybrid method）結(jié)合輔助變量的應用、預測結(jié)果不確定性分析三個方面進行闡述。1.2.1關(guān)于變異函數(shù)準確性的研究地統(tǒng)計方法以區(qū)域化變量理論為基礎(chǔ)，將土壤特性的空間變化作為一個隨機函數(shù)的實現(xiàn)。用來表征隨機變量空間結(jié)構(gòu)的變異函數(shù)是地統(tǒng)計學的核心。傳統(tǒng)的估計變異函數(shù)的方法是Matheron（1965）提出的矩量法（method of moment, MoM）。MoM法估計變異函數(shù)包括兩個階段：第一，利用

28、樣本數(shù)據(jù)對不同分離距離計算樣本變異函數(shù)；第二，對樣本變異函數(shù)進行擬合得到一個連續(xù)的變異函數(shù)。為了準確地估計變異函數(shù)，用來計算樣本變異函數(shù)值的數(shù)據(jù)量應足夠大。Webster and Oliver（1992）根據(jù)Monte Carlo方法決定的置信限顯示至少需要100150個樣本數(shù)據(jù)，利用MoM法估計的變異函數(shù)才具有可靠性。但是，在實際應用中，經(jīng)常出現(xiàn)樣本數(shù)據(jù)比較少的情況。因此，需要一種方法能夠利用比較少的樣本數(shù)據(jù)估計出比較可靠的變異函數(shù)。Pardo-Iguzquiza（1998）提出用最大似然法（maximum likelihood, ML）估計變異函數(shù)是另外一種可以用于土壤學的方法。ML是一種

29、估計空間變量協(xié)方差參數(shù)的方法。Mardia（1984）首先在回歸分析中用ML估計隨機變量的空間協(xié)方差模型。隨后，Kitanidis（1987）將ML法明確地用于地統(tǒng)計分析。Lark（2000）分別通過模擬數(shù)據(jù)和田間實測數(shù)據(jù)對MoM法和ML法進行了比較，模擬數(shù)據(jù)的結(jié)果顯示ML法在采樣點數(shù)為60的時候的結(jié)果與MoM法在采樣點數(shù)為90到120點的結(jié)果一樣，但是田間實測數(shù)據(jù)的結(jié)果卻顯示兩種方法之間沒有明顯區(qū)別。相關(guān)文獻(Cressie, 1993; Kitanidis and Lane, 1985) 指出由于ML法是同時估計趨勢參數(shù)和協(xié)方差參數(shù)，導致了ML法估計的方差參數(shù)是存在偏差的。Patterso

30、n and Thompson（1971）提出的REML法克服了ML法的缺點。REML首先對原始數(shù)據(jù)進行線性合并為廣義增量（generalized increments），將趨勢去除，然后再通過最大似然函數(shù)估計協(xié)方差參數(shù)。REML在統(tǒng)計理論上是最合理的一種方法，最近已經(jīng)在計量土壤學中引起了關(guān)注。Lark and Cullis（2004）詳細的描述了如何利用REML法估計土壤屬性的協(xié)方差參數(shù)，并建議在土壤屬性空間預測過程中使用REML法。Kerry and Oliver（2007a）通過四個不同研究區(qū)域評價了樣本點數(shù)量對REML法在估計變異函數(shù)方面的影響，結(jié)果顯示，當樣本點數(shù)小于100時，利用RE

31、ML法變異函數(shù)的預測結(jié)果比利用MoM法變異函數(shù)的預測結(jié)果要準確。Minasny and McBratney（2007）根據(jù)他們的研究結(jié)果指出，當樣本點個數(shù)小于200時，就可以考慮使用REML法分析變異函數(shù)分析。1.2.2地統(tǒng)計混合模型方法結(jié)合輔助變量在土壤屬性空間預測中的應用在地統(tǒng)計學中，土壤屬性變量被作為一種區(qū)域化變量來分析，土壤屬性的空間變異性可由兩部分來代表，即空間確定性部分和隨機部分，也稱趨勢成分和隨機成分(張仁鐸, 2005)。在普通克里金法中，假定確定性部分在空間上常量，主要是估計隨機部分，在實際應用中，這些假定可能會限制克里金法的應用。例如，由于地形地貌的影響，土壤特征可能在研究

32、區(qū)內(nèi)表現(xiàn)出明顯的趨勢，在這種情況下，確定性部分在空間上不再是常量，也就不能滿足普通克里格法對數(shù)據(jù)的要求。泛克里格法（universal kriging：UK）是一種用來分析具有空間趨勢變量的最佳線性估值方法。這種方法通常假定土壤變量的空間確定性部分（稱為空間趨勢）能夠用空間坐標的多項式函數(shù)來表示。這樣就可以首先通過空間坐標預測未采樣點的土壤趨勢部分，再利用普通克里格法預測未采樣點上的隨機部分，最后將空間趨勢部分和隨機部分進行疊加得到未采樣點的土壤預測值。但是，用空間坐標來預測土壤變量的空間趨勢也存在一定的局限性。外部趨勢克里格法（kriging with external drift： KED

33、）和回歸克里格（regression kriging：RK）作為泛克里格法的有效的推廣形式，可以利用另外一種輔助變量的函數(shù)來表示土壤變量的空間趨勢部分，彌補了UK法的缺陷。通常要求輔助變量具有以下特點，（1）獲取比較容易或者測量費用相當廉價；（2）在研究區(qū)內(nèi)信息量特別豐富；（3）與目標土壤變量存在一定的相關(guān)關(guān)系。RK法和KED法在數(shù)學推導公式上是相同的，兩者的區(qū)別在于：在KED法中用來描述趨勢的線性函數(shù)的系數(shù)在局部搜尋鄰域內(nèi)被看作常數(shù)，他們在克里格預測過程中被隱含地估計；而在RK法中，估計趨勢線性函數(shù)的系數(shù)是獨立與克里格法過程的。另外一種可以利用輔助變量的方法是多元克里格方法中的協(xié)克里格法（c

34、okriging：CK），CK法是普通克里格的拓展，考慮了空間兩個或多個變量的空間交叉相關(guān)關(guān)系。有關(guān)UK法、KED法、RK法和CK法的基礎(chǔ)理論在此不再詳述，具體參考相關(guān)文獻（Odeh et al., 1995；Hengl et al., 2004，2007；Goovaerts，1997；Webseter and Oliver，2001；Wachernagel，2003）。關(guān)于上述幾種地統(tǒng)計混合模型方法的預測能力和適用范圍，許多文獻通過不同的土壤變量和輔助變量對這幾種方法進行了比較。例如，Meul and Van Meirvenne（2003）通過不同類型的非平穩(wěn)土壤質(zhì)地數(shù)據(jù)和輔助變量高程，對O

35、K、CK、UK和RK法進行了比較，結(jié)果顯示CK法和UK法可以很好的解釋局部非平穩(wěn)性。Simbahan et al（2006）利用表面反射率、土壤電導率和相對高程等多個外部輔助變量，通過CK、KED和RK三種方法，對土壤有機碳進行空間預測，結(jié)果表明，在這幾種地統(tǒng)計混合方法，RK法是一種應用靈活且預測精度較高的方法。Baxter and Oliver（2005）以高程為輔助變量，通過CK、RK和KED三種方法，對兩個研究區(qū)的土壤礦物氮和潛在有效氮進行分析，結(jié)果表明KED法和CK法都具有比較好的預測結(jié)果。從上述文獻的結(jié)果可以看出，對于不同的土壤變量，不同的研究區(qū)域和不同的研究尺度，沒有一種最好的預測

36、方法，每種方法都有一定特點和適用范圍。下面簡要論述了每種方法的應用條件及其不足之處。CK法在理論上是一種比較完善的多元克里格法。CK法要求的條件為：（1）目標變量與輔助變量都滿足固有假定條件；（2）目標變量與輔助變量在空間上是相聯(lián)系的，即在同一位置上，目標變量的取值與輔助變量的取值是相關(guān)的（Frogbrook and Oliver, 2001）；（3）目標變量的樣本點比較少，而輔助變量的采樣點則比較多。CK法在實際應用中存在兩個方面的不足：（1）CK法要求目標變量的變異函數(shù)、輔助變量的變異函數(shù)以及二者之間的協(xié)變異函數(shù)是由共同的基本變異函數(shù)模型線性組合而成（即：協(xié)區(qū)域化線性模型），在協(xié)區(qū)域化線性

37、模型的擬合過程中受到的限制條件比較多（Wachernagel, 2003）；（2）CK法在預測過程中還存在計算量大、計算時間長的問題(Goovaerts, 1997)。因此，在應用中CK法并不是一種實用的方法。KED法的應用條件是：（1）輔助變量的空間趨勢與目標變量的空間趨勢是相關(guān)的；（2）目標變量去除趨勢后的殘差能夠用地統(tǒng)計方法模擬；（3）目標變量與輔助變量是線性相關(guān)的，輔助變量必須在所有目標變量的采樣點上是已知的，并且在預測格網(wǎng)的所有點上也是已知的。KED法的缺陷是：由于趨勢系數(shù)的估計是在克里格法估值過程中同時進行的，當輔助變量與目標變量的關(guān)系不是簡單的線性關(guān)系時，KED法的應用受到限制（

38、Deutsch and Journel, 1998）。RK法的應用條件與KED法相同，但是RK法明確地把趨勢估計與殘差插值分離，這樣就可以利用任意復雜的回歸形式，比較簡單實用，因此RK法是一種應用最廣泛的方法（Baxter and Oliver, 2005; Finke et al., 2004; Herbst et al., 2006; Simbahan et al., 2006），而且許多文獻（Hong et al., 2005; Odeh et al., 1995; Simbahan et al., 2006）也報道了RK法的預測精度要高于OK法、CK法和RK法等方法的預測精度。但是，C

39、ressie（1993）和Lark et al.（2006）指出，RK法也存在兩方面的缺陷：（1）趨勢一般通過普通最小二乘法（OLS）來估計，雖然這個估計是無偏差的，但是這個估計沒有產(chǎn)生最小方差估計；（2）利用回歸殘差估計的變異函數(shù)在理論上是有偏差的，而且偏差隨著分離距離的加長而逐漸變大。RK法的缺陷導致兩方面的結(jié)果，一方面是這種偏差將導致隨機變量總的變化被低估，另一方面是擬合的變異函數(shù)不能正確的表達隨機變量的空間相關(guān)性。鑒于RK法存在上述的缺陷，Lark et al.（2006）建議采用統(tǒng)計上更穩(wěn)健的一種方法，即以有限最大似然法為基礎(chǔ)的線性無偏最優(yōu)估計方法（BLUP-REML），來預測具有空

40、間趨勢的土壤變量。Lark (2006)指出BLUP-REML法能夠無偏的估計殘差的變異函數(shù)，能夠得到最小方差的趨勢估計，還能夠得到未采樣點上估計的趨勢成分和隨機成分的和，而且還知道預測方差。Minasny and McBratney（2007）利用了不同研究區(qū)的不同土壤屬性數(shù)據(jù)對該方法進行了檢驗，結(jié)果顯示當樣本點個數(shù)小于200時，BLUP-REML法是一種有效的方法。然而，目前有關(guān)BLUP-REML法在土壤屬性的空間預測方面的具體案例還比較少，對BLUP-REML法的預測能力及影響因素還需要進一步檢驗與評價。隨著科學技術(shù)的發(fā)展，在空間上與土壤相關(guān)的輔助變量越來越容易獲取，這些輔助變量主要包括

41、數(shù)值高程模型（DEM）、遙感影像數(shù)據(jù)和電磁波（EMI）掃描數(shù)據(jù)等。利用這些信息量豐富的數(shù)據(jù)來輔助提高土壤變量的空間預測精度，在計量土壤學土壤學中已經(jīng)引起極大的關(guān)注和研究。例如，Odeh et al.（1995）利用高程數(shù)據(jù)通過RK法提高了上層土壤沙粒和下層土壤粘粒的空間預測精度。許紅衛(wèi) 等（2007）研究了土壤養(yǎng)分(土壤全N、全P、有機質(zhì)、速效N、速效K) 與水稻各生長期水稻冠層光譜的關(guān)系, 并將光譜指數(shù)作為協(xié)因子, 進行土壤養(yǎng)分的Cokriging 插值研究，結(jié)果顯示，利用光譜指數(shù)作為協(xié)因子可以提高土壤養(yǎng)分的空間預測精度。Huang et al.（2007）利用近紅外波譜測量數(shù)據(jù)和LandS

42、at ETM影像數(shù)據(jù)，結(jié)合地形信息對土壤有機碳的空間預測進行了分析，結(jié)果證實用近紅外波譜測量數(shù)據(jù)或LandSat ETM影像數(shù)據(jù)結(jié)合地形數(shù)據(jù)對提高土壤有機碳預測精度是一種有效的方法。Liu et al.（2008）利用電導率（ECa）、遙感高光譜反射率和DEM數(shù)據(jù)對土壤滲透力（soil drainage）的空間分布進行了預測，結(jié)果表明，單獨利用地形指標作為輔助變量不能提高土壤滲透力的預測精度，但將地形指標與電導率或遙感數(shù)據(jù)結(jié)合后，就可以有效的提高土壤滲透力的預測精度。從上述可以看出，利用信息量豐富的輔助變量，可以有效地提高土壤屬性的空間預測精度。在與土壤相關(guān)的輔助變量中，地形指標變量是一類主要

43、的輔助變量，相關(guān)文獻也證實了地形指標與土壤屬性之間明顯的相關(guān)關(guān)系。More（1993）通過用地形指標作為解釋變量，發(fā)現(xiàn)坡度的傾角、地形濕度指數(shù)可以解釋土壤有機質(zhì)、pH值、有效磷和粘粒含量等土壤變量空間變化的50%；Iqbal et al.（2005）通過在逐步線性回歸分析，地形變量之間的組合能夠介紹土壤空間變化的10%62%。可以看出，地形在土壤空間分布起著重要的影響作用。1.2.3土壤屬性空間預測結(jié)果不確定性分析利用土壤采樣點數(shù)據(jù)，采用克里格等方法對未采樣點的土壤屬性值進行預測，已經(jīng)被深入研究和應用。然而，克里格估計值是有誤差的，預測結(jié)果具有一定的平滑效應，即在土壤屬性實際值比較低的點卻有一

44、個高的估計值，或者在真實值高的地方又低估了。土壤管理者如果以克里格預測結(jié)果作為施肥決策依據(jù)，就可能導致錯誤的施肥決策。例如，在土壤養(yǎng)分含量缺乏的地點應該增加施肥卻沒有進行，在不需要施肥的地方卻又增加施肥，導致資源浪費和其他環(huán)境問題。上述這種決策風險是不可避免的。如果這些風險能夠被量化，管理者將能夠?qū)@些風險進行更有效地管理。這就需要我們對未采樣點土壤屬性值進行不確定性分析，估計未采樣點土壤屬性值超過或低于某一預設(shè)閾值的概率。目前，這種不確定性分析在決策過程和模型評價方面等方面受到越來越多的關(guān)注。地統(tǒng)計學中指示克里格法在解決這些問題方面比較突出，經(jīng)常被用于土壤重金屬污染、土壤養(yǎng)分分布等的相關(guān)研究

45、中。對于某一具體位置，指示克里格法用條件累計分布函數(shù)估計某一土壤變量取值大于給定閾值的概率。李笑吟等（2006）以普通克里格法和指示克立格法為基本方法,研究土壤水分空間分布及土壤水分有效性水平在空間上的概率分布。徐英等（2006）運用非參數(shù)地質(zhì)統(tǒng)計學的多元指示克立格法, 結(jié)合單變量指示克立格法, 對黃河河套灌區(qū)的鹽漬化土地上兩個比較關(guān)鍵時期的土壤水分和鹽分進行空間分布的分析, 給出了同一時期土壤水鹽和不同時期土壤鹽分滿足一定標準的綜合概率圖。然而，指示克里格法以克里格估值作為基礎(chǔ)，使這種方法不可避免的表現(xiàn)出克里格插值方法的缺點?？死锔穹ü烙嬵A測結(jié)果的不確定性主要是通過計算預測值和相應的誤差方差

46、，然后對二者進行計算構(gòu)造高斯置信區(qū)間(Cinnirella et al., 2005)。這一過程存在三方面的不足，（1）克里格法的預測結(jié)果存在平滑效應，把目標變量空間變化的局部細節(jié)給平滑掉了，通常表現(xiàn)為把應該為高值的地方給估低了，把應該為低值的地方估高了；（2）克里格預測誤差方差與實際預測值沒有關(guān)系，它只依賴于采樣數(shù)據(jù)的空間分布和變異函數(shù)模型；（3）計算置信區(qū)間時，需要假設(shè)預測誤差為正態(tài)分布，然而土壤變量和其它環(huán)境變量經(jīng)常呈偏態(tài)分布，這一假設(shè)通常很難滿足(Bourennane et al., 2007)。由此可以看出，克里格法用估計值和估計誤差構(gòu)造置信區(qū)間并不是一個滿意的解決方案。條件模擬方法

47、可以克服克里格法在估計不確定性方面的缺陷。條件模擬的主要思想是通過產(chǎn)生一組等概率的土壤屬性空間分布圖，然后根據(jù)所模擬的這些結(jié)果圖間的差別來評價預測值的不確定性，而且這種方法還可以很好地再現(xiàn)原始變量的空間結(jié)構(gòu)。條件模擬還可以同時考慮幾個位置的空間不確定性。這一點也是指示克里格法的缺陷。所以條件隨機模擬方法在土壤及其他學科中的應用受到越來越多的關(guān)注。例如，李保國（2002）等利用序貫高斯模擬對農(nóng)田土壤表層飽和導水率進行了模擬，結(jié)果表明條件模擬與實測數(shù)據(jù)具有相同的統(tǒng)計分布特征及相同的空間結(jié)構(gòu)。Zhao et al.（2005）利用條件隨機指示模擬方法對河北省土壤有機質(zhì)碳密度的空間預測結(jié)果的不確定性進

48、行了定量地評價。史舟（2007）以杭嘉湖平原為樣區(qū)，利用普通克立格法和序貫高斯模擬方法對土壤Cu 的空間分布進行估值和模擬，并利用序貫指示模擬進行不確定性評價。結(jié)果表明由普通克里格法得到的土壤Cu的空間分布整體比較連續(xù)，具有明顯的平滑效應，估值結(jié)果數(shù)據(jù)的分布頻率趨于平緩；序貫高斯模擬結(jié)果整體分布相對離散，突出了原始數(shù)據(jù)分布的波動性，其模擬結(jié)果數(shù)據(jù)的分布頻率相對集中?？偟膩砜?，在進行空間預測和不確定性分析時隨機模擬方法比克里格方法顯示出更多的優(yōu)點，因此在土壤學科中得到越來越廣泛地應用，尤其是對于研究那些具有空間不確定性且會對環(huán)境帶來不良影響的風險性變量，更具有實際意義。1.3主要研究內(nèi)容從上述研

49、究現(xiàn)狀可以看出，目前提高土壤屬性空間預測精度的基本途徑主要有兩方面：利用更完善的預測方法和利用信息量更豐富的輔助變量。在預測方法方面，將REML法估計變異函數(shù)與地統(tǒng)計方法結(jié)合引入到土壤屬性空間預測領(lǐng)域時間并不長，相關(guān)研究還比較少。還有許多問題需要進一步的探討與研究。同樣，在利用輔助變量提高土壤空間預測方面也存在一些問題需要回答。為此，本文從預測方法和利用輔助變量兩個方面對以下主要內(nèi)容進行研究和闡述。（）評價樣本點數(shù)對變異函數(shù)準確性的影響。變異函數(shù)是對目標土壤變量空間變異結(jié)構(gòu)的一種表現(xiàn)形式。一方面不同的樣本點數(shù)影響變異函數(shù)的估計結(jié)果，另一方面，目標變量本身的內(nèi)在特性應該對變異函數(shù)具有更重要的影響

50、作用。本文將以同一研究區(qū)的土壤有機質(zhì)、土壤含水量、土壤速效鉀和土壤有效錳四個土壤屬性作為研究目標變量，通過選擇不同的樣本點數(shù)，對比分析MoM法和REML法在計算變異函數(shù)準確性方面的表現(xiàn)。（2）評價空間尺度對計算變異函數(shù)的影響。土壤的空間變化存在于不同的空間尺度上，而不同空間尺度上具有不同的主要控制因素和過程。當我們用變異函數(shù)來描述土壤變量空間變異結(jié)構(gòu)的時候，同樣需要考慮土壤變化所處的空間尺度。本文以土壤有機質(zhì)和土壤有效錳為目標變量，在保持樣本點個數(shù)不變的情況下，通過選擇不同的空間尺度樣區(qū)，對比分析MoM法和REML法在計算變異函數(shù)準確性方面的表現(xiàn)。（3）利用高程作為輔助變量提高土壤屬性空間預測

51、精度。本文將以高程數(shù)據(jù)作為輔助變量，以土壤有機質(zhì)、土壤速效鉀、土壤有效磷、土壤有效鐵等為目標變量，通過利用OK法、CK法和RK法三種方法，對比分析高程數(shù)據(jù)是否可以用來提高上述土壤目標變量的空間預測精度。探討土壤養(yǎng)分變量與高程數(shù)據(jù)之間應該滿足哪些條件，就可以明確地利用高程數(shù)據(jù)來輔助提高土壤變量的空間預測精度。（4）評價RK法和BLUP-REML法在不同外部趨勢度下的表現(xiàn)。地形在土壤的空間分布中起著重要的作用，本文將利用高程、坡度和地形濕度指數(shù)等地形指標作為輔助變量，通過組合成不同的外部趨勢模型，利用RK法和BLUP-REML法來預測土壤有機質(zhì)的空間分布。比較RK法和BLUP-REML法在不同程度

52、的外部趨勢情況下預測能力的表現(xiàn)，評價高程、坡度和地形濕度指數(shù)等地形指標作為輔助變量在土壤有機質(zhì)空間預測中的作用。（5）高程輔助變量在模擬不確定性方面的應用。本文將高程輔助數(shù)據(jù)與序貫高斯協(xié)模擬方法相結(jié)合，對土壤有機質(zhì)的空間分布進行模擬，對預測結(jié)果的準確性、局部不確定性模擬的準確性和空間不確定性模擬的準確性進行評價。探討在條件模擬計算環(huán)境下高程是否有助于提高土壤有機質(zhì)的空間預測精度，是否有助于提高預測結(jié)果不確定性模擬的準確性。1.4 技術(shù)路線本文將按照圖1-1所示的技術(shù)路線對上述主要內(nèi)容進行研究：不同土壤屬性變異函數(shù)準確性比較分析不同土壤屬性同一輔助變量同一土壤屬性不同輔助變量MoM法OK法REM

53、L法CK法RK法BLUP-REML法序貫高斯模擬序貫高斯協(xié)模擬不同樣本點數(shù)不同空間尺度利用輔助變量比較分析局部不確定性對比分析輔助變量土壤屬性空間預測精度分析預測結(jié)果不確定性分析參考方法空間不確定性對比分析土壤屬性圖1-1 本研究技術(shù)路線總體框圖Fig. 1-1 Global framework for the research中國農(nóng)業(yè)大學博士學位論文第二章 樣本點數(shù)對變異函數(shù)準確性的影響第二章 樣本點數(shù)對變異函數(shù)準確性的影響2.1 引言土壤屬性的變異函數(shù)作為其空間變異結(jié)構(gòu)的一種表現(xiàn)形式，一方面受土壤屬性本身的內(nèi)在特征決定，另一方面受樣本點數(shù)和變異函數(shù)估計方法的影響。關(guān)于對樣本點數(shù)最小需求這一

54、問題，不同的文獻有不同的結(jié)論。Journel and Huijbregts (1980)指出一般必須大于30個樣本點，而Webster and Oliver（1992）則建議至少需要100150個樣本數(shù)據(jù)MoM法估計的變異函數(shù)才具有可靠性。REML法在理論上是一種更合理的計算變異函數(shù)的方法，Kerry and Oliver（2007a）的研究結(jié)果表明，當樣本點數(shù)小于100時利用REML法變異函數(shù)的預測結(jié)果比利用MoM法變異函數(shù)的預測結(jié)果要準確，而Minasny and McBratney（2007）指出當樣本點個數(shù)小于200時，就可以考慮使用REML法來代替MoM法分析變異函數(shù)分析。然而，對于

55、REML法來說，至少需要多少個樣本點數(shù)其變異函數(shù)才具有可靠性，上述文獻也沒有明確結(jié)論。由于在土壤學中REML法還是一種比較新的方法，相關(guān)研究還比較少，仍然需要更多的土壤屬性研究實例對REML法的適用范圍和預測能力進行評價和研究。2.2 研究目的本章將以同一研究區(qū)的土壤有機質(zhì)、土壤含水量、土壤速效鉀和土壤有效錳四個土壤屬性作為研究目標變量，通過選擇不同的樣本點數(shù)，對比分析MoM法和REML法在計算變異函數(shù)準確性方面的表現(xiàn)。本章的目的主要有以下三個方面：（1）進一步評價土壤樣本點數(shù)對變異函數(shù)準確性的影響；（2）比較不同土壤變量隨樣本點數(shù)變化其變異函數(shù)準確性的變化方式；（3）比較MoM法和REML法

56、計算的變異函數(shù)的準確性隨著樣本點數(shù)變化的具體區(qū)別。2.3 理論與方法地統(tǒng)計分析的核心工具是變異函數(shù)。變異函數(shù)用來表征隨機變量的空間變異結(jié)構(gòu)。用來計算變異函數(shù)的方法有多種。目前最常用的估計變量函數(shù)的方法是MoM法。另外一種最新引進地統(tǒng)計的參數(shù)化估計方法是REML法。下面分別介紹兩種方法的主要步驟及相關(guān)數(shù)學推導公式，有關(guān)這兩種方法的詳細理論參考相關(guān)文獻(Pardo-Iguzquiza, 1997; Webster and Oliver, 2001)。2.3.1 MoM法估計變異函數(shù)變異函數(shù)與一定分離距離的兩個隨機變量取值的平方差相關(guān)。實際應用中，通過一定范圍分離距離（例如1020km）內(nèi)所有數(shù)據(jù)對

57、的平方差的平均值與各對樣本數(shù)之間的距離的平均值來計算樣本變異函數(shù)。MoM法首先根據(jù)樣本數(shù)據(jù)計算樣本變異函數(shù)值。設(shè)是土壤變量在個位置的樣本數(shù)據(jù)值系列，在固有假設(shè)條件下，即： （2-1）樣本的變異函數(shù)通過下式來計算： （2-2）式中代表樣本變異函數(shù)，代表變異函數(shù)模型，是分離距離，是在之間用來計算樣本變異函數(shù)值的樣本的對數(shù)。根據(jù)目標變量的采樣點數(shù)據(jù)利用公式（2-2）計算不同分離距離的樣本變異函數(shù)值，然后利用不同分離距離的樣本變異函數(shù)值，通過加權(quán)最小二乘法來擬合適當?shù)淖儺惡瘮?shù)模型。最常用的變異函數(shù)理論模型有球狀模型（spherical）、指數(shù)模型（exponentail）和高斯模型（gaussian）等。有關(guān)這些模型詳細介紹參考相關(guān)文獻(Goovaerts, 1997; 張仁鐸, 2005) 。變異函數(shù)模型的參數(shù)就可以直接參與到克里格計算中。在選擇分離距離是有兩種情況：（1）對于不規(guī)則的空間數(shù)據(jù)，需要以實際分離距離為中心定義一個分離距離等級，這個等級包括一定范圍內(nèi)的距離和方向。分離距離等級的選擇對樣本變異函數(shù)的影響至關(guān)重要，如果分離距離等級太短，變異函數(shù)比較嘈雜，如果分離距離太長，變異函數(shù)太光滑；（2）對于格網(wǎng)數(shù)據(jù)，格網(wǎng)間距通常被用來作為分離距離等級。2.3.2 REML法估計變異函數(shù)Patterson and Thompson (1971)提出的有限最大似然法