IDL 实现PCA算法

在多元统计分析中，主成分分析（Principal components analysis，PCA）是一种分析、简化数据集的技术。主成分分析经常用于减少数据集的维数，同时保持数据集中的对方差贡献最大的特征。【wiki】

在遥感影像解译与分类中，PCA是经常用到的降维滤噪处理技术。现在实现这个处理流程，便于熟悉和掌握IDL矩阵乘除运算操作。

PRO PCA,DATA,EIGENVALUES = egValues,EIGENVECTORS = egvec,PERCENT = PERCENT,_EXTRA=EXTRA
  GET_SZ,data,ns=ns,nl=nl,nb=nb,type = type
  
  IF ISA(DATA,/NUMBER) AND NB GT 1 AND SIZE(DATA,/N_DIMENSIONS) EQ 3 THEN BEGIN
    DATA = TRANSPOSE(REFORM(DATA,NS*NL,NB))
    corr = correlate(DATA,/covariance)
    egValues = EIGENQL( corr, EIGENVECTORS=egvec,/DOUBLE ,/ABSOLUTE)
    absEgValue = ABS(egValues)
    PERCENT = absEgValue / TOTAL(absEgValue);
    ; EGVEC FORMAT
    ;EGVEC = [  EGVEC1
    ;           EGVEC2
    ;           EGVEC3
    ;           ...
    
    DATA = TEMPORARY(REFORM(TRANSPOSE(egvec) ## TRANSPOSE(DATA),NS,NL,NB))
  ENDIF
END
;---------------------------------------
pro GET_SZ,data,ns=ns,nl=nl,nb=nb,type = type

  sz = size(data)
  
  type = sz[0] ge 1 ? sz[-2] : 0
  ns = sz[0] ge 1 ? sz[1] : 1
  nl = sz[0] ge 2 ? sz[2] : 1
  nb = sz[0] ge 3 ? sz[3] : 1
  
end

 

处理流程：

1. 获得矩阵 行NL、列NS、波段数nb。ns ,number of samples;nl,number of lines;nb,number of bands。三种缩写借鉴自Envi。

2. 矩阵变形，将3维变成2维，nb行，ns*nl列。也就是说将每一个波段的二维图像矩阵压缩成一维数组。因为求相关系数的函数correlate不支持3维矩阵。

3. 对各波段之间的相关系数矩阵求特征值和对应的特征向量。

4. 如果有需要，可以计算各主成分方差的比重。

5. 特征向量左乘原三维矩阵得到主成分结果。转置、二维变三维都是中间过程，函数用法可以查帮助，在此是次要细节略过。

说明:

• 如此，输入的三维矩阵变量data，运算完毕size不变，内容却成了各主成分。这就是IDL procedure的典型用法。若要一次处理多个变量a、b、c、d，只需定义一个pro，将abcd作为参数传入，运算完毕abcd值就都改变了。如果采用函数来return计算结果则需要定义4个function。

• function多用在连续调用，如a=foo1（foo2（foo3（“”）））；这种情形。

• 本计算与Envi菜单的PCA分析结果略有差别。因为相关系数矩阵计算结果不一致，原因暂时我也不清楚，希望大牛可以告知。不过没有太大关系，本文主要是熟悉IDL的编程，PCA的运算过程。只要流程正确、结果有效，在遥感影像解译与分类中发挥作用即可。Envi可能有预处理（如归一化），没必要一定要向它看齐。