用GEE和Landsat 8数据，5步搞定城市生态健康“体检报告”（附完整代码）

城市生态健康体检实战：用GEE和Landsat 8生成可视化评估报告

城市规划师和环保工作者常常需要快速评估城市生态状况，但传统方法耗时费力。Google Earth Engine（GEE）平台结合Landsat 8数据，为我们提供了一种高效解决方案。本文将手把手教你如何通过5个步骤完成城市生态健康"体检"，并生成专业级的可视化报告。

1. 准备工作与环境搭建

在开始之前，我们需要准备好GEE账号和基础开发环境。GEE是一个强大的地理空间分析平台，集成了PB级卫星影像数据，特别适合处理大范围、长时间的遥感分析任务。

基础环境配置步骤：

1. 访问Google Earth Engine官网注册开发者账号
2. 登录后进入代码编辑器界面
3. 新建一个空白脚本文件

提示：GEE目前对学术和非商业用途免费开放，但商业使用需要申请企业版权限

对于初次接触GEE的用户，建议先熟悉平台的基本操作：

• 左侧面板：脚本管理、资产管理和搜索功能
• 中间区域：代码编辑窗口
• 右侧面板：地图显示和图层控制
• 上方工具栏：运行、保存和分享按钮

// 测试GEE环境是否正常工作 print('Hello, Google Earth Engine!'); Map.setOptions('HYBRID'); // 设置底图为卫星混合地图

2. 定义研究区域与数据准备

确定研究区域是整个分析的第一步。在GEE中，我们可以通过多种方式定义研究区边界，最常用的是绘制多边形或导入已有的Shapefile文件。

研究区定义技巧：

• 城市尺度分析：建议范围在50×50公里以内
• 区域尺度分析：可适当扩大范围，但需考虑计算资源
• 时间范围选择：Landsat 8数据从2013年至今，建议选择完整年度的数据

// 示例：定义青岛市黄岛区为研究区域 var roi = ee.Geometry.Polygon( [[[120.121, 35.975], [120.121, 35.886], [120.257, 35.886], [120.257, 35.975]]], null, false ); Map.centerObject(roi, 12); // 将地图中心定位到研究区，缩放级别12

Landsat 8数据筛选与预处理：

Landsat 8数据在GEE中有多个版本，我们选择地表反射率产品(TOA)以保证数据质量。云覆盖是影响遥感分析的主要因素，需要特别处理。

// 定义去云函数 function removeCloud(image) { var qa = image.select('BQA'); var cloudMask = qa.bitwiseAnd(1 << 4).eq(0); var cloudShadowMask = qa.bitwiseAnd(1 << 8).eq(0); var valid = cloudMask.and(cloudShadowMask); return image.updateMask(valid); } // 加载Landsat 8数据并进行预处理 var L8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA') .filterBounds(roi) .filterDate('2018-01-01', '2019-12-31') .filterMetadata('CLOUD_COVER', 'less_than', 50) .map(function(image) { return image.set('year', ee.Image(image).date().get('year')); }) .map(removeCloud);

3. 生态指标计算与标准化

遥感生态指数(RSEI)由四个核心指标构成，分别反映城市生态的不同方面。我们需要依次计算这些指标并进行标准化处理，确保它们具有可比性。

四大生态指标说明：

// 计算四大生态指标 var L8imgCol = L8imgCol.map(function(img){ // 湿度指标计算 var Wet = img.expression( 'B*(0.1509) + G*(0.1973) + R*(0.3279) + NIR*(0.3406) + SWIR1*(-0.7112) + SWIR2*(-0.4572)',{ 'B': img.select(['B2']), 'G': img.select(['B3']), 'R': img.select(['B4']), 'NIR': img.select(['B5']), 'SWIR1': img.select(['B6']), 'SWIR2': img.select(['B7']) }); img = img.addBands(Wet.rename('WET')); // 绿度指标(NDVI)计算 var ndvi = img.normalizedDifference(['B5', 'B4']); img = img.addBands(ndvi.rename('NDVI')); // 干度指标(NDBSI)计算 var ibi = img.expression( '(2 * SWIR1 / (SWIR1 + NIR) - (NIR / (NIR + RED) + GREEN / (GREEN + SWIR1))) / (2 * SWIR1 / (SWIR1 + NIR) + (NIR / (NIR + RED) + GREEN / (GREEN + SWIR1)))', { 'SWIR1': img.select('B6'), 'NIR': img.select('B5'), 'RED': img.select('B4'), 'GREEN': img.select('B3') }); var si = img.expression( '((SWIR1 + RED) - (NIR + BLUE)) / ((SWIR1 + RED) + (NIR + BLUE))', { 'SWIR1': img.select('B6'), 'NIR': img.select('B5'), 'RED': img.select('B4'), 'BLUE': img.select('B2') }); var ndbsi = (ibi.add(si)).divide(2); return img.addBands(ndbsi.rename('NDBSI')); }); // 数据标准化处理 var img_normalize = function(img){ var minMax = img.reduceRegion({ reducer: ee.Reducer.minMax(), geometry: roi, scale: 1000, maxPixels: 10e13, }); var year = img.get('year'); var normalize = ee.ImageCollection.fromImages( img.bandNames().map(function(name){ name = ee.String(name); var band = img.select(name); return band.unitScale( ee.Number(minMax.get(name.cat('_min'))), ee.Number(minMax.get(name.cat('_max'))) ); }) ).toBands().rename(img.bandNames()).set('year', year); return normalize; }; var imgNorcol = L8imgCol.map(img_normalize);

