用Google Earth Engine玩转30年全球夜光数据从DMSP到VIIRS的完整分析指南夜光遥感数据如同一面镜子映照着人类文明的足迹。从1992年至今DMSP和VIIRS卫星持续记录着地球夜晚的灯光变化这些数据已成为研究城市化、能源消费和经济发展的重要指标。但对于许多GIS从业者和数据科学家来说如何高效处理这份跨越30年的数据集仍是一个挑战。本文将带你从零开始在Google Earth Engine平台上完成夜光数据的加载、分析和可视化全流程。1. 夜光数据基础与GEE环境准备夜光遥感数据主要来自两个卫星系统DMSP-OLS1992-2013年和VIIRS2012年至今。由于传感器差异原始数据存在明显的不连续性。而经过协调处理的Harmonized_NTL数据集解决了这一问题实现了1992-2021年间1km分辨率的数据一致性。在开始分析前你需要注册Google Earth Engine账号免费学术用途熟悉基本JavaScript语法了解GEE的基本操作界面提示GEE的代码编辑器可通过浏览器直接访问所有计算在云端完成无需本地高性能设备。2. 数据加载与初步探索在GEE中加载协调后的夜光数据集非常简单// 加载DMSP和VIIRS协调数据集 var dmsp ee.ImageCollection(projects/sat-io/open-datasets/Harmonized_NTL/dmsp); var viirs ee.ImageCollection(projects/sat-io/open-datasets/Harmonized_NTL/viirs); // 合并两个数据集 var ntlCollection dmsp.merge(viirs);数据集关键属性对比如下属性DMSP (1992-2013)VIIRS (2012-2021)协调后数据原始分辨率2.7km750m1kmDN值范围0-630-∞标准化处理过饱和问题严重较轻已校正初步可视化代码示例// 选择2010年DMSP和2020年VIIRS数据 var image2010 dmsp.filterDate(2010-01-01, 2010-12-31).first(); var image2020 viirs.filterDate(2020-01-01, 2020-12-31).first(); // 可视化参数 var visParams { min: 0, max: 30, palette: [000000, 0d0887, 7e03a8, cb4679, f89441, f0f921] }; // 添加图层到地图 Map.addLayer(image2010, visParams, 2010 DMSP); Map.addLayer(image2020, visParams, 2020 VIIRS);3. 时间序列分析与变化检测夜光数据最强大的应用之一是分析长时间序列变化。以下是计算区域夜光总量变化的完整代码// 定义分析区域以长三角为例 var roi ee.Geometry.Rectangle([118, 29, 123, 33]); // 计算年度合成 function createAnnualComposite(year) { var start ee.Date.fromYMD(year, 1, 1); var end start.advance(1, year); return ntlCollection.filterDate(start, end) .mean() .set(year, year); } // 生成1992-2021年时间序列 var years ee.List.sequence(1992, 2021); var timeSeries ee.ImageCollection.fromImages( years.map(function(year) { return createAnnualComposite(year); }) ); // 计算区域统计 var chart ui.Chart.image.series({ imageCollection: timeSeries.select(NTL), region: roi, reducer: ee.Reducer.sum(), scale: 1000, xProperty: year }).setOptions({ title: 长三角地区夜光总量变化 (1992-2021), hAxis: {title: 年份}, vAxis: {title: DN值总和}, lineWidth: 3, pointSize: 5 }); print(chart);常见分析指标包括夜光总量Total DN亮区面积Area with DNthreshold平均亮度Mean DN亮度重心Light centroid4. 城市扩张与经济发展分析夜光数据与城市化进程密切相关。以下代码演示如何提取城市建成区// 定义城市提取阈值DN7 function extractUrbanArea(image) { var urban image.gt(7).rename(urban); return urban.set(year, image.get(year)); } // 应用函数到时间序列 var urbanCollection timeSeries.map(extractUrbanArea); // 计算城市面积变化 var areaChart ui.Chart.image.series({ imageCollection: urbanCollection, region: roi, reducer: ee.Reducer.sum().multiply(1e-6), // 转换为km² scale: 1000, xProperty: year }).setOptions({ title: 长三角城市建成区面积变化, hAxis: {title: 年份}, vAxis: {title: 面积 (km²)} }); print(areaChart);实际应用中可以结合其他数据验证分析结果与官方统计年鉴中的建成区数据对比与Landsat分类结果交叉验证结合GDP、人口等社会经济数据建立回归模型5. 高级应用与技巧5.1 数据融合与增强为提高分析精度可将夜光数据与其他数据集融合// 加载人口数据集 var popData ee.ImageCollection(WorldPop/GP/100m/pop); // 计算夜光与人口的相关性 var combined timeSeries.map(function(image) { var year image.get(year); var pop popData.filter(ee.Filter.calendarRange(year, year, year)) .first() .resample(bilinear) .reproject(image.projection()); return image.addBands(pop); }); // 计算相关系数 var correlation combined.reduce(ee.Reducer.pearsonsCorrelation());5.2 异常检测与事件分析夜光数据可用于检测特殊事件如疫情封锁// 分析2020年武汉夜光变化 var wuhan ee.Geometry.Point([114.3055, 30.5928]).buffer(50000); // 计算月度变化 var monthly ntlCollection.filterDate(2019-01-01, 2021-12-31) .map(function(image) { return image.reduceRegion({ reducer: ee.Reducer.mean(), geometry: wuhan, scale: 1000 }).set(date, image.date()); }); var monthlyChart ui.Chart.feature.byFeature(monthly, date, [NTL]) .setOptions({ title: 武汉地区月度夜光变化, hAxis: {title: 日期}, vAxis: {title: 平均DN值}, lineWidth: 2 }); print(monthlyChart);5.3 数据导出与本地分析对于需要本地处理的情况可以导出数据// 导出2020年全球夜光数据 Export.image.toDrive({ image: viirs.filterDate(2020-01-01, 2020-12-31).mean(), description: Global_NTL_2020, scale: 1000, region: roi, maxPixels: 1e13 });6. 常见问题与解决方案在实际使用中可能会遇到以下典型问题数据缺失问题某些年份DMSP数据缺失解决方案使用邻近年份插值DN值解释问题不同传感器DN值含义不同解决方案参考论文中的标准化方法云平台限制GEE有计算配额限制解决方案优化代码使用clip()限制分析区域跨传感器一致性DMSP与VIIRS交接年份(2012-2014)数据差异解决方案使用协调后的数据集// 优化计算效率的示例代码 function efficientCalculation(image) { return image.clip(roi) // 先裁剪区域 .reduceRegion({ reducer: ee.Reducer.mean(), geometry: roi, scale: 1000, bestEffort: true }); }7. 案例研究粤港澳大湾区夜光演变让我们通过一个完整案例展示分析流程// 定义大湾区范围 var gba ee.Geometry.Polygon( [[[111.5, 21.5], [111.5, 24.5], [115.5, 24.5], [115.5, 21.5]]]); // 计算夜光重心 function calculateCentroid(image) { var year image.get(year); var centroid image.reduceRegion({ reducer: ee.Reducer.centroid(), geometry: gba, scale: 1000, maxPixels: 1e9 }); return ee.Feature(null, { year: year, lon: centroid.get(NTL_longitude), lat: centroid.get(NTL_latitude) }); } var centroids timeSeries.map(calculateCentroid); // 可视化重心移动轨迹 Map.addLayer(gba, {color: red}, GBA Region); Map.addLayer(ee.FeatureCollection(centroids), {color: blue}, Light Centroids); // 生成移动动画 var rgbCollection timeSeries.map(function(image) { return image.visualize(visParams).set(year, image.get(year)); }); var videoArgs { dimensions: 800, region: gba, framesPerSecond: 5, bands: [vis-red, vis-green, vis-blue], min: 0, max: 255 }; print(ui.Thumbnail(rgbCollection, videoArgs));通过这个案例我们可以观察到1992-2000年夜光重心向深圳方向移动2000-2010年多中心发展格局形成2010-2021年区域一体化趋势明显8. 扩展应用与前沿方向夜光数据的创新应用不断涌现社会经济指标估算GDP预测贫困程度评估能源消费估算环境与灾害监测战乱地区冲突评估自然灾害影响范围电力中断恢复监测公共卫生应用疾病传播与夜间活动关系疫情防控措施效果评估// 疫情封锁影响评估示例 function calculateLockdownImpact(city, radiusKm) { var roi city.buffer(radiusKm * 1000); var preCovid viirs.filterDate(2019-01-01, 2019-12-31).mean(); var duringCovid viirs.filterDate(2020-01-01, 2020-12-31).mean(); var diff duringCovid.subtract(preCovid) .clip(roi); var loss diff.lt(0).multiply(diff.abs()); var gain diff.gt(0); return { lightLoss: loss.reduceRegion(ee.Reducer.sum(), roi, 1000).get(NTL), lightGain: gain.reduceRegion(ee.Reducer.sum(), roi, 1000).get(NTL) }; }在实际项目中我发现处理长时间序列夜光数据时特别需要注意2013年前后的传感器过渡期。一个实用的技巧是对于跨传感器比较的研究最好选择1995-2010纯DMSP和2015-2021纯VIIRS两个时段分别分析避免过渡期的不确定性。
