Surfer のクリギングにおけるバリオグラムのモデル化 ── チュートリアル

このチュートリアルは、Randal Barnes 氏が Golden Software のために執筆したオリジナルの「Variogram Tutorial」を改訂したものです。本改訂版は、Surfer 17 以降の［Grid Data］ダイアログにあるバリオグラムモデリングのページの画面と用語に準拠しています。

1. はじめに

バリオグラムは、データセットの空間的な連続性、つまり粗さの度合いを特徴づけるものです。2つのデータセットの通常の1次元統計量（平均・標準偏差などの記述統計量）がほぼ同じでも、空間的な連続性はまったく異なる場合があります。バリオグラムが必要となる理由の一端は、第2節で説明します。

バリオグラム解析は、データから計算する実験バリオグラムと、それに当てはめるバリオグラムモデルとで構成されます。実験バリオグラムは、指定した離隔距離と方向にあるすべての観測値の組について、z 値の差の2乗の1/2を平均して計算し、2次元のグラフとして描きます。実験バリオグラムの計算に用いる数式の詳細は、第3節を参照してください。

バリオグラムモデルは、空間的な関係を記述する数学関数（線形・指数・球形など、候補になる関数）の中から選びます。実験バリオグラムの曲線と、候補となる各関数の曲線とを見比べ、形がよく合う関数をモデルとして選びます。

各関数の曲線の形を示すグラフは、Surfer ヘルプ（Surfer User’s Guide）のトピック「Variogram Model Graphics」を参照してください。幾何異方性（方向によって空間的な連続性が異なること）を考慮するには、データセット内の方向ごとに、実験バリオグラムとモデルバリオグラムを別々に計算できます。

2. バリオグラムは何を表すのか

2つの人工的なデータセットを考え、それぞれ A、B と呼ぶことにします。両データセットの代表的な記述統計量を Table 1.1 に示します。

	A	B
Count (データ件数)	15251	15251
Average (平均値)	100.00	100.00
Standard Deviation (標準偏差)	20.00	20.00
Median (中央値)	100.35	100.92
10 Percentile	73.89	73.95
90 Percentile	125.61	124.72

Table 1.1　2つの例題データセットの代表的な記述統計量（Count＝データ数、Average＝平均、Standard Deviation＝標準偏差、Median ＝中央値、10/90 Percentile＝10/90 パーセンタイル）

両データセットのヒストグラムを Figure 1.1 と Figure 1.2 に示します。これらを見る限り、2つのデータセットはほぼ同じです。

しかし、この2つのデータセットには、代表的な記述統計量やヒストグラムでは捉えられない大きな違いがあります。対応する等高線図（Figure 1.3、Figure 1.4）を比べると分かるように、データセット A はデータセット B よりも粗くなっています。ただし、A が B より「ばらつきが大きい」とは言えません。2つのデータセットは標準偏差が同じで、高い値・低い値の大きさも同じだからです。目に見える違いは、ばらつきではなく、きめ（テクスチャ）の違いです。

具体的には、データセット A は B よりも空間的に急激に変化しています。連続する高い値の領域（赤い部分）と低い値の領域（青い部分）は、平均すると A のほうが B よりも小さくなっています。このような違いは、サンプリング計画、サイトの特性評価、さらには空間予測全般を大きく左右することがあります。

代表的な記述統計量やヒストグラムが、この2つの例題データセットのきめの違いを定量化するどころか、見分けることすらできないのは当然です。これらの計算には、データの空間的な位置が組み込まれていないからです。

バリオグラムも定量的な記述統計量です。ただし、グラフに描くと、データセットの空間的な連続性（すなわち粗さの度合い）を特徴づけることができます。2つのデータセットのバリオグラムを Figure 1.5 と Figure 1.6 に示します。曲線の立ち上がりの傾きの違いは明らかです。

3. バリオグラムとは何か

バリオグラムの数学的な定義は次のとおりです。

ここで z(x,y) は位置 (x,y) における対象変数の値、ε は統計的期待値の演算子です。バリオグラム γ() は点どうしの離隔 (Δx,Δy) の関数であり、特定の位置 (x,y) の関数ではないことに注意してください。この数学的定義は有用な抽象化ですが、実際の観測値にそのまま当てはめるのは容易ではありません。

