【GEE】计算多波段的灰度共生矩阵纹理特征并选择波段

---

前言

---------------------------------------------新手GEE学习记录-----------------------------------------
本文介绍如何在GEE上对影像多波段分别计算灰度共生矩阵glcm及其纹理特征，并在计算后的庞大集合中进行筛选。

一、GLCM计算和纹理特征提取

这部分的重头戏是实现了从一幅影像的多个波段计算出的多种特征中，选择某个特征的所有波段。
举例：如我的影像是10波段，计算得到180个纹理波段，以下代码可以实现任意提取b1-b10的asm特征

//-------------------------------Calculate Texture Features-------------------------- 
//glcmTexture使用前影像必须转为int，为减少精度的丢失，先乘1000再除以1000
//glcmTexture参数含义请自行查询
//l8s2_fusion类型需为Image而非ImageCollection
var glcm_op = l8s2_fusion.unitScale(0,0.20).multiply(1000).toInt().glcmTexture({size: 1,kernel:null}).multiply(0.001);

//输出结果的纹理特征是18种，这里从中选择了7种

//这一步得到的是所有波段计算所得的asm特征波段，结果为List，其元素是想要筛选的波段的bandname，可以直接用在影像的select函数中
var asm_bn = glcm_op.bandNames().filter(ee.Filter.stringEndsWith('item','_asm'));
var contrast_bn = glcm_op.bandNames().filter(ee.Filter.stringEndsWith('item','_contrast'));
var corr_bn = glcm_op.bandNames().filter(ee.Filter.stringEndsWith('item','_corr'));
var varr_bn = glcm_op.bandNames().filter(ee.Filter.stringEndsWith('item','_var'));
var idm_bn = glcm_op.bandNames().filter(ee.Filter.stringEndsWith('item','_idm'));
var sent_bn = glcm_op.bandNames().filter(ee.Filter.stringEndsWith('item','_sent'));
var ent_bn = glcm_op.bandNames().filter(ee.Filter.stringEndsWith('item','_ent'));

//7种特征波段合在一起
var glcm_bn = asm_bn.cat(contrast_bn).cat(corr_bn).cat(varr_bn).cat(idm_bn).cat(sent_bn).cat(ent_bn);
print('glcm_bn',glcm_bn);
//选出这些波段组成新的Image
var glcm_select = glcm_op.select(glcm_bn);
print('glcm_select',glcm_select);

二、主成分分析（PCA）降维

这一步也是迫不得已，若不对纹理特征降维，一个10波段影像会得到180个纹理特征波段，超出GEE计算内存。上一步的选择已经降到70个波段，再做一次主成分分析，只保留贡献率最大的前几个波段即可。

1.PCA函数

参考来自：https://zhuanlan.zhihu.com/p/419712672
但是这位大佬的代码稍有问题，我做了一点修改：

1. 改成脚本函数形式
2. 增加了辅助函数getNewBandNames
3. “var sdImage = ee.Image(eigenValues.abs().sqrt())”这一行加了一个abs()，用来解决“Array: Max (NaN) cannot be less than min (NaN)”这个报错，详见：https://blog.csdn.net/qq_31988139/article/details/126467812

以下两个函数必须放在一个脚本内：

//Used in PCA function
// This helper function returns a list of new band names.
var getNewBandNames = function(prefix, bandNames) {
  var seq = ee.List.sequence(1, bandNames.length());
  return seq.map(function(b) {
    return ee.String(prefix).cat(ee.Number(b).int());
  });
};

//主成分分析函数
exports.getPrincipalComponents = function(centered, scale, region) {
    // 图像转为一维数组
    var arrays = centered.toArray();

    // 计算相关系数矩阵
    var covar = arrays.reduceRegion({
      reducer: ee.Reducer.centeredCovariance(),
      geometry: region,
      scale: scale,
      maxPixels: 1e9
    });

    // 获取“数组”协方差结果并转换为数组。
    // 波段与波段之间的协方差
    var covarArray = ee.Array(covar.get('array'));

    // 执行特征分析，并分割值和向量。
    var eigens = covarArray.eigen();

    // 特征值的P向量长度
    var eigenValues = eigens.slice(1, 0, 1);

    //计算主成分载荷
    var eigenValuesList = eigenValues.toList().flatten();
    var total = eigenValuesList.reduce(ee.Reducer.sum());
    var percentageVariance = eigenValuesList.map(function(item) {
      return (ee.Number(item).divide(total)).multiply(100).format('%.2f');
    });

    print("各个主成分的所占总信息量比例", percentageVariance) ; 

    // PxP矩阵，其特征向量为行。
    var eigenVectors = eigens.slice(1, 1);

    // 将图像转换为二维阵列
    var arrayImage = arrays.toArray(1);

    //使用特征向量矩阵左乘图像阵列
    var principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage);

    // 将特征值的平方根转换为P波段图像。
    var sdImage = ee.Image(eigenValues.abs().sqrt())
      .arrayProject([0]).arrayFlatten([getNewBandNames('sd', centered.bandNames())]);

    //将PC转换为P波段图像，通过SD标准化。
    principalComponents=principalComponents
      // 抛出一个不需要的维度，[[]]->[]。
      .arrayProject([0])
      // 使单波段阵列映像成为多波段映像，[]->image。
      .arrayFlatten([getNewBandNames('pc', centered.bandNames())])
      // 通过SDs使PC正常化。
      .divide(sdImage);
    return principalComponents;
  };

2.应用PCA降维

使用该函数会打印各个主成分的所占总信息量比例，我的数据运行结果显示前3个波段累计贡献率达到99%以上，所以只保留前三个波段：

各个主成分的所占总信息量比例：
JSON
List (70 elements)
0:
72.93
1:
22.13
2:
4.40
3:
0.37
4:
0.14

//Too many texture bands, we have to do PCA and just retain main bands
var pcImage = fp.getPrincipalComponents(glcm_select, 20, roi);
//选择保留的波段
var pcImage_output = pcImage.select(['pc1', 'pc2', 'pc3']).rename(['op_PC1','op_PC2','op_PC3']);

---