GEE遥感利器：主成分分析（PCA）原理与实战，一键实现降维处理！

GEE遥感利器：主成分分析（PCA）原理与实战，一键实现降维处理！

主成分分析（PCA）是遥感影像处理的经典方法，通过降维去除冗余信息，提取关键特征。今天我们将深入探讨PCA的数学原理，并在GEE中实现这一强大工具。

一、主成分分析：从原理到本质

核心目标：将多波段数据转换为少数几个互不相关的综合波段（主成分），最大限度保留原始信息。

核心思想：

1. 数据线性变换，寻找新坐标轴（主成分方向）
2. 第一主成分承载最大方差，后续成分依次递减
3. 各主成分之间互不相关（协方差为0）

二、PCA数学推导：步步深入

步骤0：准备原始数据

var scale = 1000;
var roi = ee.Geometry.Polygon(
        [[[98.76874122458284, 27.476046782884104],
          [98.76874122458284, 26.36424630156598],
          [100.16400489645784, 26.36424630156598],
          [100.16400489645784, 27.476046782884104]]], null, false);
Map.centerObject(roi, 8);
var hAxis = {'title': 'NDVI','viewWindow': {min: -0.8, max: 0.8}};
var vAxis = {'title': 'EVI', 'viewWindow': {min: -0.5, max: 0.5}};
var imageCol =ee.ImageCollection('MODIS/061/MOD13A1')
                       .select(['NDVI', 'EVI'])
                       .filterBounds(roi)
                       .filterDate("2020-01-01","2020-12-31")
                      .map(function (image){
                        return image.addBands(image.divide(10000), null, true);
                      });
var image = imageCol.median();
var NDVI_EVI = image.reduceRegion({
  reducer: ee.Reducer.toList(), 
  geometry: roi, 
  scale: 1000
});
var n = 100;
var x_origin = ee.Array(NDVI_EVI.get('NDVI')).slice(0, 0, n);
var y_origin = ee.Array(NDVI_EVI.get('EVI')).slice(0, 0, n);
var chart = ui.Chart.array.values(y_origin, 0, x_origin)
    .setSeriesNames(['原始数据'])
    .setOptions({
      title: '原始数据',
      hAxis: hAxis,
      vAxis: vAxis,
      pointSize: 5,
});
print(chart);

步骤1：数据标准化

对每个波段 $X_k$ 进行中心化处理：
$$Z_k=X_k-{\mu }_k,{\mu }_k=\frac{1}{n}\sum _{i=1}^n{x}_i$$

var x_mean = x_origin.reduce(ee.Reducer.mean(), [0]);
var y_mean = y_origin.reduce(ee.Reducer.mean(), [0]);
var x_centered = x_origin.subtract(x_mean.repeat(0, n));
var y_centered = y_origin.subtract(y_mean.repeat(0, n));
var chart = ui.Chart.array.values(y_centered, 0, x_centered)
    .setSeriesNames(['中心化数据'])
    .setOptions({
      title: '中心化数据',
      hAxis: hAxis,
      vAxis: vAxis,
      pointSize: 5,
});
print(chart);

步骤2：计算协方差矩阵

设标准化数据矩阵为 $Z$（m个像元 × n个波段），协方差矩阵为：
$$C=\frac{1}{m-1}{Z}^{T}Z$$

var D = ee.Array.cat([x_centered, y_centered], 1).transpose();
var covar = D.matrixMultiply(D.transpose()).divide(n-1);
print('covar',covar);

步骤3：特征值分解

解协方差矩阵的特征方程：
$$C{\mathbf{v}}_i={\lambda }_i{\mathbf{v}}_i$$
其中：

• ${\lambda }_i$：特征值（主成分方差）
• ${\mathbf{v}}_i$：特征向量（主成分方向）

var eigens = covar.eigen();
// 获取特征值
var eigenValues = eigens.slice(1, 0, 1);
print('eigenValues',eigenValues);
// 获取特征向量
var eigenVectors = eigens.slice(1, 1);
print('eigenVectors',eigenVectors);

步骤4：生成主成分

第k主成分计算公式：
$$P{C}_k=Z\cdot v_k$$

// 计算主成分
var principalComponents = eigenVectors.matrixMultiply(D);
print('principalComponents',principalComponents);
var x = principalComponents.slice(0, 0, 1).project([1]);
var y = principalComponents.slice(0, 1, 2).project([1]);
var chart = ui.Chart.array.values(y, 0, x)
    .setSeriesNames(['principalComponents'])
    .setOptions({
      title: 'principalComponents',
     hAxis: hAxis,
      vAxis: vAxis,
      // lineWidth: 1,
      pointSize: 5,
});
print(chart);