n 個の観測データ {(x₁,y₁,z₁), (x₂,y₂,z₂), …, (x_n,y_n,z_n)} を考えます。ここで (x_i, y_i) は観測値 i の位置、z_i はそれに対応する観測値です。観測値の組は重複なしで n(n-1)/2 通りあり、それぞれの組について次の離隔ベクトルを計算できます。

(Δx_i,j, Δy_i,j) = (x_i-x_j, y_i-y_j)

ある特定の離隔ベクトル (Δx,Δy) についてバリオグラムを推定したいときは、離隔ベクトルがこの対象の離隔にほぼ等しいデータの組をすべて使います。

(Δx_i,j,Δy_i,j)≈(Δx,Δy)

そのような組すべての集合を S(Δx,Δy) とします。

S(Δx,Δy)={(i,j)|(Δx_i,j,Δy_i,j)≈(Δx,Δy)}

さらに、S(Δx,Δy) に含まれる組の数を N(Δx,Δy) とします。観測データからバリオグラムを推定するには、次の実験バリオグラムの式を使います。

つまり、ある特定の離隔ベクトルに対する実験バリオグラムは、ほぼそのベクトルだけ離れたすべての観測値の組について、z 値の差の2乗の1/2を平均して計算します。

4. バリオグラムグリッド

観測データが n 個あれば、観測値の組は重複なしで n(n-1)/2 通りあります。そのため、中程度の大きさのデータセットでも膨大な数の組が生じます。たとえば n = 500 なら n(n-1)/2 = 124,750 組です〔訳注：原文の 124,745 は計算の誤り〕。これほど多くの組を処理するのは、高速なコンピューターでも時間がかかることがあります。Surfer はすべての組をあらかじめ計算し、必要な和と差をバリオグラムグリッドに格納します（注：バリオグラムグリッドは、マップの作成に使うグリッドとは形式が異なります）。クリギングで使うバリオグラムとバリオグラムグリッドを表示するには、次のようにします。

1. Home | Grid Data | Grid Data をクリックします。
2. グリッド化手法（Gridding Method）として Kriging を選びます。
3. Dataset 1 欄の横にある Browse ボタンをクリックし、Surfer の Samples フォルダーにあるデータファイル ExampleDataSetC.xls を選びます。
4. X、Y、Z の列として、それぞれ A 列、B 列、C 列を指定します。
5. 必要に応じて Filter Data ボタンをクリックし、重複データ（Duplicates）とデータ除外（Data Exclusion）のフィルタリングを指定します。
6. View Data ボタンと Statistics ボタンで、データを確認することもできます。
7. Next をクリックして、バリオグラムを表示します。

ダイアログの上部にバリオグラムのグラフが表示されます。点の付いた赤い線が全方向の実験バリオグラム、青い線が当てはめたバリオグラムモデルの最初の試案です（当てはまりはよくありませんが）。実験バリオグラムの線上の数字は、各ラグに含まれるデータの組の数です。

バリオグラムサーフェス（ダイアログ下部の円形の図）は、あらゆる方向と距離についてバリオグラム関数の値を示します〔訳注：2点の組は逆向きに見ても値の差の2乗が変わらないため、(Δx,Δy) と (−Δx,−Δy) の値は等しく、図は中心について点対称になる〕。中心は2点間の離隔距離が 0 の位置で、半径が最大ラグ距離（ラグサイズ × ラグ数）です。図の各セルはラグサイズに等しい正方形で、推定値に応じて色分けされています。組が1つもないセルは空白のままです。マウスポインターをバリオグラムサーフェスの上に置くと、そのセルのバリオグラム関数の値と、その計算に使われた組の数がツールチップに表示されます。バリオグラムのグラフは、このサーフェスから「扇形」の部分を取り出して平均したデータを示しています。扇形の向きと角度の幅は、右側のプロパティリストの Direction と Tolerance で制御します。これらの設定は、図の上のコントロールをドラッグして対話的に操作することもできます。

5. 全方向バリオグラムのモデル化

既定では、上記で紹介した最初のグラフは全方向バリオグラムです（Tolerance は 90 度）。適切なモデル（成分）の種類、部分シル、ナゲット効果は、全方向バリオグラムをもとに選びます。

5.1 バリオグラムモデルの種類を選ぶ

