当前位置：首页 > news >正文

基于 GEE 的 Sentinel-2 光谱、指数、纹理特征提取与 Sentinel-1 SAR 数据处理

news 2025/10/9 7:58:50

一、代码整体框架与核心功能

二、代码逐段解析

（一）注释与基础配置

（二）Sentinel-2 数据处理模块

1. 数据筛选与预处理

2. 指数计算（植被指数 + 土壤指数）

（1）植被指数（4 类）

（2）土壤指数（6 类）

3. 纹理特征提取（基于近红外 B8 波段）

4. 光谱波段筛选与格式统一

（三）Sentinel-1 数据处理模块

1. 数据筛选与升轨 / 降轨分离

2. 雷达数据合成与融合

（四）数据可视化模块

1. Sentinel-2 相关可视化

2. Sentinel-1 相关可视化

3. 纹理特征可视化

4. 地图中心设置

（五）数据导出模块

（六）信息打印模块

三、运行结果

四、代码关键技术要点与优化建议

（一）关键技术要点

（二）优化建议

五、代码应用场景

若觉得代码对您的研究 / 项目有帮助，欢迎点击打赏支持！需要完整代码的朋友，打赏后可在后台私信（复制文章标题发给我），我会尽快发您完整可运行代码，感谢支持！

一、代码整体框架与核心功能

本代码基于 Google Earth Engine（GEE）平台，实现了 Sentinel-2 光学影像与 Sentinel-1 雷达影像的协同处理，涵盖数据获取、指数计算、纹理特征提取及多源数据融合导出，最终生成包含多维度特征的 10 米分辨率 GeoTIFF 数据，可广泛应用于农业监测、生态评估、土地利用分类等遥感应用场景。

代码遵循 “数据输入→数据处理→特征计算→数据融合→可视化→导出” 的遥感数据处理流程，核心功能分为三大模块，各模块的职责与输出如下表所示：

核心模块	关键任务	输出结果
Sentinel-2 处理	1. 筛选低云量地表反射率数据2. 计算 4 类植被指数、6 类土壤指数3. 提取近红外波段 17 个纹理特征4. 筛选 9 个关键光谱波段	包含 “光谱波段 + 植被指数 + 土壤指数 + 纹理特征” 的光学特征集
Sentinel-1 处理	1. 筛选 IW 模式 SAR 数据2. 分离升轨 / 降轨数据3. 提取 VV/VH 双极化波段4. 多时相中值合成去噪	升轨（ASC_VV/ASC_VH）、降轨（DESC_VV/DESC_VH）雷达特征集
数据融合与导出	1. 合并光学与雷达特征2. 统一数据格式与空间范围3. 可视化关键特征4. 导出 Cloud Optimized GeoTIFF	10 米分辨率、覆盖研究区的多源融合数据集（GeoTIFF）

二、代码逐段解析

（一）注释与基础配置

代码开头的注释块是项目文档的核心，明确了代码的用途、功能边界与数据参数，便于后续维护与复用：

/** 代码名称：GEE获取Sentinel-1/2影像数据及特征计算* 描述：使用Google Earth Engine获取和处理Sentinel-1/2卫星影像数据* 功能：*   1. Sentinel-2数据处理：*      - 获取指定区域的Sentinel-2 SR数据（地表反射率数据）*      - 计算植被指数：NDVI、GNDVI、EVI、SAVI*      - 计算土壤指数：BSI、NDSI、CI、SBI、NDTI、NDRI*      - 计算纹理特征（基于近红外B8波段）：ASM、Contrast等17个特征*   *   2. Sentinel-1数据处理：*      - 获取升轨和降轨SAR数据*      - 提取VV和VH极化波段*      - 生成多时相合成影像* *   3. 数据导出：*      - 导出10m分辨率的多源遥感数据为GeoTIFF格式*      - 包含光学波段、雷达波段、植被指数、土壤指数和纹理特征* * 数据时间：2020-01-01至2020-12-31* 空间分辨率：10米* */// 定义感兴趣区域
var region = table;

关键说明：var region = table 是代码的 “空间入口”，table需提前在 GEE 平台导入（如 Shapefile、KML 等矢量边界），后续所有数据处理均基于该区域裁剪，确保空间一致性。

（二）Sentinel-2 数据处理模块

Sentinel-2 是欧空局的光学卫星，搭载多光谱成像仪，提供 13 个波段，此处重点使用地表反射率数据（SR）（去大气校正后的标准化数据，可直接用于指数计算）。

1. 数据筛选与预处理

