克里金法

<image>

克里金或克里金法（kriging，Kriging，發音為ˈ）也稱為高斯過程回歸，在統計學中（最初在地理統計中），是一種基於受先驗協方差支配的高斯過程的插值方法。在適當的先驗假設下，克里金法在未採樣位置給出最佳線性無偏預測（best linear unbiased prediction，BLUP）。基於其他準則（例如平滑樣條）的插值方法可能無法產生最佳線性無偏預測。該方法廣泛應用於空間分析和電腦實驗領域。由於諾伯特·維納和安德烈·柯爾莫哥洛夫的貢獻，該技術也被稱為維納-柯爾莫哥洛夫預測（Wiener–Kolmogorov prediction）。

該方法的理論基礎由法國數學家喬治·馬瑟隆（Georges Matheron）於1960年建立，其依據是丹尼·G·克里格（Danie G. Krige）的碩士論文。克里金是南非威特沃特斯蘭德礁複合體距離加權平均金品位的開拓性繪圖者。克里金試圖根據少數鑽孔的樣本來估算金最可能的分布。英文動詞為to krige，最常用的名詞是kriging。在文獻中，該詞有時首字母大寫為Kriging。

雖然其基本公式的計算量很大，但克里金法可以使用各種高斯過程近似擴展到更大的問題。

## 主要原理

### 相關術語與技術

克里金通過計算給定點鄰域內函數已知值的加權平均來預測該點的函數值。該方法與迴歸分析密切相關。兩種理論都基於對協方差的假設推導出最佳線性無偏估計量，利用高斯-馬爾可夫定理證明估計值與誤差的獨立性，並使用非常相似的公式。即便如此，它們分別應用於不同的框架中：克里金法用於估計隨機場的一次實現，而回歸模型則基於多變量數據集的多次觀測。

克里金估計也可以被視為再生核希爾伯特空間中的樣條函數，其再生核由協方差函數給出。與經典克里金方法的區別在於解釋方式：樣條的動機是基於希爾伯特空間結構的最小範數插值，而克里金法的動機是基於隨機模型的期望平方預測誤差。

具有「多項式趨勢面」的克里金法在數學上等同於廣義最小二乘法多項式曲線擬合。

克里金法也可以被理解為貝葉斯優化的一種形式。克里金法從函數上的先驗分布開始。該先驗採用高斯過程的形式：來自函數的$N$個樣本將服從正態分布，其中任意兩個樣本之間的協方差是在兩點空間位置上評估的高斯過程協方差函數（或核）。隨後觀測到一組值，每個值都與一個空間位置相關聯。現在，通過將高斯先驗與每個觀測值的高斯似然函數相結合，可以在任何新的空間位置預測新值。產生的後驗概率分布也是高斯的，其均值和協方差可以簡單地從觀測值、它們的方差以及從先驗導出的核矩陣中計算出來。

### 地質統計估計量

<image>

在地質統計模型中，採樣數據被解釋為隨機過程的結果。這些模型在概念化中包含不確定性，這一事實並不意味著現象（如森林、含水層、礦床）是由隨機過程產生的，而是允許人們為未觀測位置的數量空間推論建立方法論基礎，並量化與估計量相關的不確定性。

在該模型背景下，隨機過程僅僅是處理從樣本中收集的數據集的一種方式。地質統計調製的第一步是創建一個最能描述觀測數據集的隨機過程。

來自位置$x_1$（地理坐標集的通用名稱）的值被解釋為隨機變量$Z(x_1)$的一次實現$z(x_1)$。在樣本分散的空間$A$中，存在相互關聯的隨機變量$Z(x_1), Z(x_2), \ldots, Z(x_N)$的$N$次實現。

這組隨機變量構成了隨機函數，其中只有一次實現是已知的，即觀測數據集$z(x_i)$。由於每個隨機變量只有一次實現，理論上不可能確定單個變量或函數的任何統計參數。地質統計形式主義中提出的解決方案包括「假設」隨機函數具有不同程度的「平穩性」，以便使某些統計值的推論成為可能。

例如，如果基於變量分布區域$A$中樣本的同質性，假設一階矩是平穩的（即所有隨機變量具有相同的均值），那麼就假設均值可以通過採樣值的算術平均值來估計。

