为减少信息冗余,同时为了对数据进行降维,采用最佳指数因子(Optimum Index Factor,OIF)建立最优波段特征组合。其基本原理是:图像中所涵盖的信息量与其标准差成正比,标准差越大,信息量就越多;图像的独立性与波段间的相关系数成反比,其相关系数越低,信息冗余度越小,其独立性越好。此方法综合了各波段间的关联性及单波段图像的信息量,得到了广泛应用。
可以参考以下论文:
Chavez P S, GL B, LB S. Statistical method for selecting Landsat MSS ratios[J]. 1982.
(最早由Chavez提出的OIF)
赵庆展,刘伟,尹小君,张天毅.基于无人机多光谱影像特征的最佳波段组合研究[J].农业机械学报,2016,47(03):242-248+291.
裴欢,孙天娇,王晓妍.基于Landsat 8 OLI影像纹理特征的面向对象土地利用/覆盖分类[J].农业工程学报,2018,34(02):248-255.
OIF计算公式
原始论文中这样写
通常写成这种形式
选择3 个波段组合并计算其OIF 值,OIF 值越大,波段标准差越小,波段间相关性系数越小,说明波段组合数据质量最优。
matlab程序计算OIF
1)因为我既用了原始影像计算OIF,也用了不同波段之间的比值计算OIF,所以代码开头会有两个部分,但其实参与OIF运算的主变量是X。而且我把波段比值中出现的异常值(主要是0,其实还有负值没处理)赋成了NaN,如果不处理0则计算OIF的结果会出现NaN。部分函数中omitnan是忽略了NaN值。
2)我的matlab版本较低,计算标准差是自己写的过程,因为是多维数组,所以没法直接用std计算,据官网文档显示2018b以上版本可以直接计算多维数组所有元素的标准差。
而且自己写的计算标准差公式为
这里用1/(n-1)考虑到无偏估计量的概念,详情参考数理统计部分。
3)我的原始影像为int16型,在这里转为了double类型,目的是为了后续计算波段比值,如果用int16型计算比值结果不对
clear all
clc
Image = imread('20180316QFLY_GS.tif'); % 读取影像,为int16类型,如果直接把异常值换为NaN则换不了,转成double类型
Image_D = double(Image);
% 波段比值
Image_D(Image_D == 0) = NaN;
XRatio_Combination = nchoosek(1:size(Image_D,3),2); % 比值的组合
X_Combination = nchoosek(1:size(XRatio_Combination,1),3); % 从那么多比值组合里选三个
X = zeros(size(Image,1),size(Image,2),size(XRatio_Combination,1));
for m = 1:size(XRatio_Combination,1)
% 计算波段比值并按照XRatio_Combination的顺序存储到X
X(:,:,m) = Image_D(:,:,XRatio_Combination(m,1)) ./ Image_D(:,:,XRatio_Combination(m,2));
end
% 原始影像
% X = Image_D; % 原始波段计算OIF
% X_Combination = nchoosek(1:size(X,3),3); % n选3,n是波段数,这个用于原始波段计算OIF
%%
% 计算标准差和相关系数
X_StdDev = zeros(1,size(X,3));
% 不含NaN可以直接用std2
% for i = 1:size(X,3)
% X_StdDev(i) = std2((X(:,:,i))); % 各波段标准差
% end
for i = 1:size(X,3)
X_Mean1 = mean(mean(X(:,:,i),'omitnan'));
X_Sq = (X(:,:,i) - X_Mean1).^2;
numNaN = length(find(isnan(X))); % 统计X中的NaN个数,在后续计算标准差时去除这一部分
X_StdDev(i) = sqrt(sum(sum(X_Sq),'omitnan')/(size(X,1) * size(X,2) - numNaN - 1)); %这里未将NaN统计到总数中
end
XCom_StdDev = [X_Combination,zeros(size(X_Combination,1),1)];
for i = 1:size(X_Combination,1) % 注意这里的循环结尾值是组合数,也就是X_Combination的行数
XCom_StdDev(i,4) = X_StdDev(X_Combination(i,1)) + X_StdDev(X_Combination(i,2)) + ...
X_StdDev(X_Combination(i,3)); % 计算组合波段的StdDev之和
end
XCom_Corr = [X_Combination,zeros(size(X_Combination,1),1)];
for j = 1:size(X_Combination,1)
% corrcoef的'Rows','complete'用于忽略NaN值
c1 = corrcoef(double(X(:,:,X_Combination(j,1))),double(X(:,:,X_Combination(j,2))),'Rows','complete');
c2 = corrcoef(double(X(:,:,X_Combination(j,1))),double(X(:,:,X_Combination(j,3))),'Rows','complete');
c3 = corrcoef(double(X(:,:,X_Combination(j,2))),double(X(:,:,X_Combination(j,3))),'Rows','complete');
XCom_Corr(j,4) = c1(2,1) + c2(2,1) + c3(2,1);
end
%%
% 计算OIF
OIF = [X_Combination,zeros(size(X_Combination,1),1)];
OIF(:,4) = XCom_StdDev(:,4) ./ XCom_Corr(:,4); % 初步计算出OIF
X_OIF = sortrows(OIF,-4); % 按照OIF值进行降序排序
fprintf('最优组合波段比为 %d/%d、%d/%d、%d/%d,其OIF值为%.3f\n',XRatio_Combination(X_OIF(1,1),:),...
XRatio_Combination(X_OIF(1,2),:),XRatio_Combination(X_OIF(1,3),:),X_OIF(1,4)); % 波段比
% fprintf('最优组合波段为 %d、%d、%d,其OIF值为%.3f\n',X_OIF(1,:)); % 原始波段