// 1. 加载Sentinel-2地表反射率数据集（HARMONIZED版本确保不同传感器数据一致性）
var s2Collection = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED').filterDate('2020-01-01', '2020-12-31') // 时间范围：2020全年.filterBounds(region) // 空间筛选：仅保留覆盖研究区的影像.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 5)) // 质量筛选：云量<5%（减少云干扰）.map(function(image) {return image.clip(region); // 逐影像裁剪：提前裁剪减少后续计算量});// 2. 多时相中值合成（核心去噪步骤）
var s2Median = s2Collection.mean().clip(region);

技术细节：
- mean()（均值合成）：相较于中值（median()），均值更适合连续反射率数据，能平滑季节波动，同时保留地物整体特征（如作物生长季的平均植被覆盖度）；
- COPERNICUS/S2_SR_HARMONIZED：是 Sentinel-2 的 “harmonized” 版本，解决了 Sentinel-2A/B 卫星传感器差异导致的波段偏移问题，确保跨时间序列数据的可比性。

2. 指数计算（植被指数 + 土壤指数）

通过normalizedDifference()（归一化差值）和expression()（自定义公式）计算 10 类指数，覆盖植被生长状态与土壤属性监测需求，各指数的计算逻辑与应用场景如下：

（1）植被指数（4 类）

用于表征植被覆盖度、生物量、光合能力等，核心是利用 “植被对近红外（NIR）强反射、对红光（RED）强吸收” 的光谱特性：

var indices = {// 1. NDVI（归一化植被指数）：最常用植被指数，范围[-1,1]，正值越大植被越茂盛NDVI: s2Median.normalizedDifference(['B8', 'B4']).clip(region), // 2. GNDVI（绿度归一化植被指数）：用绿光（B3）替代红光，对早期植被（如幼苗）更敏感GNDVI: s2Median.normalizedDifference(['B8', 'B3']).clip(region), // 3. EVI（增强植被指数）：引入蓝光（B2）修正大气散射，减少高植被覆盖区的饱和问题（如森林）EVI: s2Median.expression('2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {'NIR': s2Median.select('B8'), // 近红外波段（842nm）'RED': s2Median.select('B4'), // 红光波段（665nm）'BLUE': s2Median.select('B2') // 蓝光波段（490nm）}).clip(region),// 4. SAVI（土壤调整植被指数）：引入土壤调节系数0.5，减少土壤背景（如裸地、稀疏植被区）对指数的干扰SAVI: s2Median.expression('((NIR - RED) * (1 + 0.5))/(NIR + RED + 0.5)', {'NIR': s2Median.select('B8'),'RED': s2Median.select('B4')}).clip(region),// ...后续土壤指数
};

（2）土壤指数（6 类）

利用 “土壤对短波红外（SWIR）强反射、对可见光 / 近红外弱反射” 的特性，表征土壤湿度、有机质含量、盐碱化程度等：

var indices = {// ...前文植被指数// 1. BSI（裸土指数）：正值越大裸土比例越高，用于区分裸地与植被BSI: s2Median.expression('((SWIR1 + RED) - (NIR + BLUE)) / ((SWIR1 + RED) + (NIR + BLUE))', {'SWIR1': s2Median.select('B11'), // 短波红外1（1610nm）'RED': s2Median.select('B4'),'NIR': s2Median.select('B8'),'BLUE': s2Median.select('B2')}).clip(region),// 2. NDSI（归一化差值土壤指数）：反映土壤湿度，值越高土壤越湿润NDSI: s2Median.normalizedDifference(['B11', 'B8']).clip(region), // 3. CI（粘土指数）：SWIR1/SWIR2比值，值越高粘土含量越高（粘土对SWIR2吸收更强）CI: s2Median.expression('(SWIR1 / SWIR2)', {'SWIR1': s2Median.select('B11'),'SWIR2': s2Median.select('B12') // 短波红外2（2190nm）}).clip(region),// 4. SBI（土壤亮度指数）：值越高土壤亮度越高（如盐碱土亮度高于普通土壤）SBI: s2Median.normalizedDifference(['B12', 'B8']).clip(region), // 5. NDTI（归一化差值土壤湿度指数）：专门监测土壤表层湿度，值越高湿度越大NDTI: s2Median.normalizedDifference(['B11', 'B12']).clip(region), // 6. NDRI（归一化差值红边指数）：用红光（B4）和绿光（B3），对土壤有机质敏感NDRI: s2Median.normalizedDifference(['B4', 'B3']).clip(region)
};