與二階矩相關的平穩性假設定義如下：兩個隨機變量之間的相關性僅取決於它們之間的空間距離，而與它們的位置無關。因此，如果$\mathbf{h} = x_2 - x_1$且$h = |\mathbf{h}|$，則：$$C\big(Z(x_1), Z(x_2)\big) = C\big(Z(x_i), Z(x_i + \mathbf{h})\big) = C(h)$$ $$\gamma\big(Z(x_1), Z(x_2)\big) = \gamma\big(Z(x_i), Z(x_i + \mathbf{h})\big) = \gamma(h)$$為簡單起見，我們定義$C(x_i, x_j) = C\big(Z(x_i), Z(x_j)\big)$和$\gamma(x_i, x_j) = \gamma\big(Z(x_i), Z(x_j)\big)$。

該假設允許人們推論出這兩種度量——變異函數和協方差函數：$$\gamma(h) = \frac{1}{2|N(h)|} \sum_{(i,j)\in N(h)} \big(Z(x_i) - Z(x_j)\big)^2$$ $$C(h) = \frac{1}{|N(h)|} \sum_{(i,j)\in N(h)} \big(Z(x_i) - m(h)\big)\big(Z(x_j) - m(h)\big)$$其中：$$m(h)=\frac{1}{2|N(h)|} \sum_{(i,j)\in N(h)} \left( Z(x_i) + Z(x_j) \right)$$ $N(h)$表示觀測對$i$，$j$的集合，使得$|x_i - x_j| = h$，而$|N(h)|$是集合中對的數量。

在該集合中，$(i,;j)$和$(j,;i)$表示相同的元素。通常使用「近似距離」$h$，並結合一定的容差來實現。

### 線性估計

數量$Z \colon \mathbb{R}^n \to \mathbb{R}$在未觀測位置$x_0$的空間推論或估計，是根據觀測值$z_i = Z(x_i)$和權重$w_i(x_0)$；$i = 1, \ldots, N$的線性組合計算得出的：$$\hat{Z}(x_0) = <pre> \begin{bmatrix} w_1 & w_2 & \cdots & w_N \end{bmatrix} \begin{bmatrix} z_1 \\ z_2 \\ \vdots \\ z_N \end{bmatrix} = \sum_{i=1}^N w_i(x_0) Z(x_i).$$權重$w_i$旨在總結空間推論過程中的兩個極其重要的步驟：

- 反映樣本與估計位置$x_0$的結構「接近度」。
- 同時，它們應該具有去聚集效應，以避免由最終的樣本「聚類」引起的偏差。

在計算權重$w_i$時，地質統計形式主義有兩個目標：「無偏性」和「最小估計方差」。

如果將實際值$Z(x_0)$的點雲與估計值$\hat{Z}(x_0)$進行對比，全局無偏性、該場的「本徵平穩性」或廣義平穩性的標準意味著估計值的均值必須等於實際值的均值。

第二個標準指出平方偏差$\big(\hat{Z}(x) - Z(x)\big)$的均值必須最小，這意味著當估計值與實際值的點雲更加分散時，估計量就更加不精確。

## 方法

根據隨機場的隨機性質和假設的不同程度的平穩性，可以推導出計算權重的不同方法，即應用不同類型的克里金法。經典方法包括：

- 普通克里金法（Ordinary kriging）假設僅在$x_0$的搜索鄰域內存在恆定的未知均值。
- 簡單克里金法（Simple kriging）假設一階矩在整個域內平穩且均值已知：$E{Z(x)} = E{Z(x_0)} = m$，其中$m$是已知均值。
- 泛克里金法（Universal kriging）假設一個通用的多項式趨勢模型，例如線性趨勢模型$\textstyle E{Z(x)} = \sum_{k=0}^p \beta_k f_k(x)$。
- IRFk-克里金法假設$E{Z(x)}$為$x$的一個未知多項式。
- 指標克里金法（Indicator kriging）使用指示函數而不是過程本身，以便估計轉移機率。
- 多指標克里金法（Multiple-indicator kriging）是指標克里金法的一個版本，處理一系列指標。最初，MIK作為一種能更準確估算整體全球礦床濃度或品位的新方法，展現了相當大的前景。然而，由於所使用的區塊尺寸本質上很大，且缺乏採礦規模的解析度，這些優勢已被建模中其他固有的實用性問題所抵消。在這種情況下，條件模擬正迅速成為被接受的替代技術。
- 析取克里金法（Disjunctive kriging）是克里金法的非線性推廣。
- 對數正態克里金法通過對數對正數據進行插值。
- 潛在克里金法假設在非線性混合效應模型的潛在水平（第二階段）進行各種克里金處理，以產生空間功能預測。當分析空間功能數據${ (y_i, x_i, s_i) }{i=1}^n$時，該技術非常有用，其中$y_i = (y{i1} ,y_{i2}, \cdots, y_{iT_i})^\top$是跨越$T_i$周期的時間序列數據，$x_i = (x_{i1}, x_{i2}, \cdots, x_{ip})^\top$是包含$p$個協變量的向量，而$s_i = (s_{i1}, s_{i2})^\top$是第$i$個受試者的空間位置（經度、緯度）。
- 協同克里金法（Co-kriging）表示對來自多個來源且不同數據源之間存在關係的觀測數據進行聯合克里金處理。協同克里金法也可以在貝葉斯方法中實現。
- 貝葉斯克里金法（Bayesian kriging）偏離了對未知係數和超參數的優化，從貝葉斯角度來看，優化被理解為最大似然估計。相反，係數和超參數是從它們的期望值中估計出來的。貝葉斯克里金法的一個優點是，它允許量化克里金仿真器的證據和不確定性。如果使用仿真器來傳播不確定性，可以通過比較仿真器不確定性與總不確定性來評估克里金仿真器的質量（另見貝葉斯多項式混沌）。貝葉斯克里金法也可以與協同克里金法混合使用。

