Discussion :

[Cuve9]

Bonjour,

Je souhaite coder un programme sur du Krigeage, il s'agit d'une méthode d'interpolation.
J'ai écrit tout mon algo, lorsque je lance la fonction, il n'y a pas d'erreur mais scilab me renvoie des valeurs bizarres... Je pense que l'erreur pourrais provenir de la résolution de mon système linéaire par la méthode de remontée (après avoir procédé à une décomposition QR), ou peut être de la construction des matrices A et A1...

Voici mon code:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
function Fp_krigeage=krigeage(Xp,X,F)   
 
//Estimation de F(Xp) par la méthode de krigeage, où Xp est un vecteur contenant les coordonnées du point dont on souhaite déterminer la valeur de la fonction F. 
//X est une matrice contenant les coordonnées des points dont on connait déjà la valeur de la fonction F, donc la matrice X possède 2 ligne (la première correspondant à l'abscisses des points et la deuxième correspondant aux ordonnées de nos points). 
//Enfin F est un vecteur (une seule ligne) contenant la valeur de la fonctions F pour les différents points contenus dans la matrice X
 
//Construction des vecteurs dist et n contenant respectivement les distances entre chaque paires de points et les nombres de paires de 
//points séparés par la même distance.
 
dist=[] //Initialisation de 'dist' par un vecteur vide. 'dist' contient toutes les valeurs différentes des distances h_k entre deux points
        //Xi(xi,yi) et Xj(xj,yj)
n=[]    
compteur=0
 
    for i=1:(size(X,2)-1) //Détermination du vecteur n qui contient, pour chaque indice k, le nombre de paires de points séparés par une
                          //même distance h_k.
        for j=i+1:size(X,2)
            d=sqrt((X(1,i)-X(1,j))^2+(X(2,i)-X(2,j))^2)
            trouve=%F;
            for k=1:(size(dist,2)) //'dist' et 'n' sont des vecteurs lignes
                if d==dist(1,k) then
                    n(1,k)=n(1,k)+1
                    trouve=%T;
                    break;
                end
            end
            if ~trouve then
                compteur=compteur+1
                dist(1,compteur)=d
                n(1,compteur)=1
            end
        end
    end
 
 
//Construction du vecteur semi-variance contenant les semi-variances de chaque distance différente séparant une paire de points
 
semivar=[]
 
    for k=1:size(dist,2) 
        S=0
        for l=1:n(1,k)
            S=S+(X(1,l)-X(2,l))^2
        end
        semivar(k,1)=(0.5/n(1,k))*S //'semivar' est un vecteur colonne
    end
 
 
//Construction de la matrice A contenant les semi-variances de h_ij
 
m=size(X,2)
semivar_h0=(0.5/m)*S*ones(m,1) //Semi-variance pour une distance h nulle (n(h)=m car h=0 pour m paires de points : les m paires Xi_Xi)
A1=zeros(m,m)+diag(semivar_h0) //Initialement, la matrice A1 ne contient que semi-var_h0 sur la diagonale
 
    for i=1:size(X,2)
        for j=1:size(X,2) 
             if (i~=j) then
                 h=sqrt((X(1,i)-X(1,j))^2+(X(2,i)-X(2,j))^2)
                     for k=1:size(dist,2)
                           if h==dist(1,k) then
                                A1(i,j)=semivar(k,1) //Construction de la matrice A1 contenant les semi-variances de chaque paires de pts
                            end
                     end
              end
         end
    end
A=A1 //Construction de la matrice A finale
A(:,m+1)=1
 
A(m+1,:)=1
A(m+1,m+1)=0
 
 
//Interpolation linéaire
 
S_x=0
S_y=0
 
    for i=1:size(X,2)
        S_x=S_x+X(1,i)
        S_y=S_y+X(2,i)
    end
 
x_moy=S_x/size(X,2)
y_moy=S_y/size(X,2)
 
S1=0
S2=0
 
    for i=1,size(X,2)
        S1=S1+X(i,1)*X(i,2)
        S2=S2+X(i,1)*X(i,1)
    end
 
S1_interpo=size(X,2)*S1
S2_interpo=size(X,2)*S2
 
a_interpol=(S1_interpo-x_moy*y_moy)/(S2_interpo-x_moy*x_moy)
 
b_interpol=y_moy-a_interpol*x_moy
 
 
//Construction du vecteur B contenant les semi-variances de h_ip
 
semivarbis=[] 
h_p=[]
 
    for i=1:m
        h_p(1,i)=sqrt((X(1,i)-Xp(1))^2+(X(2,i)-Xp(2))^2) //h_p est un vecteur ligne
        semivarbis(i,1)=a_interpol*h_p(1,i)+b_interpol //semivarbis est un vecteur colonne 
    end
 
B=semivarbis
B(m+1,1)=1
 
 
//Résolution du sytème A.W=B par la méthode QR
 
[Q,R]=qr(A)
Qt=zeros(length(B),length(B))
 
    for i=1:length(B)
       for j=1:length(B)
            Qt(i,j)=Q(j,i)   //  Construction de la matrice Q transposée. Ceci est nécessaire par la suite pour calculer le produit Qt*B 
                             //lors de la résolution du système
        end
    end
 
Somme=0
    for i=1:length(B)
        for j=1:length(B)
            Somme=Somme+Qt(i,j)*B(j) // Construction du vecteur E=Q*B. On a donc le système : R*W=E
        end
        E(i)=Somme
    end
 
       //méthode de remontée résoudre le système
 
W=zeros(size(X,2)) // Construction du vecteur W contenant les poids. La dimension de ce vecteur est égale au nombre de points contenus
                   //dans X
 
W(size(X,2))=E(size(X,2))/R(size(X,2),size(X,2)) 
 
   for i=size(X,2)-1:-1:1
       sum=0
        for j=i+1:size(X,2)
            sum=sum+R(i,j)*W(j)
        end
        W(i)=(E(i)-sum)/R(i,i)
    end
 
 
 
//Calcul de Fp_kriegage
 
Fp_krigeage=0
 
    for i=1:size(X,2)
        Fp_krigeage=Fp_krigeage+W(i)*F(i)
    end
 
endfunction

Merci d'avance pour votre aide,

Bien cordialement.

[__dardanos__]

Salut,
Le calcul des semi-variances me semble bizarre.
Je ne comprends pas la différence entre X[2,:] et F[:].