共性操作：所有指数均通过clip(region)确保空间范围与研究区一致，避免无效数据干扰。

3. 纹理特征提取（基于近红外 B8 波段）

纹理特征反映影像中像素灰度值的空间分布规律（如均匀性、对比度），可补充光谱特征的不足（如区分 “茂密森林” 与 “稀疏森林”，二者光谱相似但纹理差异大）。代码中基于灰度共生矩阵（GLCM） 提取 17 个纹理特征：

// 1. 预处理：近红外波段（B8）转换为8位无符号整数（GLCM要求输入为整数型）
var nir_int = s2Median.select('B8').multiply(255).toUint8() // 反射率（0-1）×255映射到0-255灰度值.clip(region);// 2. 计算GLCM纹理矩阵
var glcm = nir_int.glcmTexture({size: 3, // 窗口大小：3×3像素（常用窗口，平衡细节与计算效率）kernel: ee.Kernel.square({radius: 1, // 核半径1像素（对应3×3窗口）units: 'pixels',normalize: true // 归一化：消除窗口内像素数量对纹理值的影响})
}).clip(region);// 3. 纹理特征重命名（便于后续识别与使用）
var texture = {'NIR_asm': glcm.select('B8_asm'), // 角二阶矩：反映纹理均匀性，值越高越均匀'NIR_contrast': glcm.select('B8_contrast'), // 对比度：反映像素灰度差异，值越高边界越清晰'NIR_corr': glcm.select('B8_corr'), // 相关性：反映像素间灰度相关性，值越高纹理越有规律'NIR_var': glcm.select('B8_var'), // 方差：反映灰度值离散程度，值越高纹理越粗糙// ...其余13个特征（idm、savg、ent等，均为GLCM标准输出）
};

关键原理：GLCM 通过统计 3×3 窗口内 “相邻像素灰度值对” 的出现频率，计算出 17 个量化指标，每个指标从不同维度描述纹理（如ent（熵）反映纹理复杂度，值越高纹理越复杂）。

4. 光谱波段筛选与格式统一

筛选 Sentinel-2 中10 米 / 20 米分辨率的关键波段（20 米波段会自动重采样至 10 米，与其他数据匹配），并转换为Float32格式（确保数值精度，避免后续计算溢出）：