### 普通克里金法

未知值$Z (x_0)$被解釋為位於$x_0$的隨機變量，鄰近樣本的值$Z(x_i),\ i = 1, \ldots, N$也是如此。估計量$\hat{Z}(x_0)$也被解釋為位於$x_0$的隨機變量，是變量線性組合的結果。

克里金法旨在最小化在無偏條件下估計$Z(x_0)$的以下誤差的均方值：$$\epsilon(x_0) = \hat{Z}(x_0) - Z(x_0) = \begin{bmatrix} W^T & -1 \end{bmatrix} \cdot \begin{bmatrix} Z(x_1) & \cdots & Z(x_N) & Z(x_0) \end{bmatrix}^T = \sum^N_{i=1} w_i(x_0) \times Z(x_i) - Z(x_0)$$前面提到的兩個質量標準現在可以用新隨機變量$\epsilon(x_0)$的均值和方差來表示：

#### 無偏性

由於隨機函數是平穩的，$E[Z(x_i)] = E[Z(x_0)] = m$，權重之和必須為1，以確保模型是無偏的。可以看作如下：$$E[\epsilon(x_0)] = 0 \Leftrightarrow \sum^N_{i=1} w_i(x_0) \times E[Z(x_i)] - E[Z(x_0)] = 0$$ $$\Leftrightarrow m \sum^N_{i=1} w_i(x_0) - m = 0 \Leftrightarrow \sum^N_{i=1} w_i(x_0) = 1 \Leftrightarrow \mathbf{1}^T \cdot W = 1$$最小方差

兩個估計量可以都有$E[\epsilon(x_0)] = 0$，但圍繞其均值的分散程度決定了估計量質量的差異。為了找到具有最小方差的估計量，我們需要最小化$E[\epsilon(x_0)^2]$。$$\begin{align} <pre> \operatorname{Var}(\epsilon(x_0)) &= \operatorname{Var}\left(\begin{bmatrix} W^T & -1 \end{bmatrix} \cdot \begin{bmatrix} Z(x_1) & \cdots & Z(x_N) & Z(x_0) \end{bmatrix}^T\right) \\ &= \begin{bmatrix} W^T & -1 \end{bmatrix} \cdot \operatorname{Var}\left(\begin{bmatrix} Z(x_1) & \cdots & Z(x_N) & Z(x_0) \end{bmatrix}^T\right) \cdot \begin{bmatrix} W \\ -1 \end{bmatrix} \end{align}$$有關詳細解釋，請參閱協方差矩陣。$$\operatorname{Var}(\epsilon(x_0)) = \begin{bmatrix} W^T & -1 \end{bmatrix} \cdot \begin{bmatrix} \operatorname{Var}_{x_i} & \operatorname{Cov}_{x_ix_0} \\ \operatorname{Cov}_{x_i x_0}^T & \operatorname{Var}_{x_0} \end{bmatrix} \cdot \begin{bmatrix} W \\ -1 \end{bmatrix}$$其中字母$\left\{\operatorname{Var}_{x_i}, \operatorname{Var}_{x_0}, \operatorname{Cov}_{x_ix_0}\right\}$代表$\left\{\operatorname{Var}\left(\begin{bmatrix} Z(x_1) & \cdots & Z(x_N) \end{bmatrix}^T\right), \operatorname{Var}\big(Z(x_0)\big), \operatorname{Cov}\left(\begin{bmatrix} Z(x_1) & \cdots & Z(x_N) \end{bmatrix}^T, Z(x_0)\right)\right\}$

