用300万张卫星图解码城市生长：GEE实战WHU不透水面数据集全流程指南
站在城市的天台上俯瞰，那些由混凝土、沥青和玻璃构成的"皮肤"正以惊人速度蔓延。这种被称为不透水面的城市表皮，不仅改变了地表水文循环，更重塑着人类与自然的互动方式。武汉大学团队发布的30米分辨率全球不透水面数据集（WHU/GISA），就像一台时间显微镜，让我们能精确观测近半个世纪来城市扩张的每一个细节。本文将带你用Google Earth Engine这把数字手术刀，亲手解剖城市生长的DNA序列。
1. 理解不透水面数据的科学价值
不透水面就像城市的第二层皮肤，其扩张程度直接反映人类活动的强度。传统测绘手段需要耗费数月才能完成的城市扩张分析，现在通过WHU/GISA数据集只需几行代码就能实现时空维度的全面观测。这套数据产品的独特之处在于：
- 多时序光谱特征：融合Landsat、MODIS和ESA_CCI数据源，构建连续42年（1978-2019）的时间序列
- 机器学习验证：基于12万个验证点的评估显示，错检率仅0.82%，漏检率控制在5.16%
- 地理网格化建模：全球按2°格网划分，局部建立独立分类模型，避免"一刀切"带来的精度损失
// 数据集元数据结构示例 var metadata = { id: 'WHU/GISA', resolution: '30m', time_span: '1978-2019', bands: ['B1'], // 单波段存储不透水面概率(0-100) valid_range: [0, 100], no_data_value: 255 };城市规划者常用该数据监测"城市蔓延指数"，环境科学家则用它计算"地表径流系数"。在深圳前海新区的规划中，研究团队发现2015-2019年间不透水面增长导致地表径流增加了37%，这一发现直接影响了该区域排水系统的设计标准。
2. GEE平台数据获取实战
2.1 初始化数据环境
在GEE代码编辑器中，首先需要建立时空过滤条件。不同于普通影像数据集，WHU/GISA采用年度合成方式存储，时间筛选需要特别注意跨年问题：
// 正确的时间筛选方式（整年数据） var collection = ee.ImageCollection('WHU/GISA') .filter(ee.Filter.calendarRange(2015, 2019, 'year')) .select('B1'); // B1波段存储不透水面数据 // 常见错误示例（会导致数据缺失） var wrongFilter = collection.filterDate('2015-06-01', '2019-06-01');提示：该数据集每年生成一张全球合成图，建议使用calendarRange而非filterDate进行年度筛选
2.2 空间范围精准裁剪
当研究特定城市区域时，需要结合行政区划矢量数据进行裁剪。以下示例展示如何获取北京市五环内不透水面变化：
// 加载北京市边界 var beijing = ee.FeatureCollection('users/shapefiles/Beijing'); var fifthRing = beijing.filter(ee.Filter.eq('name', '5th_ring')); // 裁剪并计算年度均值 var clipped = collection.map(function(image) { return image.clip(fifthRing); }); // 时序统计计算 var stats = clipped.reduceColumns({ reducer: ee.Reducer.mean(), selectors: ['B1'] });3. 多维数据分析技巧
3.1 变化检测可视化
通过差值计算可直观呈现城市扩张热点区域。以下代码生成上海2010-2020年不透水面变化热力图：
var shanghai = ee.FeatureCollection('users/shapefiles/Shanghai'); var baseline = ee.ImageCollection('WHU/GISA') .filter(ee.Filter.eq('year', 2010)) .first(); var current = ee.ImageCollection('WHU/GISA') .filter(ee.Filter.eq('year', 2020)) .first(); var change = current.subtract(baseline).clip(shanghai); // 可视化参数 var visParams = { min: -30, // 减少区域 max: 30, // 增加区域 palette: ['blue', 'white', 'red'] }; Map.addLayer(change, visParams, 'Shanghai Change 2010-2020');3.2 统计指标自动化计算
城市扩张分析常需要量化指标，下表列出常用统计指标及其GEE实现方法：
| 指标名称 | 计算公式 | GEE实现方法 | 应用场景 |
|---|---|---|---|
| 不透水面占比 | 不透水像元数/总像元数×100% | reduceRegion+比例计算 | 城市密度评估 |
| 年增长率 | (本年值-上年值)/上年值×100% | 时序差分+百分比计算 | 扩张速度监测 |
| 空间自相关指数 | Moran's I | 邻域分析+空间统计 | 扩张模式识别 |
| 重心迁移距离 | 年度重心坐标差 | 质心计算+距离公式 | 扩张方向分析 |
// 计算杭州市不透水面占比示例 var hangzhou = ee.FeatureCollection('users/shapefiles/Hangzhou'); var stats = image.reduceRegion({ reducer: ee.Reducer.mean(), geometry: hangzhou, scale: 30, maxPixels: 1e9 }); // 转换为百分比 var percentage = ee.Number(stats.get('B1')).multiply(100); print('Impervious surface percentage:', percentage);4. 典型应用场景解析
4.1 城市边界动态监测
结合夜间灯光数据(NPP/VIIRS)和WHU/GISA数据集，可以构建城市边界自动识别算法。成都市的案例显示，该方法识别的城市边界与政府公布的规划边界吻合度达到92%：
- 数据融合：将不透水面概率>50%的区域与夜间灯光值>10nW/cm²/sr的区域取并集
- 形态学处理：使用closing操作填充细小孔洞，opening操作消除孤立噪点
- 边界提取：采用Canny边缘检测算法提取连续边界线
// 城市边界提取代码片段 var urban = image.gt(50).or(nightlight.gt(10)); urban = urban.focal_max(3).focal_min(3); // 形态学处理 var edges = urban.convolve(ee.Kernel.canny(100, 0.5)); Map.addLayer(edges.selfMask(), {palette: 'FF0000'}, 'Urban Edge');4.2 热岛效应关联分析
不透水面比例与地表温度(LST)存在显著相关性。在广州中心城区的分析中，发现不透水面每增加10%，夏季日间地表温度平均上升1.2℃：
// 温度-不透水面回归分析 var samplePoints = urbanArea.randomPoints(1000); var pairedData = lst.addBands(impervious).sampleRegions({ collection: samplePoints, scale: 30 }); // 计算相关系数 var correlation = pairedData.reduceColumns({ reducers: ee.Reducer.pearsonsCorrelation(), selectors: ['LST', 'B1'] }); print('Temperature-Impervious Correlation:', correlation);5. 疑难问题解决方案库
5.1 数据缺失处理
当遇到特定年份数据缺失时，可采用时空插值方法补全：
时间插值：使用前后年份数据的线性组合
var filled = collection.map(function(img) { var year = ee.Date(img.get('system:time_start')).get('year'); return img.set('year', year); }).interpolate({ propertyName: 'B1', propertyTarget: 'year', inputSeries: [2010, 2015, 2020] });空间插值：对云覆盖区域使用邻近像元填充
var filled = image.focal_mean(3, 'circle', 'meters');
5.2 大规模计算优化
处理省级以上区域时，需采用分块计算策略：
// 分块计算示例 var grid = ee.FeatureCollection('users/grids/China_Grid'); var results = grid.map(function(feature) { return image.reduceRegion({ reducer: ee.Reducer.mean(), geometry: feature.geometry(), scale: 30, maxPixels: 1e8 }).set('grid_id', feature.id()); }); // 合并结果 var finalResult = ee.FeatureCollection(results).aggregate_array('mean');注意：当处理超大城市群时，建议使用Export将中间结果导出至Google Drive，避免浏览器内存溢出
6. 进阶分析工作流
6.1 机器学习辅助分类
结合Sentinel-2数据，可以构建更高精度的城市用地分类模型：
特征工程：构建NDVI、NDBI等光谱指数
var addIndices = function(image) { var ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI'); var ndbi = image.normalizedDifference(['B11', 'B8']).rename('NDBI'); return image.addBands([ndvi, ndbi]); };样本生成：基于WHU/GISA生成训练样本
var samples = impervious.gt(70).stratifiedSample({ numPoints: 1000, classBand: 'B1', region: roi });模型训练：使用随机森林算法
var classifier = ee.Classifier.smileRandomForest(50) .train({ features: samples, classProperty: 'B1', inputProperties: ['B2', 'B3', 'B4', 'NDVI', 'NDBI'] });
6.2 三维城市扩张模拟
将时序数据导入Cesium等平台可实现城市生长动画：
// 生成年度序列 var timeline = ee.List.sequence(1978, 2019).map(function(year) { return ee.ImageCollection('WHU/GISA') .filter(ee.Filter.eq('year', year)) .first() .visualize({min: 0, max: 100, palette: ['white', 'black']}) .set('year', year); }); // 导出为GIF Export.video.toDrive({ collection: ee.ImageCollection(timeline), description: 'UrbanGrowth', dimensions: 800, framesPerSecond: 5, region: roi });在重庆市的案例中，这种三维模拟清晰展现了"两江四岸"城市发展轴的演变过程，帮助规划者理解山地城市扩张的特殊模式。