var spectralBands = s2Median.select(['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'B11']).toFloat().clip(region);
// 波段说明：
// B2(蓝)、B3(绿)、B4(红)、B8(近红外)：10米分辨率，核心光学波段
// B5(红边1)、B6(红边2)、B7(红边3)、B8A(近红外2)：20米分辨率，对植被生长阶段敏感
// B11(短波红外1)：20米分辨率，对土壤湿度敏感

（三）Sentinel-1 数据处理模块

Sentinel-1 是欧空局的合成孔径雷达（SAR）卫星，不受云、雨、光照影响，可全天候监测地物。代码重点处理IW 模式（干涉宽幅模式） 的双极化（VV/VH）数据，用于补充光学数据的不足（如多云地区监测）。

1. 数据筛选与升轨 / 降轨分离

// 1. 加载Sentinel-1地距影像（GRD：经过多视处理的存档数据，可直接使用）
var s1Collection = ee.ImageCollection('COPERNICUS/S1_GRD').filterDate('2020-10-25', '2020-12-31') // 时间范围：仅冬季（可能用于监测冬季作物/土壤冻结）.filterBounds(region).filter(ee.Filter.eq('instrumentMode', 'IW')) // 模式筛选：IW模式（覆盖范围大，分辨率10米）.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV')) // 保留VV极化.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH')); // 保留VH极化// 2. 分离升轨（ASCENDING）与降轨（DESCENDING）数据
var s1Ascending = s1Collection.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'));
var s1Descending = s1Collection.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'));

极化方式说明：
- VV（垂直发射 - 垂直接收）：对地表粗糙度敏感（如裸地 VV 值高于植被）；
- VH（垂直发射 - 水平接收）：对植被结构敏感（如森林 VH 值高于农田）；
- 升轨 / 降轨分离：由于 SAR 的 “侧视成像” 特性，同一地物在升轨 / 降轨影像中的散射特征可能不同，保留两种轨道数据可提高监测准确性。

2. 雷达数据合成与融合

// 1. 多时相中值合成（SAR数据去噪核心：减少雷达斑点噪声）
var s1AscMedian = s1Ascending.mean().clip(region);
var s1DescMedian = s1Descending.mean().clip(region);// 2. 条件融合：仅当升轨/降轨数据存在时，才添加到最终影像（避免空波段报错）
if(s1AscMedian.bandNames().length().getInfo() > 0) {indexImage = indexImage.addBands(s1AscMedian.select('VV').rename('ASC_VV')) // 升轨VV极化.addBands(s1AscMedian.select('VH').rename('ASC_VH')); // 升轨VH极化
}
if(s1DescMedian.bandNames().length().getInfo() > 0) {indexImage = indexImage.addBands(s1DescMedian.select('VV').rename('DESC_VV')) // 降轨VV极化.addBands(s1DescMedian.select('VH').rename('DESC_VH')); // 降轨VH极化
}

技术细节：
- getInfo()：用于获取影像波段数量的 “客户端值”，判断数据是否存在（GEE 中大部分操作是 “服务器端延迟执行”，getInfo()可强制获取即时结果）；
- 雷达数据单位：Sentinel-1 GRD 数据的原始值是 “分贝（dB）”，范围通常为 [-30, 0]，代码中后续可视化设置min: -25, max: 0是基于该范围的合理调整。

（四）数据可视化模块

通过Map.addLayer()将关键数据在 GEE 地图上可视化，便于实时检查数据质量（如是否存在云干扰、裁剪是否正确），各可视化图层的配置逻辑如下：

1. Sentinel-2 相关可视化

// 1. 真彩色合成（RGB）：模拟人眼视觉，用于直观识别地物（如耕地、水体、建筑）
Map.addLayer(s2Median.clip(region), {bands: ['B4', 'B3', 'B2'], // 红（B4）、绿（B3）、蓝（B2）波段组合min: 0, max: 3000, // 反射率范围（Sentinel-2 SR数据的典型有效值）gamma: 1.2 // 伽马校正：增强图像对比度（1.2为常用值，避免过亮或过暗）
}, 'Sentinel-2 RGB');// 2. 植被指数可视化（以NDVI为例）：绿色越深植被越茂盛
Map.addLayer(indices.NDVI.clip(region), {min: -1, max: 1, palette: ['white', 'green']}, 'NDVI');// 3. 土壤指数可视化（以BSI为例）：棕色越深裸土比例越高
Map.addLayer(indices.BSI.clip(region), {min: -1, max: 1, palette: ['white', 'brown']}, 'BSI');

2. Sentinel-1 相关可视化

// 升轨VV极化可视化：值越高（越亮），地表后向散射越强（如建筑、裸地）
Map.addLayer(s1AscMedian, {bands: ['VV'], min: -25, max: 0}, 'S1 Ascending VV');
// 升轨VH极化可视化：植被区VH值通常低于VV值（表现为更暗）
Map.addLayer(s1AscMedian, {bands: ['VH'], min: -25, max: 0}, 'S1 Ascending VH');

3. 纹理特征可视化

// 对比度（NIR_contrast）：红色越深，纹理对比度越高（如道路边缘、建筑边界）
Map.addLayer(texture.NIR_contrast.toFloat().clip(region), {min: 0, max: 1000, palette: ['white', 'red']}, 'NIR Contrast');
// 均匀性（NIR_asm）：蓝色越深，纹理越均匀（如茂密森林、平整农田）
Map.addLayer(texture.NIR_asm.toFloat().clip(region), {min: 0, max: 1, palette: ['white', 'blue']}, 'NIR ASM');

4. 地图中心设置

Map.centerObject(region, 9); // 自动将地图中心定位到研究区，缩放级别9（10米分辨率数据的合理缩放级）

（五）数据导出模块

将融合后的多源数据导出为Cloud Optimized GeoTIFF（COG） 格式（云优化格式，支持云端高效读取，无需下载完整文件），导出参数配置如下：

Export.image.toDrive({image: indexImage, // 待导出的融合影像（包含所有波段/指数/纹理/雷达特征）description: 'Sentinel2_S1_2020_Indices_Texture', // 导出文件名称（清晰标识数据内容）folder: 'GEE_Export', // 导出到Google Drive的文件夹（需提前创建）scale: 10, // 空间分辨率：10米（与Sentinel-2核心波段一致）region: region.geometry(), // 导出范围：研究区矢量边界的几何对象maxPixels: 1e13, // 最大像素数：1e13（支持超大范围导出，避免像素超限报错）fileFormat: 'GeoTIFF', // 文件格式：GeoTIFF（遥感数据标准格式，支持GIS软件读取）formatOptions: {cloudOptimized: true // 启用COG格式：优化云端存储与访问效率}
});

导出注意事项：
- Google Drive 存储空间：需确保 Drive 有足够空间（10 米分辨率、100km² 的区域，导出文件大小约为 100-200MB）；
- 导出时间：数据量越大，导出时间越长（可在 GEE “Tasks” 面板查看进度）。

（六）信息打印模块

通过print()函数在 GEE 控制台输出关键信息，用于数据质量检查与日志记录：

// 1. 打印所有波段名称：确认融合后的数据包含所有预期特征（如光谱、指数、纹理、雷达）
print('所有波段和指数:', indexImage.bandNames());
// 2. 打印Sentinel-2影像数量：判断数据覆盖度（如全年云量<5%的影像数量是否足够）
print('影像数量:', s2Collection.size());
// 3. 打印Sentinel-2影像日期：确认影像时间分布是否均匀（避免集中在某一季节）
print('影像日期:', s2Collection.aggregate_array('system:time_start'));

三、运行结果

S1 Ascending VV：哨兵 - 1 卫星升轨飞行时的 VV 极化雷达影像，和降轨 VV 极化影像配合，能为地物监测提供更多角度的信息。

S1 Ascending VH：哨兵 - 1 卫星升轨飞行时的 VH 极化雷达影像，与降轨 VH 极化影像结合，可更全面地反映地物的散射特征。

S1 Descending VV：哨兵 - 1 卫星降轨飞行时的 VV 极化雷达影像，VV 极化对地表的后向散射特性有独特表现，可用于识别建筑、裸地等。

S1 Descending VH：哨兵 - 1 卫星降轨飞行时的 VH 极化雷达影像，VH 极化对植被等具有特定的散射响应，可用于监测植被结构、地表粗糙度等。

NIR Contrast：近红外波段对比度，属于纹理特征，反映近红外波段影像中像素灰度值的差异程度，值越高说明对比度越强，地物边界等可能更清晰。

NIR ASM：近红外波段角二阶矩，也是纹理特征，体现近红外波段影像纹理的均匀性，值越高说明纹理越均匀。

NIR Correlation：近红外波段相关性，属于纹理特征的一种，反映近红外波段影像中像素间灰度值的相关程度，可辅助分析地物的空间结构等信息。

Sentinel-2 RGB：哨兵 - 2 卫星的真彩色合成影像（红、绿、蓝波段组合），可直观呈现地物的自然色彩，用于识别各类地物（如耕地、森林、水体、建筑物等）。

NDVI：归一化植被指数，是应用最广泛的植被指数之一，用于反映植被的覆盖度、生长状况等，数值越高通常表示植被越茂盛。

NDSI：归一化差值土壤指数，可用于识别土壤、区分土壤与其他地物（如植被、水体等），也能在一定程度上反映土壤的性质。

SBI：土壤亮度指数，能体现土壤的亮度特征，可辅助区分不同类型的土壤或者土壤的不同状态（如盐碱化程度等）。

NDTI：归一化差值土壤湿度指数，主要用于反映土壤的湿度状况，帮助了解土壤水分含量情况。

NDRI：归一化差值红边指数，可用于监测植被生长状态、区分不同植被类型等，对植被的细微变化较为敏感。

四、代码关键技术要点与优化建议

（一）关键技术要点

多源数据协同：通过 “光学（Sentinel-2）+ 雷达（Sentinel-1）” 融合，弥补了光学数据受云干扰、雷达数据缺乏光谱信息的不足，提高了数据适用性；
去噪策略：Sentinel-2 用 “均值合成” 平滑季节波动，Sentinel-1 用 “均值合成” 减少斑点噪声，纹理特征用 “3×3 窗口 GLCM” 平衡细节与效率；
数据一致性：所有数据均通过clip(region)裁剪至研究区，统一为 10 米分辨率、Float32 格式，确保后续分析（如机器学习分类）的兼容性。

（二）优化建议

时间范围细化：Sentinel-1 当前时间范围（2020.10.25-12.31）过短，可调整为与 Sentinel-2 一致的全年，避免季节不匹配；
云量筛选优化：对于多云地区，可放宽云量阈值（如 < 20%），并结合cloudMask函数（基于 Sentinel-2 QA60 波段）进一步去除云像素，而非直接过滤整景影像；
纹理窗口可选：可增加窗口大小参数（如 5×5、7×7），通过print()对比不同窗口的纹理结果，选择最适合研究区地物的窗口；
导出分块：若研究区范围超大（如 > 1000km²），可通过ee.Geometry.split()将区域分块导出，避免单次导出时间过长。