关键点：特征值大小反映主成分的信息量占比。

三、为什么必须用协方差矩阵？

PCA的数学本质

PCA的首要目标是找到使标准化数据$Z$投影方差最大化的方向。
假设这个方向为单位向量$w$（$w^{T}w=1$）,那么每个数据在这个方向下的坐标值为：$w^{T}Z$，于是有方差：
$$\begin{array}{c}D(x)=\frac{1}{m-1}{\sum }_{i=1}^m{\left(w^T{x}_i\right)}^2\\ =\frac{1}{m-1}\left(w^{T}Z\right){\left(w^{T}Z\right)}^T\\ =\frac{1}{m-1}{w}^{T}Z{Z}^{T}w\\ =w^T\left(\frac{1}{m-1}Z{Z}^T\right)w\\ =w^{T}Cw\end{array}$$
在约束条件 $w^{T}w=1$ 下求极值，构造拉格朗日函数：
$$\mathcal{L}(w,\lambda )=w^{T}Cw-\lambda (w^{T}w-1)$$
对 $w$ 求导并令导数为0：
$$\frac{\partial \mathcal{L}}{\partial w}=2Cw-2\lambda w=0$$
得到关键方程：
$$\boxed{Cw=\lambda w}$$
这正是协方差矩阵 $C$ 的特征值方程！
代入方差$D(x)$中：
$$D(x)=w^{T}Cw=w^T\lambda w=\lambda$$
因此所求方差$D(x)$最大值就是协方差矩阵$C$的最大特征值$\lambda$，而此时的方向就是最大特征值所对应的特征向量$w$，那么次大方向自然就是第二大特征值对应的特征向量，依此类推，直到我们找出了n个，以这n个特征向量作为新的坐标系，然后将原始数据投影到这个坐标系下即可得到降维后的数据。

三、GEE实战：将PCA应用到影像

完整代码示例

var scale = 1000;
var roi = ee.Geometry.Polygon(
        [[[98.76874122458284, 27.476046782884104],
          [98.76874122458284, 26.36424630156598],
          [100.16400489645784, 26.36424630156598],
          [100.16400489645784, 27.476046782884104]]], null, false);
Map.centerObject(roi, 8);
var visualization = {
  min: 0,
  max: 0.9,
  palette: [
    'ffffff', 'ce7e45', 'df923d', 'f1b555', 'fcd163', '99b718', '74a901',
    '66a000', '529400', '3e8601', '207401', '056201', '004c00', '023b01',
    '012e01', '011d01', '011301'
  ],
};
var imageCol =ee.ImageCollection('MODIS/061/MOD13A1')
                       .select(['NDVI', 'EVI'])
                       .filterBounds(roi)
                       .filterDate("2020-01-01","2020-12-31")
                      .map(function (image){
                        return image.addBands(image.divide(10000), null, true);
                      });
var image = imageCol.median().clip(roi);
var bandNames = image.bandNames();
Map.addLayer(image.select('NDVI'), visualization, 'NDVI');   
Map.addLayer(image.select('EVI'), visualization, 'EVI');   
// 计算均值
var meanDict = image.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: roi,
    scale: scale,
    maxPixels: 1e13
});
// 中心化
var means = ee.Image.constant(meanDict.values(bandNames));
var centered = image.subtract(means);

// 转换为数组
var arrays = centered.toArray();
// 计算协方差矩阵
var covar = arrays.reduceRegion({
    reducer: ee.Reducer.centeredCovariance(),
    geometry: roi,
    scale: scale,
    maxPixels: 1e13
});
// 特征值分解
var covarArray = ee.Array(covar.get('array'));
print('covarArray', covarArray);
var eigens = covarArray.eigen();
print('eigens', eigens);
var eigenValues = eigens.slice(1, 0, 1);
print('eigenValues', eigenValues);
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);   
// 获取特征向量
var eigenVectors = eigens.slice(1, 1);
print('eigenVectors', eigenVectors);
// 计算主成分
var arrayImage = arrays.toArray(1);

var principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage).arrayProject([0]).arrayFlatten([['PC1', 'PC2']]);
Map.addLayer(principalComponents.select('PC1'), {min:-1,max:1}, 'PC1');  
Map.addLayer(principalComponents.select('PC2'), {min:-1,max:1}, 'PC2'); 

关键函数解析

ee.Reducer.centeredCovariance()：

• 计算协方差矩阵

eigen()：

• 执行矩阵的特征分解

输出结果：