interp2

内容

1 メニュー情報
2 概要
3 追加の情報
4 コマンドラインでの使用法
5 変数
6 説明
7 アルゴリズム
8 参考文献
9 関連のXファンクション

メニュー情報

XYからZを補間

概要

8つの手法を使用して、与えられたXYデータにおけるZデータを探すためにXYZデータの2D補間を実行します。

追加の情報

必要なOriginのバージョン:9.1SR0

コマンドラインでの使用法

1.interp2 ixy:=[XYZRandomGaus]"XYZ Random Gaussian"!(D,E) ir:=[XYZRandomGaus]"XYZ Random Gaussian"!(A,B,C) method:=rk pts:=3 smooth:=4;

変数

表示名	変数名	I/O とデータ型	デフォルト値	説明
補間するXY値	ixy	入力 XYRange	<unassigned>	与えられたXY値のために入力XY範囲を指定
入力	ir	入力 XYZRange	<active>	XYZデータのグループのための入力XYZ範囲を指定
手法	method	入力 int	0	使用するスムージングの手法を指定します。8つの方法が利用できます。オプションリスト nr:近傍 rk:ランダムKriging法 rc:ランダム Renka Cline rs:ランダム Shepard tps:ランダム TPS sp:スプライン ta:三角 wa:加重平均詳細は、アルゴリズムのセクションにあります。
検索半径	radius	入力 double	2	検索半径を指定します。手法として、ランダムKriging または、加重平均を指定した時に使用可能です。
最小ポイント	pts	入力 int	2	最小ポイント数を指定します。手法として、ランダムKriging を指定した時に使用可能です。
スムージング因子	smooth	入力 double	<unassigned>	Kriging またはスプラインのスムージン因子を指定します。手法として、ランダムKriging または、スプラインを指定した時に使用可能です。
重みの入力	weight	入力 vector		重みとしてベクトルを指定します。手法として、スプラインを指定した時に使用可能です。
重み関数局所因子	nw	入力 int	9	重み関数局所要因を指定します。手法として、ランダムShepard を指定した時に使用可能です。
二次補間局所因子	nq	入力 int	18	二次補間局所要因を指定します。手法として、ランダムShepard を指定した時に使用可能です。
補間の結果	oz	出力 vector	<new>	補間の結果を出力する場所を指定します。

説明

このXファンクションは、与えられた(x,y)の値に対するzの値を検索するために、既存のXYZデータのグループを補間するために使用され、、近傍、ランダムKriging、ランダムRenka Cline、ランダムShepard、ランダムTPS、スプライン、三角、加重平均の8つの手法が利用可能です。

xyからzを補間するために、 interp2 Xファンクションを使用します。

スクリプトウィンドウで、以下を実行します。
```
interp2 -d;
	
```
すると interp2 ダイアログボックスが開きます。
補間に使用するXYの値として、ワークシートから(x,y)データを選択し、補間のための入力としてxyzデータの範囲を選択します。
手法ドロップダウンリストから手法を選択します。
補間の結果として、目的の列を選択します。

アルゴリズム

Interp2 は、周囲の位置における観測値から未計測の位置における変数の推定を提供します。XファンクションInterp2では、8つの手法を利用でき、下表に各手法の文書とアルゴリズムがあります。

手法	アルゴリズム
近傍	最近傍補間
ランダム Kriging法	通常 Kriging 補間
ランダム Renka Cline法	三角法ベースの Renka Cline 補間
ランダム Shepard法	修正Shepard補間の手法
ランダム TPS法	ランダム TPS補間
スプライン	双三次B-スプライン補間
三角	三角法ベースの Renka Cline 補間
加重平均法	逆距離加重法補間

(ランダム)通常Kriging補間

Originには、関数ocmath_2d_kriging_scatがあり、通常Kriging補間を実行します。通常Krigingの基本的な手順は、以下のようになります。

$x_1, x_2, ..., x_n$ は観測されたランダムポイントとし、 $z(x_1), z(x_2), ..., z(x_n)$ は一連のポイントそれぞれの計測された値とします。計測されていない与えられたポイント $x_0$ は、通常Krigingに従い、この未計測ポイントの推定値 $z^*(x_0)$ は、以下のように記述できる計測されたポイントの加重線形組み合わせによって生成されます。

$z^*(x_0) = \sum_{i=1}^nw_iz(x_i)$

ここで $w_i$ は対応するポイントの重みです。

未計測ポイント $x_0$ の推定値 $z^*(x_0)$ を取得するために、計算に寄与するそれぞれの計測されたポイントの重みを理解する必要があります。

