matlab计算遥感影像最“佳”指数因子OIF

为减少信息冗余，同时为了对数据进行降维，采用最佳指数因子（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,:));  % 原始波段