一旦定義了在$Z(x)$的所有分析領域中有效的協方差模型或變異函數（$C(\mathbf{h})$或$\gamma(\mathbf{h})$），我們就可以寫出任何估計量的估計方差表達式，它是樣本之間協方差以及樣本與待估計點之間協方差的函數：$$\begin{cases} <pre> \operatorname{Var}\big(\epsilon(x_0)\big) = W^T \cdot \operatorname{Var}_{x_i} \cdot W - \operatorname{Cov}_{x_ix_0}^T \cdot W - W^T \cdot \operatorname{Cov}_{x_ix_0} + \operatorname{Var}_{x_0}, \\ \operatorname{Var}\big(\epsilon(x_0)\big) = \operatorname{Cov}(0) + \sum_i \sum_j w_i w_j \operatorname{Cov}(x_i,x_j) - 2 \sum_iw_i C(x_i,x_0). \end{cases}$$從該表達式可以得出一些結論。估計方差：

- 一旦假設了均值和空間協方差（或變異函數）的平穩性，就無法量化任何線性估計量。
- 當樣本與待估計點之間的協方差減小時，估計方差增大。這意味著，當樣本距離$x_0$越遠，估計效果越差。
- 隨變量$Z(x)$的先驗方差$C(0)$增大而增大；當變量的分散程度較小時，區域$A$內任何點的方差都較低。
- 不取決於樣本的值，這意味著相同的空間配置（樣本與待估計點之間具有相同的幾何關係）在區域$A$的任何部分總是產生相同的估計方差；通過這種方式，方差無法衡量由局部變量產生的估計不確定性。

#### 方程組

$$W = \underset{\mathbf{1}^T \cdot W = 1}{\operatorname{arg,min}}\left( W^T \cdot \operatorname{Var}{x_i} \cdot W - \operatorname{Cov}{x_ix_0}^T \cdot W - W^T \cdot \operatorname{Cov}{x_ix_0} + \operatorname{Var}{x_0} \right)$$解決此優化問題（參見拉格朗日乘數）會產生「克里金系統」：$$\begin{bmatrix}\hat{W}\\\mu\end{bmatrix} = \begin{bmatrix} \operatorname{Var}_{x_i}& \mathbf{1}\\ \mathbf{1}^T& 0 \end{bmatrix}^{-1}\cdot \begin{bmatrix} \operatorname{Cov}_{x_ix_0}\\ 1\end{bmatrix} = \begin{bmatrix} \gamma(x_1,x_1) & \cdots & \gamma(x_1,x_n) &1 \\ \vdots & \ddots & \vdots & \vdots \\ \gamma(x_n,x_1) & \cdots & \gamma(x_n,x_n) & 1 \\ 1 &\cdots& 1 & 0 \end{bmatrix}^{-1} \begin{bmatrix}\gamma(x_1,x^*) \\ \vdots \\ \gamma(x_n,x^*) \\ 1\end{bmatrix}$$附加參數$\mu$是一個拉格朗日乘子，用於在最小化克里金誤差$\sigma_k^2(x)$時滿足無偏條件。====

### 簡單克里金法

<image>

簡單克里金法在數學上是最簡單的，但通用性最差。它假設隨機場期望值已知，並依賴於協方差函數。然而，在大多數應用中，期望值和協方差都無法預先得知。

應用「簡單克里金法」的實際假設是：

- 場的廣義平穩性（方差平穩）。
- 各處期望值為零：$\mu(x) = 0$。
- 已知的協方差函數$c(x, y) = \operatorname{Cov}\big(Z(x), Z(y)\big)$。

協方差函數是一個關鍵的設計選擇，因為它規定了高斯過程的屬性，從而決定了模型的行為。協方差函數編碼了諸如平滑度和周期性之類的資訊，這些資訊反映在生成的估計值中。一種非常常見的協方差函數是平方指數函數，它非常利於平滑函數估計。因此，在許多現實世界的應用中，它可能會產生較差的估計，特別是當真實的底層函數包含不連續性和快速變化時。

#### 方程組

「簡單克里金法」的克里金權重沒有無偏條件，由「簡單克里金方程組」給出：$$\begin{pmatrix} w_1 \\ \vdots \\ w_n \end{pmatrix} = <pre> \begin{pmatrix} c(x_1, x_1) & \cdots & c(x_1, x_n) \\ \vdots & \ddots & \vdots \\ c(x_n, x_1) & \cdots & c(x_n, x_n) \end{pmatrix}^{-1} \begin{pmatrix} c(x_1,x_0) \\ \vdots \\ c(x_n,x_0) \end{pmatrix}$$這類似於$Z(x_0)$對其他$z_1, \ldots, z_n$的線性回歸。

#### 估計