第一に、ポイントのためのランダム関数は、結果の普遍性のために変動しないと仮定します。つまり、下式で予測された値はゼロである必要があります。

$E[z^*(x_0)-z(x_0)] = E[\sum_{i=1}^nw_iz(x_i)-z(x_0)] = \sum_{i=1}^nw_iE[z(x_i)]-E[z(x_0)] =0$

$E[z(x_i)] = E[z(x_0)] = m$ （ $m$ は平均）ですので、 $\sum_{i=1}^nw_i = 1$ を取得します。

次に、エラーのための最適の分散を満たす必要があります。これは、以下の式で表される分散の最小化です。

$\sigma_0^2 = Var[z^*(x_0)-z(x_0)] = E[((z^*(x_0)-z(x_0))-E[z^*(x_0)-z(x_0)])^2] = E[(z^*(x_0)-z(x_0))^2] = E[(\sum_{i=1}^nw_iz(x_i)-z(x_0))^2] = E[(z(x_0))^2] - 2E[z(x_0)\sum_{i=1}^nw_iz(x_i)] + E[(\sum_{i=1}^nw_iz(x_i))^2] = C(x_0, x_0) - 2\sum_{i=1}^nw_iC(x_i, x_0) + \sum_{i=1}^n\sum_{j=1}^nw_iw_jC(x_i, x_j) = minimum$

上述の式による最小化は、Lagrange(言語)パラメータ $\mu$ があり、以下のようにします。

$F = \sigma_0^2 - 2\mu (\sum_{i=1}^nw_i - 1)$

$w_i$ と $\mu$ に関する上式の偏第一微分を計算し、取得します。

