% Fuzzy Logic Determination of Lithoogy from Well Log Data: Application To

function[]=flc(test)

clc, close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Surveys in Geophysics

% Fuzzy logic determination of lithoogy from well log data: application to

% the KTB project dataset (Germany).

% David Bosch, Juanjo Ledo, Pilar Queralt.

% Departament de Geodinàmica i Geofísica, Universitat de Barcelona.

% C/ Martí i Franquès, s/n. 08028 Barcelona, Spain.

% email:

% last revision date: 29 January 2013

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Note: required fields are delimited with < >. Delete these symbols before

% running the program.

% Input parameters

% The "test" input parameter is a number for executing the program on

% 1: the training well dataset; 2: the test well dataset and 3: the

% synthetic test dataset (if available).

% Membership function creation (histogram normalization):

nbins=<30>; % Number of bins to be computed in the histograms

lits=<1:9>; % Lithologies to be considered

a=dlmread(<training well dataset filename>);

% The training well dataset file has to be organized as follows: the

% first column is depth, the following columns contain the property

% values (one column for each log property) and the last column

% contains the reference lithology number (a previous numbering is

% required).

props=1:8; % The vector of properties to be used in classification.

for prop=props

maxprop=max(a(:,prop+1));

minprop=min(a(:,prop+1));

x0=(maxprop-minprop)/nbins;

x(prop,:)=minprop+x0/2:x0:maxprop-x0/2;

for j=lits

lit{j}=a(a(:,10)==j,1:9);

c{j,prop}=hist(lit{j}(:,prop+1),x(prop,:));

mf{j,prop}=c{j,prop}/max(c{j,prop});

end

end

% End of membership function creation

%Reading of the test dataset file

if test==2

a=dlmread(<test well dataset filename>);

elseif test==3

a=dlmread(<synthetic well dataset filename>);

end

% Start of the classification process

clas=zeros(<9>,<10>); % Classification matrix variable. 9 rows (for nine

% lithologies) and 10 columns (for 10 lithologies plus the Unknown

% class).

for i=1:length(a(:,1)) % For each data point...

for j=lits % ... and for each considered lithology...

for prop=props % ... and property..

% ... each degree of membership (DoM) is obtained by

% interpolation:

if a(i,prop+1)>x(prop,end)

k=length(x(prop,:));

else

k=find(x(prop,:)>=a(i,prop+1),1,'First');

end

if k>1 & k<length(x(prop,:))

if mf{j,prop}(k)-mf{j,prop}(k-1)>0

DoM(j,prop)=mf{j,prop}(k-1)+(a(i,prop+1)-x(prop,k-1))*(mf{j,prop}(k)-mf{j,prop}(k-1))/(x(prop,k)-x(prop,k-1));

else

DoM(j,prop)=mf{j,prop}(k)+(x(prop,k)-a(i,prop+1))*(mf{j,prop}(k-1)-mf{j,prop}(k))/(x(prop,k)-x(prop,k-1));

end

elseif k==1

if a(i,prop+1)>=x(prop,1)

DoM(j,prop)=mf{j,prop}(1)*a(i,prop+1)/(0.5*(x(prop,2)-x(prop,1))); %Es pot canviar, per una extrapolacio pex.

else

DoM(j,prop)=0;

end

else

if a(i,prop+1)<=x(prop,end)

DoM(j,prop)=mf{j,prop}(end)*(x(prop,end)+0.5*(x(prop,2)-x(prop,1))-a(i,prop+1))/(0.5*(x(prop,2)-x(prop,1)));

else

DoM(j,prop)=0;

end

end

end

% Then for each property the combination operator is applied

% to obtain the final lithology DoMs:

% Arithmetic mean

% litDoM(j)=mean(DoM(j,props));

% Geometric mean

% [aa,bb]=size(props);

% litDoM(j)=(prod(DoM(j,props)))^(1/bb);

% Harmonic mean

% litDoM(j)=8/(sum(1./DoM(j,props)));

% Normalized Euclidean norm

% [aa,bb]=size(props);

% litDoM(j)=norm(DoM(j,props))/sqrt(bb);

% Fuzzy AND

% litDoM(j)=min(DoM(j,props));

% Fuzzy OR

% litDoM(j)=max(DoM(j,props));

% Fuzzy Algebraic Sum

% litDoM(j)=1-prod(1-DoM(j,props));

% Fuzzy Gamma

gamma=<0.9>;

litDoM(j)=((1-prod(1-DoM(j,props)))^gamma)*(prod(DoM(j,props)))^(1-gamma);

end

% Lithology assignation:

if sum(litDoM==max(litDoM))>1

% There has been a tie between two or more lithology DoMs,

% so the sample is assigned to the "Unknown" class.

clas(a(i,end),length(clas))=clas(a(i,end),length(clas))+1;

else

% There is only one maximum lithology DoM.

clas(a(i,end),find(minim==max(minim)))=clas(a(i,end),find(minim==max(minim)))+1;

end

% Simplified three-class system:

SlitDoM=[max(litDoM(1:3)) max(litDoM(4:6)) max(litDoM(7:9))];

class(i,:)=ceil(max(SlitDoM)/3); % Simplified classification

DoMs(i,:)=SlitDoM; % Retaining the simplified lithology DoMs.

model(i,:)=ceil(a(i,end)/3); % Simplified reference column

end

%Displaying the classification matrix

disp(' ')

disp(num2str(clas))

disp(' ')

% Writing the simplified model, classification and lithology DoMs:

dlmwrite(['model' num2str(test) '.txt'],model,'delimiter','\t','precision',10)

dlmwrite(['class' num2str(test) '.txt'],class,'delimiter','\t','precision',10)

dlmwrite(['DoMs' num2str(test) '.txt'],DoMs,'delimiter','\t','precision',10)

clear all

close all

% End of the program