クリギングに関する備忘録

知ったかぶりを減らしたいんや

クリギングってアレでしょ、セミバミ…バリ…オグラムとかいうので空間自己相関を考慮みたいな、それで内挿をいい感じにやっちゃうみたいな雰囲気のやつでしょ、という程度のゆるふわな理解のまま生きてきてしまったので、いちおう導出ぐらい追っておこうと思いまして。

と思ってネットを漁るとわかったようなわからないような概説はたくさん出てくるのですが、なんでもっと必要十分に書いてくれないの？ねえなんで？という気持ちが迸ったのでまとめておきます。

なお詳細は文献[1][2]などに詳しいですので、もっと真面目に背景を追いたい方は、そちらをご覧ください。

セミバリオグラム

空間座標 x を変数にもつ関数 z(x) があるとします。まあ標高マップみたいなものをイメージしてください。画素が x, 対応する標高が z(x) です。

このとき、x から距離 h だけ離れた点での関数値 z (x+h) が、z (x) からどれだけ異なっているかの指標として、誤差二乗の期待値を使うことにします。とりあえず方向性はないものとし、1次元で考えておきましょう。

この期待値を「バリオグラム」といい、その半分である γ(h) を「セミバリオグラム」と言います[3]。バリオグラムあるいはセミバリオグラムは、もとの点から離れるごとに値が平均的にどのようにばらついていくかを表す関数です。

γ(h) はどんな形になりそうでしょうか。

• 空間的に近い点同士は、値も似ていそう（空間自己相関）

• なんならh→0では、もとの値に一致するはず

• ある程度以上空間的に離れてしまうと、もとの点の情報との関連はなくなりそう（空間自己相関が消え、γはマップ全体の分散に一致する）

などを考慮すると、h→0 で γ(h) →0, そこから単調増加し、ある程度のところで上限（飽和）値をとる、ということで

などのモデル（の実数倍）でよく表現できそうです。

ただし、場合によっては h が0に近づいてもγ(0)が0にならず値をもつ（ナゲット効果）ケースもあるため、フィッティングモデルは一般に b0 + b1×g(h) という形が使われます。切片である γ(0) = b0をナゲット、飽和に達する h の値をレンジ、飽和値 γ(∞) をシルなどと呼ぶようです。

セミバリオグラムと空間分散・共分散との関係

あとで使うので見ておきます。

ただし、z の空間平均（期待値）を μ、空間分散をVar[z(x)]、h 離れた点同士の空間共分散をCov[z(x), z(x+h)]と書いています。なお任意点から h 離れた点の座標 x+h を新たな変数としても、原点が変わるだけなので空間平均と空間分散は同じとすれば E[z(x)] = E[z(x+h)] = μ、Var[z(x)] = Var[z(x+h)]と考えてよいでしょう。

h 離れた点同士の空間共分散 Cov [z(x), z(x+h)] はけっきょく h の関数なので C(h)と簡素に書くことにし、定義式から分散 Var[z(x)] は C(h) の h→0 での値 C(0) とも書けることに注意すれば

となることが分かります。上式からγ(0) = C(0) - C(0) = 0 が要請されますが、前述のとおり、γ(0) が0にならずオフセットを持つような場合も考慮し

としておくことにしましょう（γ(0)はオフセット補正量 i.e., ナゲット）。

でクリギングの話はいつ出てくるんや

今からです。お待たせしました。関数値が未知の点 x0 があり、このまわり（実際上はレンジ程度の範囲）に関数値が既知の点 xi (i = 1,2, ..., n) が分布しているとします。このとき、既知点の関数値 z(xi) の重みづけ平均によって未知の値 z(x0) を推定することにします。つまり以下のような内挿問題です。

ハットは推定量を表します。wi は重みで、w1+w2+...+wn = 1を満たすものとします（不偏条件）。近傍データからの線形補間ならTIN、多項式補間ならスプライン、距離の逆数ならIDWですが、今回は以下のように、推定値と真値の二乗誤差の期待値（空間平均）が最小となる条件を求めることにします。

左辺を少々展開しますと、

ただし xj, xk 間の距離を hjk としました。C(h) は先述の通り距離 h 離れた点同士の空間共分散で、C(0) は単なる空間分散に一致します。

ここで「セミバリオグラムと空間分散・共分散との関係」を用いれば、

とも書けます。この式が、w1 + w2 + ... + wn = 1 の制約条件のもと最小化される条件を求めればよいです。ラグランジュの未定乗数法により、以下の関数 L(w, λ) の w, λについての偏微分が0となる条件を求めます[4]。

右辺第三項の wi (i = 1, 2, ... , n)での微分は、wiが順列で2回登場することに注意して

Lのλでの微分は制約条件そのものになります。以上を行列にまとめれば、

となり、γに関する行列の逆行列をもとめて両辺左からかければ重みwiが求まります。念のためですが、γ(hjk) は点 j, k の距離 hjk におけるセミバリオグラム値です。

ちなみにセミバリオグラム γ(h) を用いた式と共分散 C(h) を用いた式は、ちょうど全体の符号が反転した書き方になり、極値をとる条件としては結局同じになるので、C(h) を使って重みを求めることもできます。結果の γ という文字を Cという文字に書き換えればよいだけです。

これで自信をもって「クリギングってアレでしょ、セミバリオグラムみたいなやつで以下略」と言えるようになりました。めでたいです。

ちなみにここで述べたのは「普通クリギング（ordinary kriging）」と呼ばれるタイプで、いくつかの変種もあるようです。たとえば、普通クリギングだけだと異方性を考慮できませんが、近傍の空間範囲にトレンド成分を追記するような「普遍クリギング（universal kriging）」もあります。特徴量として z(xi) だけでなく別の f(xi) も考慮し、^z(x0)=Σ wi z(xi) + vi f(xi) とおくこともでき、この場合は「コクリギング（cokriging）」と言います。また、z(x0)が連続値ではなくカテゴリ量 (0 or 1) であるような場合は「指示クリギング（indicator kriging）」というようです。

[1] 正路, 小池. 2007, 講座「地球統計学」ヴァリオグラム―データの空間的連続性の解析. 日本地熱学会誌, 29(3), 125-140.

[2] 正路, 小池. 2007, 講座「地球統計学」クリギング: 誤差を考慮した空間データの補間. 日本地熱学会誌, 29(4), 183-194.

[3] わざわざ用語を分けるほどのものなのかとも思いますが、もしかして深遠な理由があるのかもしれません。。

[4] ラグランジュの未定乗数法自体については、PCAの回で概要を説明しています。