Introduction à l'Interpolation par Krigeage
L'interpolation par krigeage est une technique géostatistique avancée utilisée pour estimer des valeurs à des emplacements non échantillonnés, basée sur des observations à des points connus. Contrairemetn aux méthodes d'interpolation plus simples, le krigeage prend en compte la corrélation spatiale des données, modélisée par un variogramme. Ce document explore différentes approches pour implémenter le krigeage en utilisant MATLAB.
- Krigeage avec la Boîte à Outils de Cartographie (Mapping Toolbox)
La Boîte à outils de Cartographie de MATLAB offre des fonctions spécifiques pour le calcul de variogrammes et l'interpolation par krigeage, notmament variogram et gkrige. Cette approche est robuste et intégrée pour les applications géospatiales.
% 1. Préparation des données d'observation
% Génération de points d'observation synthétiques avec des valeurs associées.
% Ces données simulent, par exemple, des relevés de température ou d'altitude.
latitude_obs = [10 12 15 18 20 22 25 28 30 32];
longitude_obs = [5 7 10 12 15 17 20 22 25 27];
valeurs_observées = [5.2 6.1 7.5 8.9 10.3 9.8 8.5 7.0 6.5 5.8] + randn(1, length(latitude_obs)) * 0.5;
% Ajout d'un léger bruit pour simuler des données réelles.
% 2. Calcul du variogramme expérimental
% La fonction 'variogram' analyse la structure de corrélation spatiale des données.
% Elle retourne un objet variogramme contenant le variogramme expérimental.
[obj_variogramme, distances_lag, semi_variance_exp] = variogram(latitude_obs, longitude_obs, valeurs_observées);
% 3. Définition de la grille cible pour l'interpolation
% Création d'une grille fine de points où les valeurs seront estimées.
pas_grille = 0.5;
[grille_long_interp, grille_lat_interp] = meshgrid(min(longitude_obs):pas_grille:max(longitude_obs), ...
min(latitude_obs):pas_grille:max(latitude_obs));
% 4. Exécution de l'interpolation par krigeage
% La fonction 'gkrige' utilise le variogramme calculé pour effectuer l'interpolation.
% Elle retourne la grille des valeurs interpolées.
[~, ~, ~, valeurs_interpolées_krigeage] = gkrige(latitude_obs, longitude_obs, valeurs_observées, ...
obj_variogramme, grille_lat_interp, grille_long_interp);
% 5. Visualisation des résultats du krigeage
figure;
surf(grille_long_interp, grille_lat_interp, valeurs_interpolées_krigeage, 'EdgeColor', 'none');
hold on;
% Affichage des points d'observation d'origine.
plot3(longitude_obs, latitude_obs, valeurs_observées, 'o', 'MarkerSize', 7, 'MarkerFaceColor', 'red', 'DisplayName', 'Observations');
title('Interpolation par Krigeage (avec Mapping Toolbox)');
xlabel('Longitude');
ylabel('Latitude');
zlabel('Valeur Estimée');
colormap jet; % Choix d'une palette de couleurs.
colorbar;
view(3); % Vue en 3D
grid on;
axis tight;
legend('show');
- Approche Conceptuelle du Krigeage Ordinaire (Implémentation Simplifiée)
Pour comprendre les mécanismes sous-jacents du krigeage, il est utile d'examiner une implémentation simplifiée du krigeage ordinaire. Cette section démontre comment un modèle de variogramme peut être défini et utilisé pour calculer les poids d'interpolation pour chaque point de la grille.
% 1. Données d'entrée pour l'exemple conceptuel
% Coordonnées et valeurs de points échantillonnés.
x_donnees = [1 2 3 4 5 6 7 8 9 10];
y_donnees = [1 4 2 5 3 6 4 7 5 8];
z_donnees = [10 12 15 13 16 14 17 15 18 16] + randn(1, length(x_donnees)) * 0.3; % Valeurs avec bruit
% 2. Définition de la grille d'interpolation
% Points de la grille où les valeurs doivent être prédites.
pas_grille_manuel = 0.2;
[grille_X_pred, grille_Y_pred] = meshgrid(min(x_donnees):pas_grille_manuel:max(x_donnees), ...
min(y_donnees):pas_grille_manuel:max(y_donnees));
points_a_predire = [grille_X_pred(:), grille_Y_pred(:)];
valeurs_predites_manuel = zeros(size(points_a_predire, 1), 1);
% 3. Définition du modèle de variogramme (exemple: modèle Sphérique)
% Le modèle sphérique est une fonction qui décrit la semi-variance en fonction de la distance (h).
% Paramètres: C0 (effet de pépite), C (palier - pépite), A (portée).
effet_pepite = 0.1; % Nugget
contribution_spatiale = 5; % Sill - Nugget
portee_modele = 3; % Range
% Fonction de variogramme sphérique
modele_variogramme_sph = @(h) (effet_pepite + contribution_spatiale * (1.5 * (h / portee_modele) - 0.5 * (h / portee_modele).^3)) .* (h <= portee_modele) + ...
(effet_pepite + contribution_spatiale) * (h > portee_modele);
% 4. Fonction utilitaire pour calculer la distance euclidienne
calculer_distance = @(p1, p2) sqrt(sum((p1 - p2).^2));
% 5. Boucle d'interpolation par Krigeage Ordinaire
% Pour chaque point de la grille, on calcule les poids de krigeage.
num_observations = length(x_donnees);
points_observes = [x_donnees(:), y_donnees(:)];
for idx_pred = 1:size(points_a_predire, 1)
point_actuel_pred = points_a_predire(idx_pred, :);
% Calcul des distances entre le point à prédire et tous les points observés
distances_h_obs_pred = arrayfun(@(j) calculer_distance(point_actuel_pred, points_observes(j,:)), 1:num_observations);
% Calcul des valeurs du variogramme pour ces distances
gamma_obs_pred = modele_variogramme_sph(distances_h_obs_pred)';
% Construction de la matrice de variogramme K (entre points observés)
K_variogramme = zeros(num_observations, num_observations);
for r = 1:num_observations
for c = 1:num_observations
dist_entre_obs = calculer_distance(points_observes(r,:), points_observes(c,:));
K_variogramme(r, c) = modele_variogramme_sph(dist_entre_obs);
end
end
% Ajout des contraintes pour le krigeage ordinaire (somme des poids = 1)
% On construit un système linéaire augmenté (Kriging Matrix)
matrice_krigeage_étendue = zeros(num_observations + 1, num_observations + 1);
matrice_krigeage_étendue(1:num_observations, 1:num_observations) = K_variogramme;
matrice_krigeage_étendue(1:num_observations, num_observations + 1) = 1;
matrice_krigeage_étendue(num_observations + 1, 1:num_observations) = 1;
matrice_krigeage_étendue(num_observations + 1, num_observations + 1) = 0; % Pour le multiplicateur de Lagrange
% Vecteur du second membre pour la résolution (k)
vecteur_second_membre = [gamma_obs_pred; 1];
% Résolution du système linéaire pour obtenir les poids de krigeage (lambda)
% Utilisation de la pseudo-inverse pour plus de robustesse si la matrice est quasi-singulière.
poids_lambda_avec_multiplicateur = pinv(matrice_krigeage_étendue) * vecteur_second_membre;
% Les poids pour l'interpolation sont les premiers 'num_observations' éléments
poids_krigeage = poids_lambda_avec_multiplicateur(1:num_observations);
% Calcul de la valeur prédite au point actuel
valeurs_predites_manuel(idx_pred) = sum(poids_krigeage .* z_donnees');
end
% 6. Visualisation des résultats de l'implémentation conceptuelle
figure;
surf(grille_X_pred, grille_Y_pred, reshape(valeurs_predites_manuel, size(grille_X_pred)), 'EdgeColor', 'none');
hold on;
% Affichage des points d'observation d'origine.
plot3(x_donnees, y_donnees, z_donnees, 'o', 'MarkerSize', 7, 'MarkerFaceColor', 'blue', 'DisplayName', 'Observations');
title('Krigeage Ordinaire - Implémentation Conceptuelle');
xlabel('Coordonnée X');
ylabel('Coordonnée Y');
zlabel('Valeur Estimée');
colormap parula;
colorbar;
view(3);
grid on;
axis tight;
legend('show');