最近的幾篇頭條文章，分别從多光譜與高光譜遙感的實際應用出發，對影像前期處理與相關算法、反演操作等加以詳細介紹。而通過遙感手段獲取了豐富的各類地表信息數據後，如何對數據加以良好的數學處理與科學分析，同樣是我們需要重視的問題。因此，準備由這一篇文章入手，新建一個專欄，逐篇地對地學計算方面的内容加以初步總結。

  那麼首先，我們就由地學計算的幾個基本概念入手，對相關理論方面的内容加以一定了解。

  需要注意的是，以下内容如果單獨來看或許有些不好理解，但一旦将其與實際應用結合，便會豁然開朗。其中，具體的實際應用部分我将會在後面的文章中涉及。

  空間數據的獲取是進行空間分析的基礎與起源。為了提高研究結論的精度，我們亦總是希望能夠獲取研究區域内更多、更全面的精确空間屬性數據信息。然而，在實際研究、工作中，由于人力、成本、資源等外部條件限制，我們不可能對全部未知區域加以采樣與測量，而往往隻能得到研究區域内有限數量的采樣點及其相關屬性數據。因此，往往可以考慮選取合适的空間采樣點，利用一定數學模型，依據已知采樣點各自對應屬性數據對研究區域所有位置的未知屬性信息加以預測。

  空間插值（Spatial Interpolation）即可以實現這一需求。其是一種将離散采樣點測量數據轉換為連續數據曲面的常用方法，包括内插（Interpolation）和外推（Extrapolation）兩種應用形式。一般地，對樣本點範圍以内（即所有采樣點最大外接矩形内部）的空間進行插值，才可稱作“内插”（部分文獻亦直接用“插值”代替“内插”）；反之則稱“外推”或“預測”，往往認為外推的結果誤差較大。

  空間插值理論及其方法基于著名的“地理學第一定律（Tobler's First Law of Geography）”，即一般地，距離越近的地物具有越高的相關性。這一至今已産生深遠影響的地學定律最早由美籍瑞士地理學家Waldo R. Tobler教授于1970年提出。

  在各方法所對應的數學計算原理層面，空間插值一般可以分為确定性插值方法（Deterministic Interpolation）與地統計插值方法（Geostatistics，亦稱非确定性插值方法）兩種。其中，确定性插值方法基于研究區域内各信息點之間相似程度或整個曲面的平滑程度，從而創建連續的拟合曲面；其依據插值計算時納入考慮的采樣點分布範圍，又可進一步分為整體插值法與局部插值法。地統計插值方法則是基于研究區域内各信息點的綜合統計學規律，以變異函數（Variogram）理論與結構分析為基礎，實現其屬性的空間自相關性定量化，從而創建得出連續插值曲面。

  在所創建連續插值表面通過全部采樣點的與否層面，空間插值一般又可以分為精确性插值與非精确性插值兩種。其中，前者預測樣點的屬性數值與其各自實測值相等，即其采樣點屬性數據全部落入預測結果曲面；後者預測樣點的屬性數值則往往不與其各自實測值相等，即其采樣點屬性數據一般不會落入預測結果曲面。因此，使用非精确性插值方法往往可以避免在預測表面中出現明顯的波峰或波谷，整體呈現出平緩态勢。

  地學計算中，幾個重要假設可以說是關鍵中的關鍵；而其往往也比較難以理解，不怕，我們慢慢往下看。

  平穩假設（Stationary Assumption）是指，一組觀測值的均值是始終固定的，其與觀測值所在位置無關；将既定的某個點集由某一研究區域内某處移動至另一處時，随機函數的性質保持不變。即随機函數的分布規律不因位置的改變而改變，具有嚴格的平穩性。平穩性假設的公式表達為：

  其中，F_(x_1,⋯,x_n ) (z_1,⋯,z_n )表示位置在(x_1,⋯,x_n )上的點集(z_1,⋯,z_n )對應的随機函數。

  二階平穩性假設（Second Stationary Assumption）又稱弱平穩假設，其與協方差函數（Covariance Function）有關。這一假設認為，随機函數的均值為一常數，且任意兩個随機變量之間的協方差僅僅依賴于其二者之間的距離與方向，而與其具體位置無關。

  将上述兩個條件用公式分别表達為：

  其中，E[Z(x)]為區域化變量Z(x)的數學期望， Cov[Z(x),Z(x h)]為區域化變量Z(x)與Z(x h)所對應的協方差函數，C(h)為僅與h有關的協方差取值，m為一常數，h為滞後距。

  上述二階平穩性假設是針對整個研究區域範圍而言。若區域化變量僅僅在整個研究區域内的某個有限區域中滿足上述條件，即條件僅在局部區域成立，則稱此區域化變量滿足準二階平穩性假設（Quasi Second Stationary Assumption）。準二階平穩性假設可以視作一種折中方案，既考慮到平穩範圍的大小，又顧及到有效數據的數量。

  本征假設（Intrinsic Hypothesis）又稱内蘊假設，其與變異函數有關。這一假設認為，區域化變量的增量滿足以下兩個條件：在整個研究區域内，區域化變量增量的數學期望為0；且其方差函數存在，并隻依賴于滞後距，而與所處位置無關。

  将上述兩個條件用公式分别表達為：

  當E(x)存在時，上述第一個公式可以寫作：

  其中，Var[Z(x)-Z(x h)]為區域化變量Z(x)與Z(x h)所對應的方差函數，γ(h)為區域化變量在滞後距為h時的變異函數，m為一常數，h為滞後距，其它符号同前述意義。

  同樣的，上述本征假設亦是針對整個研究區域範圍而言。若區域化變量僅僅在整個研究區域内的某個有限區域中滿足上述條件，則稱此區域化變量滿足準本征假設（Quasi Intrinsic Hypothesis）。與準二階平穩性假設類似，準本征假設亦可視作一種折中方案，其同樣既考慮到了本征假設對應範圍的大小，又顧及到了有效數據的數量。

  此外，本征假設是地統計學中對随機函數的基本假設。

  結合上述二階平穩性假設與本征假設相關原理，可以看到兩種假設的讨論對象具有一定區别：二階平穩性假設更多是讨論某一特定研究區域内區域化變量自身【即Z(x)】的特征，而本征假設則是研究區域化變量所對應增量【即Z(x)-Z(x h)】的特征。

  一般認為，二階平穩性假設對區域化變量的要求較之本征假設更為嚴格，即若某個研究區域内的某一區域化變量滿足二階平穩性假設，則其一定滿足本征假設；反之則反，若僅知區域化變量滿足本征假設，則其不一定滿足二階平穩性假設。

  再結合平穩假設，上述三種假設的嚴格程度由大至小依次排列為：平穩假設、二階平穩性假設與本征假設。

  此外，結合二階平穩性假設的兩個條件，還可以推出協方差函數與變異函數之間的關系：

  其中，γ(h)為區域化變量所對應變異函數，C(0)為距離為0時此區域化變量所對應的協方差取值【即基台值γ(∞)】，C(h)為距離為h時此區域化變量所對應的協方差取值。

  由這一關系可知，用以衡量某一區域化變量在相距為h的兩空間位置點分别取值的自相關性的指标——協方差函數與變異函數之間具有相互關系。因此，在滿足二階平穩性假設的條件下，若協方差函數平穩，則可知變異函數平穩，即其取值隻與滞後距h有關。

  克裡格插值法需要借助空間數據的試驗變異函數及其散點圖特點，因此變異函數的計算在克裡格插值過程中發揮着重要作用；變異函數及其模型拟合對克裡格插值結果精度具有較大影響。

  變異函數（Variogram），又稱為半變異函數、半方差函數（Semi-variogram）等，其用以描述區域化變量的空間變化特征與強度，被定義為區域化變量增量平方的數學期望。

  在一維條件下，直接将區域化變量Z在位置(x)與(x h)處的取值Z(x)與Z(x h)之差的方差定義為變異函數，其因變量為距離h；而在二維或三維條件下，可以将上述一維中具有單一方向的距離h進一步引申為在任意方向α上的距離|h|。具體公式表達為如下。

  其中，γ(x,h)即為變異函數。由于公式中在其前具有一個系數“2”，因此其亦被稱作半變異函數。結合（準）二階平穩性假設、（準）本征假設等地學基本假設，變異函數取值與區域化變量樣點所處位置x無關，僅僅與樣點之間的距離h有關，則變異函數可以寫作：

  其中，〖γ(h)〗^#為區域化變量Z(x)的變異函數，N(h)為區域化變量樣點集中距離為h的點對數量，x_i為距離為h時所對應的第i個點。在這裡，距離h亦被稱作滞後距（Lag Distance）。

  一般地，區域化變量變異函數圖像往往呈現出“先快速上升，再增速減緩，後趨于平穩”的曲線特征。其具有三個十分重要的相關概念，分别為塊金常數（Nugget）、基台值（Sill）與變程（Range）。

  塊金常數代表區域化變量的随機性大小。由理論角度，在間距為0（即滞後距為零）時，區域化變量采樣點數值應當相等；而在間距無限趨近于0時，對應變異函數數值應當亦向0趨近。但是，在實際研究中，試驗變異函數在滞後距為0時，其取值并不為0，而是一個大于0的數值。這一數值便稱為塊金常數。一般地，上述塊金效應的産生可以歸因于測量誤差，或小于采樣間隔距離處的空間變化。

  基台值用以衡量區域化變量變化幅度的大小。當滞後距無限增大并到達某一程度後，試驗變異函數若趨于平穩，則這一平穩水平所對應的數值即為基台值。然而，并不是所有的區域化變量均具有基台值——如無基台值模型對應的變異函數。

  變程用以衡量區域化變量自相關範圍的大小。當滞後距無限增大并到達某一程度後，試驗變異函數若趨于平穩，則此時對應的滞後距即為變程。其中，小于變程的距離所對應的樣本位置與空間自相關，而大于變程的距離所對應樣本位置不存在空間自相關。

  此外，變異函數還有其它相關指标，如基台值與塊金常數的差值——偏基台值（Partial Sill），用以衡量空間變異性程度的塊金常數與基台值的比值——塊金系數等。

  基于不同區域化變量對應的變異函數特征，可以将其分為不同類别。依據變異函數基台值的有無，可以将模型分為有基台值模型、無基台值模型與孔穴效應模型。

  其中，有基台值模型可以依據變異函數的特征進一步分為純塊金效應模型（Pure Nugget Effect Model）、球狀模型（Spherical Model）、指數模型（Exponential Model）、高斯模型（Cubic Model或Gaussian Model）與線性有基台模型（Linear with Sill Model）等。上述幾種模型中，較為常用的模型包括球狀模型、指數模型與高斯模型。

  此外，同樣依據變異函數的特征，無基台值模型還可進一步分為線性無基台值模型、幂指數模型與對數模型等。同樣的，孔穴效應模型可分為基台值模型和無基台值模型。

  同時，針對某種區域化變量而言，其在不同方向、不同滞後距情況下可能受到不同因素影響；套合結構可以很好解決這一問題。套合結構可以表示為多個變異函數之和，每一個變異函數均代表着某種方向或某一尺度中的變異性，從而對區域化變量的特征加以更好概括。

  克裡格插值法（Kriging Method）又稱為空間局部插值法，是以上述變異函數理論及其結構分析為基礎，在有限區域内對區域化變量進行線性無偏最優估計（Best Linear Unbiased Prediction，BLUP）的一種方法，在地統計學中也被稱為空間最優無偏估計器（Spatial BLUP）。

