GEE主成分分析全解析

0.背景

主成分分析作為數據降維的重要方法，目前中文網站上沒有完整的GEE代碼與教程。而我的畢業論文也使用到了主成分法，因此和它很有感情，就寫下了這篇博客。

1.介紹

主成分分析是將眾多具有相關性的數據指標，重新組合成一組新的指標,新形成的指標互不相關，並且前幾個主成分能代表原始數據的大部分信息。

在GEE中，可能會遇到波段數非常多的情況，這時就可以考慮使用主成分分析法只生成兩、三個主成分，減少后續工作量。

2.代碼思路

3.實操

3.1 數據篩選與預處理

這一步主要是選擇研究區(四川省_資陽市_樂至縣)的哨兵影像，並對篩選的數據按照行政邊界進行鑲嵌與裁剪

//篩選數據
var sentImages = ee.ImageCollection(Sentinel)
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10))
.filterDate("2021-08-01", "2021-08-08")
.filterBounds(LeZhiXian);

//鑲嵌與裁剪數據
var sentmosaic = sentImages.mosaic();
var sentImage = sentmosaic.clip(LeZhiXian);

//加載研究區影像圖層
Map.addLayer(sentImage, trueColor, "樂至縣真彩色");

3.2 數據篩選與預處理

選擇需要進行主成分分析的原始影像波段，並且設置主成分分析影像的分辨率等。在進行主成分分析之前，進行預處理（協方差縮減等）。

//需要進行主成分分析的波段選擇
var bands=["B1","B2","B3","B4","B5","B6","B7","B8","B9","B11","B12"]
sentImage =sentImage.select(bands)

// 輸入到主成分函數的參數設置
var region = LeZhiXian;
var image =  sentImage.select(bands);
var scale = 30;
var bandNames = image.bandNames();

//數據平均
var meanDict = image.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: region,
    scale: scale,
    maxPixels: 1e9
});
var means = ee.Image.constant(meanDict.values(bandNames));
var centered = image.subtract(means);

3.3 主成分分析

這部分的代碼是從Google earth engine官方文檔中copy過來的，針對輸入影像，進行轉為數組、計算正交矩陣、計算主成分載荷等操作。

//主成分分析函數
var 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.sqrt())
      .arrayProject([0]).arrayFlatten([getNewBandNames('sd')]);
      
    //將PC轉換為P波段圖像，通過SD標准化。
    principalComponents=principalComponents
      // 拋出一個不需要的維度，[[]]->[]。
      .arrayProject([0])
      // 使單波段陣列映像成為多波段映像，[]->image。
      .arrayFlatten([getNewBandNames('pc')])
      // 通過SDs使PC正常化。
      .divide(sdImage);
    return principalComponents
  };
//進行主成分分析，獲得分析結果
var pcImage = getPrincipalComponents(centered, scale, region);
// 主要成分可視化
Map.addLayer(pcImage, {bands: ['pc3', 'pc2', 'pc1'], min: -2, max: 2}, 'Sentinel 2 - PCA');

3.4 數據導出

選擇需要導出的主成分波段，這里導出的是前三個波段，因為前三個波段累計貢獻率超過95%，完全夠用了。

//選擇導出的波段
var pcImage_output =pcImage.select(['pc1', 'pc2', 'pc3'])
//導出函數
Export.image.toDrive({
  image: pcImage_output,
  description: 'LeZhiXian_Sentinel_PAC',
  folder:'LeZhiXian',
  scale: 10,
  region:LeZhiXian,
  maxPixels: 1e10
});

點擊Run進行數據下載，一個縣的研究區大概10分鍾就能得到預期影像。

將下載的影像，導入到arcgis或envi中，進行你所需要的分析。
這里稍微提一下，在深度學習中，也可以用主成分分析法處理多波段影像，獲得三個波段，用於訓練與預測。

4.測試鏈接

https://code.earthengine.google.com/8cdedf188b5e1660c65d84c38767b818