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

2022-01-06 21:32:00 阅读：255 来源： 互联网

标签：Combination OIF Image XRatio 波段遥感 matlab size

为减少信息冗余，同时为了对数据进行降维，采用最佳指数因子（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;
    X_StdDev(i) = sqrt(sum(sum(X_Sq),'omitnan')/(size(X,1) * size(X,2) - 1));
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,:));  % 原始波段

标签：Combination,OIF,Image,XRatio,波段,遥感,matlab,size
来源： https://blog.csdn.net/weixin_43637490/article/details/122292097

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9

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