バリオグラムモデルは無数に考えられます。Surfer では、用意されている11種類の成分を組み合わせて、何千通りものバリオグラムモデルを構成できます。ナゲット効果と組み合わせれば、ほとんどのデータセットには、線形モデル、指数モデル、球形モデルの3つのいずれかで十分です。この3つのモデルの例を Figure 5.1 に示します。

実験バリオグラムがいつまでも横ばいにならない場合は、通常、線形モデルが適しています。横ばいにはなるものの、そこまでずっと「曲線的」に上がっていく場合は、指数モデルを検討します。最初は直線的に上がり、その後急に曲がって横ばいになる場合は、まず球形モデルを選ぶのがよいでしょう。

ExampleDataSetC.xls のデータでは、曲線が増加したあと横ばいになっているように見えるため、球形モデルが適していそうです（指数モデルを試してもかまいません）。Variogram Component #1 セクションで、Component type を Linear から Spherical に変更します。

5.2 バリオグラムモデルのスケールと長さのパラメーターを選ぶ

次に、球形モデルの Partial sill（部分シル）と Range（レンジ）を、反復的な方法（つまり試行錯誤）で設定します。Partial sill は、バリオグラムが横ばいになる y 軸上の高さです。グラフを見ると 400〜450 が妥当そうなので、425 と入力します。球形モデルの Range は、バリオグラムが横ばいになるラグ距離です。これもグラフから 30〜40 が妥当そうなので、35 と入力します。グラフはすぐに更新されます。

最初の推測としては悪くありませんが、描き直された曲線をよく見ると、モデル（青い線）が実験バリオグラム（赤い線と点）の右側にあることから、Range が少し長すぎるようです。Range を 30 に設定し直します。これでもまだ少し長いので、Range を 29 にしてみます。これでバリオグラムによく当てはまります。

5.3 バリオグラムのナゲット効果を選ぶ

実験バリオグラムが縦軸上で 0 でない切片をもつように見える場合は、モデルにナゲット効果の成分が必要かもしれません。［Grid Data］のバリオグラムモデリングのページにある Display Statistics のレポートに含まれる NN Gamma Z の値は、ほとんどの場合、ナゲット効果の定量的な上限の目安になります。

Surfer では、ナゲット効果を Error variance（誤差分散）と Micro variance（微小分散）の2つの副成分に分けます。どちらも負でない値で、2つの和が見かけ上の 0 でない切片に等しくなるようにします。Error variance は観測の再現性を表し、サンプリングの誤差と分析（測定）の誤差の両方を含みます。Error variance は、重複サンプル間の差の分散を計算して決めるのが最善です。Micro variance は、一般的なサンプル間隔よりも短い離隔距離における未知のバリオグラムの代わりとなるものです。これは、実験バリオグラムの見かけ上の 0 でない切片から誤差分散を差し引いて決めるのが最善です。

この例のモデルは縦軸と 0 で交わるように見えるので、ナゲット効果は適用しません。

6. バリオグラムの異方性のモデル化

実験バリオグラムは、方向によってレンジの長さが異なることがよくあります。これを幾何異方性（geometric anisotropy）と呼びます。幾何異方性は、線形バリオグラムモデルでは方向による傾きの違いとして、球形モデルでは方向による見かけの Range の違いとして現れます。バリオグラムの異方性を調べるには、Experimental セクションのプロパティを編集します。全方向の実験バリオグラムは、すべての方向の挙動を平均したものです。特定の方向を見るには、Direction と Tolerance のパラメーターを変更します。まず Tolerance を 30 度に設定します。バリオグラムサーフェスの図に、0 度（このモデルでは東）を中心とし、角度の許容幅が ±30 度の「扇形」が表示されることに注目してください。

ダイアログボックスを横にドラッグして、バリオグラムのグラフが見えるようにします。すべての方向をすばやく見渡すには、バリオグラムサーフェスの図で扇形の中心線（0° の線）をクリックし、360° 一周ドラッグします。これを何回か繰り返してください。方向を変えるにつれて、グラフ上の実験バリオグラムの形と組の数が更新されることに注目してください。

90°〜180° の方向では、全方向モデルが実験バリオグラムの上昇部分よりも一貫して右側にあるように見えます。逆に 0°〜90° の方向では、全方向モデルが実験バリオグラムの上昇部分よりも一貫して左側にあります。これは幾何異方性があることをはっきり示しています。

