Geoscience Reference
In-Depth Information
Yg = reshape(Xg2,[],1);
We then allocate memory for the kriging estimates
Zg
and the kriging
variance
s2_k
by:
Zg = Xg * NaN;
s2_k = Xg * NaN;
We now krige the unknown values at each grid point:
for k = 1 : length(Xg)
DOR = ((x - Xg(k)).^2 + (y - Yg(k)).^2).^0.5;
G_R = (nugget + sill*(1 - exp(-3*DOR/range))).*(DOR>0);
G_R(n+1) = 1;
E = G_inv * G_R;
Zg(k) = sum(E(1:n,1).*z);
s2_k(k) = sum(E(1:n,1).*G_R(1:n,1))+E(n+1,1);
end
h e i rst command computes the distance between the grid points
(Xg,Yg)
and the observation points
(x,y)
. We then create the right-hand-side
vector of the kriging system using the variogram model
G_R
and adding
one to the last row. We next obtain the matrix
E
with the weights and the
Lagrange multiplier. h e estimate
Zg
at each point
k
is the weighted sum
of the observations
z
. Finally, the kriging variance
s2_k
of the grid point is
computed and we can plot the results. We i rst create a grid of the kriging
estimate and the kriging variance:
r = length(R);
Z = reshape(Zg,r,r);
SK = reshape(s2_k,r,r);
A subplot on the let presents the kriged values:
subplot(1,2,1)
h = pcolor(Xg1,Xg2,Z);
set(h,'LineStyle','none')
axis equal
ylim([0 200])
title('Kriging Estimate')
xlabel('x-Coordinates')
ylabel('y-Coordinates')
colormap(jet)
colorbar
h e let hand subplot presents the kriging variance:
subplot(1,2,2)
h = pcolor(Xg1,Xg2,SK);
set(h,'LineStyle','none')
axis equal