$\left\{\begin{matrix}\frac{\partial F}{\partial w_i} = 2\sum_{j=1}^nw_jC(x_i, x_j) - 2C(x_i, x_0) - 2\mu = 0 \\ \frac{\partial F}{\partial \mu} = -2(\sum_{i=1}^nw_i-1) = 0 \end{matrix}\right.$

単純化されたものを取得します。

$\left\{\begin{matrix}\sum_{j=1}^nw_jC(x_i, x_j) - \mu = C(x_i, x_0) \\ \sum_{i=1}^nw_i = 1 \end{matrix}\right.$

そして $\sigma_0^2 = C(x_0, x_0) - \sum_{i=1}^nw_iC(x_i, x_0) + \mu$ となります

次のようにします。

$C = \begin{bmatrix}C(x_1, x_1) & C(x_1, x_2) & \cdots & C(x_1, x_n) & 1 \\ C(x_2, x_1) & C(x_2, x_2) & \cdots & C(x_2, x_n) & 1 \\ \vdots & \vdots & \: & \vdots & \vdots \\ C(x_n, x_1) & C(x_n, x_2) & \cdots & C(x_n, x_n) & 1 \\ 1 & 1 & \cdots & 1 & 0 \end{bmatrix}, W = \begin{bmatrix}w_1 \\ w_2 \\ \vdots \\ w_n \\ -\mu\end{bmatrix}, D = \begin{bmatrix}C(x_1, x_0) \\ C(x_2, x_0) \\ \vdots \\ C(x_n, x_0) \\ 1\end{bmatrix}$

偏第一微分は拡張でき、行列形式 $C\cdot W = D$ として再記述できます。

重みを解くには、行列形式の両側に $C$ 行列の逆上列である $C^{-1}$ をかけ、 $W = C^{-1}\cdot D$ を取得します。

最後のステップは、共分散の決定です。上述のすべての計算は、ランダム関数内のランダム変数はすべて同じ平均と分散を持つという仮定をベースにしています。このような仮定は、モデルバリオグラムとモデルの共分散間の関係を、以下のように導きます。

$\gamma (x_i, x_j) = \frac{1}{2}E[(z(x_i)-z(x_j))^2] = \frac{1}{2}E[(z(x_i))^2] + \frac{1}{2}E[(z(x_j))^2] - E[z(x_i) \cdot z(x_j)] = E[z^2]-E[z(x_i) \cdot z(x_j)] = (E[z^2] - m^2)-(E[z(x_i) \cdot z(x_j)] - m^2) = \sigma^2-C(x_i, x_j)$

この関係式を偏第一微分に置くと、他の形式(行列形式は同様)を取得できます。

$\left\{\begin{matrix}\sum_{j=1}^nw_j\gamma (x_i, x_j) + \mu =\gamma (x_i, x_0) \\ \sum_{i=1}^nw_i = 1 \end{matrix}\right.$

Originでは、パワーモデルがバリオグラムの推定に使用され、

$\gamma (h) = \left\{\begin{matrix}0 \;\;\;\;\;\;\;\;\;\;\;\;\;\; |h| = 0 \\ a + b|h|^\lambda \;\;\;\; |h| > 0 \end{matrix}\right.$

ここで $a = 0, b = 1$ と $0 < \lambda < 2$ になります。 $h$ は $x_i$ と $x_j$ 間の距離を参照し、 $\gamma (h) = \gamma (x_i, x_j)$ となります。

この手法が使用された時の、このXファンクションのダイアログでは、ocmath_2d_kriging_scat 関数の $\lambda$ はスムージングファクター、検索半径と最小ポイント(補間を含む最小のポイント) は、radius と noctMin パラメーターに対応します。

ランダム Renka Cline法

ランダム Renka Cline法のアルゴリズムはは、NAGにより提供された関数nag_2d_scat_interpolant (e01sac) をベースにしています。しかし、NAGライブラリによると、Advice on Replacement Calls for Withdrawn/Superseded Functions, the nag_2d_scat_interpolant (e01sac) は nag_2d_shep_interp (e01sgc) または、nag_2d_triang_interp (e01sjc)で置き換えられます。つまり、ランダム Renka Cline法は、 nag_2d_triang_interp をベースにした三角法に似たアルゴリズムをもちますが、後半は、最小の最大補間誤差が生成される可能性があるので推奨されます。

ランダム Shepard法

詳細はランダム Shepardを参照してください。

ランダム TPS法

詳細はランダム TPSを参照してください。

スプライン

nag_2d_spline_fit_scat (e02ddc) をベースにした手法で、散布データのセットに3次スプライン近似を計算します。

三角

nag_2d_triang_interp (e01sjc)のアルゴリズムをベースにた手法で、 Renka Clineによる手法を使用した $r =1,2,\cdots ,m$ のための $m$ 散布データポイント $(x_r,y_r,f_r)$ のセットを通した補間曲面 $F(x,y)$ で構築されます。構築された曲面は連続した曲線で、連続した一次導関数を持ちます。

加重平均

加重平均は、ocmath_gridding_weighted_averageをベースにした手法です。未知の値 $P(x_p, y_p, z_p)$ は、指定された検索半径 r 内の真の値 $Q_i$ の周囲から計算されます。補間値は次式となります。

$Z_p=\frac{\sum_{i=1}^{n}\frac{Z_i}{3d_i^2+1}}{\sum_{i=1}^{n}\frac{1}{3d_i^2+1}}$

$d_i$ は $P$ と $Q_i$ の間の距離です。検索半径内に値がない場合、Zp は欠損値になります。検索半径を拡大するということは、各ポイントは隣り合うポイントとの内部的な関係を強め、詳細な部分を失って、より滑らかな曲面を作成するということになります。

比較

8つの手法を比較してみます。6つの関数F1からF6（指数ピーク、急斜面の曲面、山のピーク、緩やかなピーク、急勾配のピーク、球形曲面）に、8つの手法を適用しました。手法は、Ref. 2由来のものです。1089ポイントの補間は、100ノードの一様分布をベースにしています。平均誤差と最大誤差は、それぞれ1089ポイントの実際の値と補間値間の差の平均と最大です。下図は、平均誤差と最大誤差の棒グラフです。

関数式：

$F1: 0.75e^{-\frac{(9x-2)^2+(9y-2)^2}{4}}+0.75e^{\frac{-(9x+1)^2}{49}-\frac{9y+1}{10}}+0.5e^{-\frac{(9x-7)^2+(9y-3)^2}{4}}-0.2e^{-(9x-4)^2-(9y-7)^2}$

$F2: \frac{tanh(9y-9x)+1}{9}$

$F3: \frac{1.25+cos(5.4y)}{6+6(3x-1)^2}$

$F4: e^{\frac{-(81/16)((x-0.5)^2+(y-0.5)^2)}{3}}$

$F5: e^{\frac{-(81/4)((x-0.5)^2+(y-0.5)^2)}{3}}$

$F6: \frac{\sqrt{64-81((x-0.5)^2+(y-0.5)^2}}{9}-0.5$

参考文献

Isaaks and Srivastava.An Introduction to Applied Geostatistics, Chapter 12.Oxford University Press (1989)
R. J. Renka and A. K. Cline.A Triangle-Based $C^1$ Interpolation Method, Rocky Mountain Journal of Mathematics.14, 223-237 (1984)
P. Dierckx.An Algorithm for Surface-Fitting with Spline Functions, IMA Journal of Numerical Analysis.1, 267-283 (1981)