其中，上述“線性”是指克裡格插值法對未知點屬性數值的估計采用線性估計，其公式如下：

  其中，(z_0 ) ̂是區域化變量在點(x_0,y_0 )位置處的預測值，λ_i為第i個已知點的權重系數，z_i為第i個已知點的實測值。

  上述權重λ_i可以使得各點處的預測值與實測值之間的方差最小，而求取這一權重則為克裡格插值的主要内容。

  上述“無偏”是指區域化變量在各點上估計量的數學期望等于其在同一位置上的真實值，公式如下：

  結合本征假設，無偏性可以表示為：

  上述“最優”是指區域化變量在各點上估計量與其在同一位置上的真實值的方差最小，公式如下：

  其中，上述方差被稱作估計方差或估值方差，是對估值準确程度的一種定量表示；而在克裡格插值方法中，又可以稱為克裡金方差。後續将克裡金方差記作σ_k^2。

  經過統計學相關推導，可以将克裡金方差寫作：

  由此轉換為在無偏條件約束下的最小值求解問題。引入拉格朗日乘子φ，構造拉格朗日函數：

  對權值與拉格朗日乘子求一階偏微分：

  偏微分求解結果為：

  将上述求和部分展開：

  可将上式進一步寫為矩陣相乘形式，即化為：

  其中，A代表在原有變異函數矩陣中額外添加全1行與全1列（交界處1換為0）後的矩陣，λ代表各權重組成的列向量，φ代表前述分析引入的拉格朗日乘子，B為各位置與待求解位置對應距離的變異函數值組成的列向量，且在列尾增加一個1。

  由此，即将上述函數轉化為(n 1)個未知數、(n 1)個表達式組成的方程組；通過矩陣求逆，求解方程組即可得到待求解位置與其它已知點的權重。對每一個待插點進行同樣操作，完成克裡格插值。

  正如上述分析，普通克裡格（Ordinary Kriging）方法更多依靠所選采樣點數據對空間加以插值，其插值效果較依賴于采樣點的個數、密度，以及其數據的準确性；另一方面，許多空間屬性往往受到其它環境變量的影響。例如，土壤有機碳含量與氣溫、海拔、降水量、坡度等多種環境因子在0.05水平顯著相關[3, 4]；近地面氣溫與海拔、海陸距離、NDVI等環境因子在0.01水平顯著相關[5]。由此觀之，若簡單地忽略環境要素對待插值空間屬性的影響，可能會降低最終插值結果精度。

  基于這種考慮，可以使用回歸克裡格（Regression Kriging）方法對環境因素加以考慮。應用回歸克裡格插值，首先需要确定目标變量與輔助環境變量之間的相關性；确認相關關系符合要求後，建立目标變量與環境變量的一元或多元回歸模型，從而剔除空間趨勢項。随後，依據采樣點實測數據與回歸模型計算得出的對應位置數值，求得目标變量的确定性趨勢項。運用普通克裡格方法，将殘差進行插值，并最終将回歸預測的趨勢項與普通克裡格的插值結果相加，從而得到目标變量估測值。

  最終空間插值結果的公式表達為：

  其中，Z(x_0 )為回歸克裡格估值在x_0位置的結果，m ̂(x_0 )為确定性趨勢項在這一位置的取值，e ̂(x_0 )為殘差普通克裡格插值在這一位置的取值。對研究區域内每一位置執行同樣操作，完成回歸克裡格插值。

  一般地，對于所受外界影響較大的空間屬性，回歸克裡格插值效果優于普通克裡格插值[6]。

  另一方面，如前所述，盡管回歸克裡格方法對于殘差項的估計仍采用普通克裡格插值計算，但由于其對最終目标變量插值結果的影響較之由環境因子得到的空間确定性趨勢項而言仍較低，因此回歸克裡格方法所得結果較之普通克裡格方法結果将會粗糙、破碎一些。此外，由方法考慮的因素這一角度來看，回歸克裡格方法與協同克裡格（Cokriging）方法具有一定相同之處，即二者均利用相關環境因子對目标變量的空間分布加以估測；而其不同之處在于，回歸克裡格的輔助數據為面狀分布，如栅格圖層；而協同克裡格的輔助數據為點狀分布。但無論采取何種方法，最終所得插值結果均為一個完整的表面。