バリオグラムサーフェスの図で Direction のコントロールをドラッグし、実験バリオグラムの見かけの Range と Partial sill が最大になる方向（たとえば 30°）に合わせます。Variogram component #1 セクションで、この方向別バリオグラムによく当てはまるように Range パラメーターを変更します。40 では長すぎ、30 では短すぎ、37 でよく当てはまるようです。Anisotropy Angle を 30° に変更します（この時点ではグラフは変わりません）。

Experimental セクションで Direction を 120° に変更します（これは、先に求めた最大長の方向 30° に直交する方向です。120° = 30° + 90°）。再び Variogram component #1 セクションで Anisotropy Ratio を変更し、現在の方向別実験バリオグラムにモデルを当てはめます。Ratio の値は 1.6 がよく当てはまるようです。

モデルを最終的に目で確かめるため、「一回転」させてみましょう。バリオグラムサーフェスの図で Direction のコントロールを 360° 一周ドラッグします。実験バリオグラムの形の変化に、バリオグラムモデルが追従することに注目してください。モデルのパラメーターを微調整して実験バリオグラムにさらに合わせることもできますが、過剰にモデル化しないでください。期待できるのは有効数字1〜2桁までです。

適切な初期パラメーターを選んだら、AutoFit ボタンをクリックしてモデルを微調整します。OK をクリックして既定の設定を受け入れると、AutoFit により Anisotropy Ratio 1.742、Angle 38.56° が得られます。

7. 実務上の目安

データをよく知ること。実験バリオグラムを計算する前に、空間を考慮しない通常の統計量を計算しておきます。バリオグラムモデリングのページにある Display Statistics ボタンを使うと、データの最小値、最大値、中央値、平均、標準偏差、分散、歪度を表示できます。Surfer でポストマップや分類ポストマップ（classed post map）を作成すれば、散布図として表示できます。ヒストグラムや累積度数グラフの作成には Grapher を使います。

• 過剰にモデル化しないこと。実験バリオグラムの重要な特徴を再現できる最も単純なモデルが、最良のモデルです。
• 迷ったら、グリッド化には既定のバリオグラムモデルを使います。単純な線形バリオグラムモデルでも、たいていは許容できるグリッドが得られます。初期のデータ解析では特にそうです。ただし、既定のバリオグラムで生成したクリギング標準偏差グリッドには意味がないことに注意してください〔訳注：クリギングの推定値はバリオグラムを定数倍しても変わらないが、標準偏差はバリオグラムの大きさに応じて変わる。既定のバリオグラムは大きさを実験バリオグラムに合わせていないため、標準偏差の値に根拠がない〕。
• 明確で疑問の余地のない物理的な根拠がない限り、3対1 を超える異方性比は使わないでください。実験バリオグラムが 3対1 を超える異方性を示しているように見えても、そのような強い異方性に明確な物理的根拠がない場合は、データにトレンドがある可能性があります。バリオグラム解析の前に、データのトレンド除去を検討してください。
• モデルを試し、得られたグリッドを等高線図で確認します。候補となるバリオグラムモデルが複数ある場合は、それぞれからグリッドと等高線図を作成します。大きな違いがなければ、最も単純なバリオグラムモデルを選びます。
• バリオグラムのレンジは、Z 値の空間的な変動に現れる物理的な異常域の平均的な大きさに近いことがよくあります。信頼できる実験バリオグラムがない場合は、この目安を使ってバリオグラムのレンジを仮定できます。
• 一定の値の周りで変動する実験バリオグラムは、異常なバリオグラムではありません。一般的なサンプル間隔のスケールでは Z 値どうしに相関がない、ということを示しています。このような場合、どのグリッド化手法を使っても等高線図はデータを確かに表すものにはならず、局所的な特性を詳しく把握するには、より短いサンプル間隔でデータを追加する必要があります。
• 次の3つの条件を満たす場合、標本分散はバリオグラムのシルの妥当な近似値になります。
  1. データが対象領域全体に均等に分布している（ポストマップで確認できます）。
  2. 対象領域全体にわたって、データに大きなトレンドがない。
  3. 対象領域の大きさが、バリオグラムの有効レンジの3倍を超えている。

記事原文

• Variogram modelling for kriging in Surfer – a tutorial