Matlab K-means聚類算法對多光譜遙感圖像進行分類

作者： 白藝亭

測試了下matlab自帶kmeans函數，作者編寫函數，以及ENVI下的Kmeans方法，對比其效果，代碼及結果圖展示見下。(K均值聚類的matlab代碼，對其算法進行實現。)

1. K-Means無監督聚類算法：

ENVI中，在主菜單上，選擇Classification>>Unsupervised>>K-Means,在Classification Input File中選擇TM影像，單擊OK。打開K-means Parameters對話框中，設置以下的參數：

1)分類數量(Number of Classes)：一般為最終輸出分類數量的2-3倍。

2)變換閾值(Change Threshold):5。每當一類的變化像元數小于閾值時，結束迭代過程。這個值越小得到的結果越精確，運算量也越大。

3)最大迭代次數(Maximum Iterations)：30(我選的).迭代次數越大，精度越高。

4)距離類別的值的最大誤差(Maximum Stdev From Mean)。此數可選。

5)允許的最大距離誤差：(Maximum Distance Error);可選

6)選擇路徑，OK執行。

以上參考東騰的博客。

2. 自編matlab函數實現

2.1被調function

function [new_class_label] = Kmeans_of_muldim(data,k,change_threshold,iteration)

% 功能：實現多光譜遙感數據非監督分類算法之K-means聚類算法

%Author: Mr. BAI

% 輸入：data是s*fl*b的矩陣，s為列數(sample)，fl為行數(fileline)，b為波段數(band)；

% k 為類別數,如果有背景值，背景值會歸到某一地類中去，到時再用矢量邊界圖形裁剪一下即可。我考慮過將出現次數最多

% 的背景值單獨劃歸一類，但是程序設計時不好判斷，取數組中元素出現次數最多的像元為一類，有點大膽，因為無背

% 景的圖像像元值也可能出現這種情況；

% change_threshold變化閾值，ENVI中默認為0.05；

% iteration為最大迭代次數，ENVI中默認為1

% 輸出：new_class_label為聚類后的矩陣，賦予每個行列號一個類別標簽，之后可在GIS或者ENVI中出圖

% Reference：https://www.cnblogs.com/dongteng/p/5415071.html

[fl,s,b] = size(data);

%original_seed為迭代前的種子，存放一個k行，b個波段數值列的數組

old_seed = zeros(k,b);

%newseed為迭代后的新種子，存放一個k行，b個波段數值列的數組

new_seed = zeros(k,b);

%-------------------------------------------------------------------------------------------------

% 產生k個隨機種子作為遙感圖像各地物類別的種子像元

%-------------------------------------------------------------------------------------------------

index_record = zeros(1,k);

for i = 1:k

index_i = round(rand()*fl*s);

judge = find(index_record == index_i);

%如果已經有這個值了，那么重新循環取值

if isempty(judge) == 0

i = i-1;

continue;

end

index_record(i) = index_i;

%計算取到的隨機值對應圖像的行列號

fl_index = floor(index_i/s);%行號

sample_index = index_i - fl_index*s;%列號

%將該種子像元的b個波段值存入

old_seed(i,:) = data(fl_index,sample_index,:);

end

%--------------------------------------------------------------------------------------

% 下面進行迭代，如果本次分別所有類新得到的像元數目變化在change_threshold內，則認為分類完畢。

%--------------------------------------------------------------------------------------

n = 1;

new_class_label = zeros(fl,s);

while n

distance_matrix = zeros(fl,s,k);

for kind = 1:k

sum = 0;

for i=1:b

temp = power(abs(data(:,:,i)-old_seed(kind,i)),2);

sum = sum+temp;

end

%每個像元與初始7個類別中心的歐式距離

ou_distance = sqrt(sum);

distance_matrix(:,:,kind) = ou_distance;

end

%給給各類別賦值類別標注

for i=1:fl

for j=1:s

currentpixel_vector = distance_matrix(i,j,:);

currentpixel_class = find(currentpixel_vector == min(currentpixel_vector));

new_class_label(i,j) = currentpixel_class(1);

end

end

%計算新的各類別中心

for i=1:k

id = find(new_class_label==i);

for j=1:b

temp1 = data(:,:,j);

temp2 = temp1(id);

new_seed(i,j)= mean(temp2(:));

end

end

new_class_pixcel_number = zeros(1,k);

for i=1:k

new_class_pixcel_number(i) = length(find(new_class_label(:)==i));

end

%Change threshold:0.05

if n == 1

old_class_pixcel_number = ones(1,k);

end

%size(new_class_pixcel_number)

if max(abs((new_class_pixcel_number-old_class_pixcel_number)./old_class_pixcel_number)) < change_threshold || n>iteration

break;

end

n=n+1;

if max(abs((new_class_pixcel_number-old_class_pixcel_number)./old_class_pixcel_number)) >change_threshold

%old_class_label = new_class_label;

old_class_pixcel_number = new_class_pixcel_number;

old_seed = new_seed;

continue;

end

end

end

2.2 main函數調用上述function

clc;

clear;

t0 = cputime;

cd 'E:\MATLAB\'

data=imread('nantong_city_landsat8.tif');%讀取純數據

[multi_data,r]=geotiffread('nantong_city_landsat8.tif'); % read the geo information

info=geotiffinfo('nantong_city_landsat8.tif'); % read the geo information

class_result = Kmeans_of_muldim(data,5,0.05,30);

geotiffwrite('K-means_class.tif',class_result,r,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);

t1 = cputime;

during = t1 - t0;

disp('耗時：');

disp(during);

耗時：

139.59375

時間長是因為程序中Kmeans_of_muldim()函數中使用了逐像元循環，用了139秒，后面附上改進的算法， MKmeans_of_muldim()，Matlab K-means聚類算法改進對多光譜遙感圖像進行分類(二)。

3 matlab庫中的kmeans()

clc;

clear;

%matlab自帶的kmeans使用方法:

t0 = cputime;

cd 'E:\MATLAB\'

data=imread('nantong_city_landsat8.tif');%讀取純數據

[multi_data,r]=geotiffread('nantong_city_landsat8.tif'); % read the geo information

info=geotiffinfo('nantong_city_landsat8.tif'); % read the geo information

[fl,s,b] = size(data);

dat = zeros(fl*s,b);%matlab K-means算法要求輸入矩陣是一個列向量組成的矩陣，列數為波段數，每一列為fl*s的像元值

for i=1:b

dat(:,i) = reshape(data(:,:,i),fl*s,1);

end

class_result = kmeans(dat,5);只能傳兩參

out_data = reshape(class_result,fl,s);

geotiffwrite('K-means_class_matlab.tif',out_data,r,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);

t1 = cputime;

during = t1 - t0;

disp('耗時：');

disp(during);

耗時：

38.796875

4 結果展示

4.1ENVI操作結果展示：

0.05,30.

4.2 Kmeans_of_muldim,作者的matlab代碼結果展示：

0.05,30

4.3K-means matlab自帶函數結果展示：

效果相當，下面對Kmeans_of_muldim進行改進。

見下一篇：

Matlab K-means聚類算法改進對多光譜遙感圖像進行分類(二)

后來，我們都學會了勇敢，不遺憾，不留戀。417看完了劉若英的線上演唱會，1500w人同時在線的，很壯觀。

與50位技術專家面對面20年技術見證，附贈技術全景圖

總結

以上是生活随笔為你收集整理的光谱分类算法 matlab,Matlab K-means聚类算法对多光谱遥感图像进行分类（一）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。