通過簡單克里金法進行的插值由下式給出：$$\operatorname{Var}\big(\hat{Z}(x_0) - Z(x_0)\big) = <pre> \underbrace{c(x_0,x_0)}_{\operatorname{Var}\big(Z(x_0)\big)} - \underbrace{\begin{pmatrix} c(x_1,x_0) \\ \vdots \\ c(x_n,x_0) \end{pmatrix}' \begin{pmatrix} c(x_1,x_1) & \cdots & c(x_1,x_n) \\ \vdots & \ddots & \vdots \\ c(x_n,x_1) & \cdots & c(x_n,x_n) \end{pmatrix}^{-1} \begin{pmatrix} c(x_1,x_0) \\ \vdots \\ c(x_n,x_0) \end{pmatrix}}_{\operatorname{Var}\big(\hat{Z}(x_0)\big)}$$克里金誤差由下式給出：$$\operatorname{Var}\big(\hat{Z}(x_0) - Z(x_0)\big) = \underbrace{c(x_0,x_0)}_{\operatorname{Var}\big(Z(x_0)\big)} - \underbrace{\begin{pmatrix} c(x_1,x_0) \\ \vdots \\ c(x_n,x_0) \end{pmatrix}' \begin{pmatrix} c(x_1,x_1) & \cdots & c(x_1,x_n) \\ \vdots & \ddots & \vdots \\ c(x_n,x_1) & \cdots & c(x_n,x_n) \end{pmatrix}^{-1} \begin{pmatrix} c(x_1,x_0) \\ \vdots \\ c(x_n,x_0) \end{pmatrix}}_{\operatorname{Var}\big(\hat{Z}(x_0)\big)},$$由此得到高斯-馬爾可夫定理的廣義最小二乘版本（Chiles & Delfiner 1999, p. 159）：$$\operatorname{Var}\big(Z(x_0)\big) = \operatorname{Var}\big(\hat{Z}(x_0)\big) + \operatorname{Var}\big(\hat{Z}(x_0) - Z(x_0)\big)$$貝葉斯克里金法

另見貝葉斯多項式混沌

### 屬性

- 克里金估計是無偏的：$E[\hat{Z}(x_i)] = E[Z(x_i)]$。
- 克里金估計尊重實際觀測值：$\hat{Z}(x_i) = Z(x_i)$（假設沒有發生測量誤差）。
- 如果假設成立，克里金估計$\hat{Z}(x)$是$Z(x)$的最佳線性無偏估計量。然而（例如Cressie 1993）：
- 與任何方法一樣，如果假設不成立，克里金法可能效果不佳。
- 可能存在更好的非線性及/或有偏方法。
- 當使用錯誤的變異函數時，無法保證任何屬性。然而，通常仍能獲得「良好」的插值。
- 「最佳」並不一定意味著「好」：例如，在沒有空間依賴性的情況下，克里金插值僅與算術平均值一樣好。
- 克里金法提供$\sigma_k^2$作為精確度的衡量標準。然而，該衡量標準依賴於變異函數的正確性。

## 應用

儘管克里金法最初是為地質統計學應用而開發的，但它是一種通用的統計插值方法，可以應用於滿足適當數學假設的隨機場採樣數據的任何學科。它可以用於已收集空間相關數據（二維或三維）的情況，並且需要估算實際測量值之間位置（空間間隙）的「填充」數據。

迄今為止，克里金法已用於多種學科，包括：

- 環境科學
- 水文地質學
- 採礦業
- 自然資源
- 遙感
- 房地產估價
- 集成電路分析與優化
- 微波器件建模
- 天文學
- 頁岩油井產油曲線預測

### 電腦實驗的設計與分析

在工程學中，另一個非常重要且快速增長的應用領域是對作為確定性電腦模擬響應變量的數據進行插值，例如有限元法（FEM）模擬。在這種情況下，克里金法被用作元建模工具，即在設計的電腦實驗集上建立的黑箱模型。在許多實際工程問題中，例如金屬材料成形過程的設計，單次FEM模擬可能長達數小時甚至數天。因此，設計並運行有限次數的電腦模擬，然後使用克里金插值器快速預測任何其他設計點的響應，效率會更高。因此，克里金法經常作為所謂的代理模型，實現在優化程序中。在混合整數輸入的情況下，也可以使用基於克里金法的代理模型。

## 相關

- 貝葉斯線性統計
- 高斯過程
- 多元插值
- 非參數回歸
- 徑向基函數插值
- 空間映射
- 空間依相關性
- 變異函數
- 梯度增強克里金法（GEK）
- 代理模型
- 資訊場論
- 反距離加權