4. 主成分分析与RSEI计算

主成分分析(PCA)是RSEI计算的核心步骤，它能将多个相关指标转换为少数几个不相关的主成分，保留原始数据的大部分信息。

PCA在生态评价中的优势：

• 消除指标间的冗余信息
• 自动确定各指标权重
• 简化复杂生态系统评估过程
• 结果具有明确的数学意义

// 主成分分析实现 var pca = function(img){ var bandNames = img.bandNames(); var region = roi; var year = img.get('year'); // 数据中心化 var meanDict = img.reduceRegion({ reducer: ee.Reducer.mean(), geometry: region, scale: 1000, maxPixels: 10e13 }); var means = ee.Image.constant(meanDict.values(bandNames)); var centered = img.subtract(means).set('year', year); // 计算协方差矩阵和特征值分解 var arrays = centered.toArray(); var covar = arrays.reduceRegion({ reducer: ee.Reducer.centeredCovariance(), geometry: region, scale: 1000, bestEffort: true, maxPixels: 10e13 }); var covarArray = ee.Array(covar.get('array')); var eigens = covarArray.eigen(); // 计算主成分 var eigenVectors = eigens.slice(1, 1); var arrayImage = arrays.toArray(1); var principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage); return principalComponents .arrayProject([0]) .arrayFlatten([['PC1', 'PC2', 'PC3', 'PC4']]) .set('year', year); }; var PCA_imgcol = imgNorcol.map(pca); Map.addLayer(PCA_imgcol.first().select('PC1'), {min: -2, max: 2}, 'First Principal Component');

RSEI计算与可视化：

第一主成分(PC1)通常包含了原始数据的大部分信息，我们可以用它来计算最终的遥感生态指数。

// RSEI计算与归一化 var RSEI_imgcol = PCA_imgcol.map(function(img){ img = img.addBands(ee.Image(1).rename('constant')); var rsei = img.expression('constant - pc1', { constant: img.select('constant'), pc1: img.select('PC1') }); rsei = img_normalize(rsei); return img.addBands(rsei.rename('RSEI')); }); // 可视化参数设置 var visParam = { min: 0, max: 1, palette: [ 'd7191c', 'fdae61', 'ffffbf', 'a6d96a', '1a9641' ] }; Map.addLayer(RSEI_imgcol.first().select('RSEI'), visParam, 'RSEI');

5. 结果解读与应用实践

获得RSEI结果后，如何正确解读并将其应用于实际工作中是关键。RSEI值范围在0-1之间，值越大表示生态状况越好。

RSEI分级标准参考：

结果应用场景：

• 城市规划：识别生态敏感区，优化城市空间布局
• 环境评估：监测生态变化趋势，评估政策效果
• 生态修复：确定优先治理区域，指导修复工程
• 学术研究：支持城市生态系统的长期监测研究

// 统计不同生态等级的面积比例 var areaStats = function(image) { var area = ee.Image.pixelArea().addBands(image); var stats = area.reduceRegion({ reducer: ee.Reducer.sum().group({ groupField: 1, groupName: 'level', }), geometry: roi, scale: 30, maxPixels: 1e13 }); return stats; }; // 定义分级函数 var classifyRSEI = function(image) { return image.select('RSEI').gt(0.8).multiply(1) .add(image.select('RSEI').gt(0.6).multiply(1)) .add(image.select('RSEI').gt(0.4).multiply(1)) .add(image.select('RSEI').gt(0.2).multiply(1)) .rename('RSEI_class'); }; var classified = RSEI_imgcol.map(classifyRSEI); var stats = classified.map(areaStats); print('Area statistics by RSEI level:', stats);

报告生成技巧：

1. 使用GEE的导出功能将结果保存为GeoTIFF或CSV格式
2. 在QGIS或ArcGIS中进行进一步制图
3. 结合行政区划数据生成分区统计报告
4. 制作时间序列动画展示生态变化过程