基于GEE的NDVI与NDBI长时序城市扩张监测实战复盘

去年我在做一个多城市对比的遥感监测项目时最头疼的不是算法本身而是如何在统一坐标系下快速处理八年、七个城市的影像时序数据。当时需要用到NDVI和NDBI来识别城市扩张范围本地下载Landsat影像再逐步处理的话光是云掩膜和大气校正就能耗掉大半个月。后来我把整套流程迁移到GEE上用NDVI和NDBI两套指数在几个小时内完成了土耳其七个主要城市、2018到2025年的年度合成与建成区提取。这篇文章就是那次实战的完整复盘内容包括指数原理、实验设计、GEE代码实现、结果解读方法以及我在跑数据时踩过的几个关键坑。如果你正准备做类似的长时序城市扩张分析或者刚开始接触GEE上的NDVI、NDBI计算这篇文章可以直接作为参考框架来用。1. 为什么城市扩张监测偏偏要看NDVI和NDBI1.1 两个指数的光谱逻辑先讲清楚这两个指数各自在做什么。NDVI归一化差值植被指数的核心是植被在红光波段强烈吸收、在近红外波段强烈反射这个特性。公式是NDVI (NIR - Red) / (NIR Red)叶绿素含量越高红光吸收越强近红外反射越强NDVI值越大。所以NDVI高值区域基本对应茂密植被低值区域对应水体、裸土或硬化地表。NDBI归一化差值建筑指数的思路刚好反着来。城市里的屋顶、柏油路、水泥地这类不透水面在短波红外SWIR1波段的反射率通常高于近红外波段而植被正好相反。公式是NDBI (SWIR1 - NIR) / (SWIR1 NIR)于是NDBI高值区域更可能对应建筑密集区。Zha等人在2003年提出这个指数时主要就是利用它对不透水面的敏感性。1.2 为什么单独用一个不够要组合起来城市扩张监测本质上回答一个非常朴素的遥感问题哪个地方从非城市状态变成了城市状态。如果只用NDVI你能看到的只是植被覆盖减少。但植被减少的原因很多可能是裸地、采石场、翻耕后的农田不一定变成城市。如果只用NDBI麻烦更大因为裸土在SWIR1波段的反射率也很高NDBI很容易把大片裸地判成建成区尤其像安纳托利亚高原这种半干旱区这个问题会直接导致面积统计虚高。所以我在项目中用的核心判据是一个很经典的经验规则当NDBI大于NDVI时该像元被认为是城市像元。逻辑也不复杂——建筑地表的SWIR1反射特征压过了植被特征而裸土的NIR反射通常会比SWIR1更平衡一些不容易同时在NDVI上低到和城市一样。当然这个规则并不是万能的后面会专门讲它在哪些场景下失效。1.3 用GEE做这个项目的三个直接好处以前做这种跨年、跨城市的对比常规流程是去USGS下载影像进ENVI做大气校正算指数再拿进ArcGIS逐市统计最后导出Excel做图表。八个年份七个城市光数据量就够喝一壶。GEE解决三个痛点第一数据就绪。Landsat Collection 2 Level 2已经是地表反射率产品不需要自己跑大气校正。第二弹性计算。ImageCollection的map操作会在服务器端并行处理七城市八年的年度合成和指数计算跑完基本是几十秒到几分钟的量级。第三统计方便。reduceRegion可以直接按行政边界统计像元面积不用导出栅格再在桌面GIS里裁剪求和。换句话说GEE把“数据准备算法计算区域统计”压缩到了同一个工作流里这才是它能做这类长时序多城市分析的根本原因。2. 实验设计土耳其七大城市怎么选、看什么窗口2.1 七个城市的选取逻辑土耳其的城市扩张研究选哪几个城市其实有点讲究。如果只看伊斯坦布尔虽然最能代表超大城市的问题但分析结论没法推广到其他发展阶段的区域。我在项目里选了七个分布在不同地理区、人口规模都是百万级的城市。城市所在区域地形与气候特点城市扩张特征伊斯坦布尔马尔马拉区跨博斯普鲁斯海峡丘陵、地中海气候超大城市填充式扩张与边缘蔓延并存安卡拉中安纳托利亚高原内陆半干旱气候首都沿主要交通走廊向北部、西北扩展伊兹密尔爱琴海区沿海平原地中海气候港口城市沿爱琴海岸线延伸布尔萨马尔马拉区山前平原气候温和工业城市工业园区带动边缘扩张安塔利亚地中海区滨海平原典型地中海气候旅游城市沿海岸线东西两侧蔓延科尼亚中安纳托利亚内陆平原干燥少雨农业区与城市边缘转换裸土干扰明显阿达纳地中海东部冲积平原气候炎热农业平原上的城市农田与建设用地交织这个组合的好处是覆盖了沿海和内陆、工业驱动和旅游驱动、超大城市和中等城市可以在后面的结果对比里看出不同扩张模式在NDVI、NDBI变化上的差异。2.2 时间窗口为什么定在6到9月我建的是2018到2025年的年度序列每年取6月到9月的数据做合成。这个时间窗口不是随手定的有两点考虑。第一土耳其大部分地区属于地中海气候夏季炎热干燥云量少有效观测多。第二夏季是植被生长季NDVI的对比度最强城市和植被在光谱上拉得开阈值判断的可靠性更高。合成方式我选的是中值合成median composite不是均值合成。中值对云的残留、传感器异常值更稳健。均值容易被个别高值或低值拉偏可能在NDBI上产生虚假的高值像元。2.3 数据源选择Landsat 8为主Sentinel-2为辅当时也考虑过直接用Sentinel-2。10米分辨率确实比Landsat的30米精细但有两个问题让我放弃了它作为主数据源。一个是2018年Sentinel-2可用的高质量影像数量不稳定不同城市之间时相覆盖不均匀每年合成结果的可比性会打折扣。另一个是Landsat 8在2013年就开始运行2018到2025年整个时间范围内波段设置完全一致NDVI和NDBI的计算不需要处理传感器差异。所以主流程用Landsat 8 OLI地表反射率产品30米分辨率对城市尺度研究完全够用。Sentinel-2后来只用来做交叉验证比如对某一年的提取结果用10米影像目视抽检看30米的结论靠不靠谱。3. GEE实现全过程从年度合成到面积统计3.1 最容易忽略的一步Landsat C2 L2的缩放系数这一步如果漏了后面所有NDVI和NDBI都是错的但我看到很多GEE教程里根本没提。Landsat Collection 2 Level 2产品里的SR波段存的是整数型数值需要乘以0.0000275再减去0.2才能换算成实际的地表反射率。如果不做这个缩放直接拿来算NDVI虽然分子分母同时缩放会抵消一部分影响但在浮点运算和后续阈值判断上仍然会导致数值偏移尤其在做多时相合成时误差会被中值计算进一步放大。我封装了一个缩放函数放在处理流程最前面function applyScaleFactors(image) { var optical image.select(SR_B.).multiply(0.0000275).add(-0.2); return image.addBands(optical, null, true); }这里addBands(optical, null, true)的true表示用缩放后的波段覆盖原始波段避免影像上同时存在两套SR波段导致后面波段选择混乱。3.2 云掩膜与年度中值合成Landsat C2 L2产品自带一个像素质量波段QA_PIXEL云和云阴影对应固定的bit位。我用的掩膜逻辑是提取bit3云和bit4云阴影把这两个标志位出现的位置遮掉。function maskCloudAndShadow(image) { var qa image.select(QA_PIXEL); var cloudBitMask 1 3; var shadowBitMask 1 4; var mask qa.bitwiseAnd(cloudBitMask).eq(0) .and(qa.bitwiseAnd(shadowBitMask).eq(0)); return image.updateMask(mask); }然后按年份过滤、掩膜、缩放到逐年中值合成function annualMedian(year, region) { var start ee.Date.fromYMD(ee.Number(year), 6, 1); var end ee.Date.fromYMD(ee.Number(year), 9, 1); var imgs ee.ImageCollection(LANDSAT/LC08/C02/T1_L2) .filterDate(start, end) .filterBounds(region) .map(maskCloudAndShadow) .map(applyScaleFactors); var median imgs.median(); return median.set(year, year); }关于窗口边界filterDate是左闭右开区间即6月1日到9月1日之间不含9月1日当天的影像。如果你想把9月初的影像也拿进来可以直接把end改成从YMD(2025, 9, 30)。我建议窗口至少给到三个半月以上因为有些年份局部云量偏高可用像元不足后面面积统计会明显偏小。3.3 指数计算与城市像元判定年度合成影像出来以后计算NDVI、NDBI和MNDWI。MNDWIModified NDWI用的是绿波段和短波红外公式如下MNDWI (Green - SWIR1) / (Green SWIR1)它用来掩膜水体因为土耳其几个沿海城市里伊斯坦布尔的博斯普鲁斯海峡、伊兹密尔海湾这类水体如果不管很容易在NIR和SWIR1的比值关系上产生假信号。用MNDWI0识别水体像元并掩膜掉是后面减少误判的关键一步。function addIndices(image) { var ndvi image.normalizedDifference([SR_B5, SR_B4]).rename(NDVI); var ndbi image.normalizedDifference([SR_B6, SR_B5]).rename(NDBI); var mndwi image.normalizedDifference([SR_B3, SR_B6]).rename(MNDWI); return image.addBands([ndvi, ndbi, mndwi]); }城市像元判定我用的是NDBI大于NDVI这条规则function builtUpMask(image) { var idx addIndices(image); var waterMask idx.select(MNDWI).lt(0); var built idx.updateMask(waterMask) .expression((b(NDBI) b(NDVI)) ? 1 : 0); return built; }updateMask(waterMask)会把MNDWI大于0的水体像元全部遮掉再去做NDBI NDVI的判断这样海峡、海湾和水库不会进来凑数。3.4 逐年逐市的面积统计与导出面积统计的思路是把城市像元掩膜乘以ee.Image.pixelArea()得到每个像元的实际地表面积然后按行政边界求和得到平方米再除以1e6转成平方千米。function builtUpAreaKm2(year, region) { var img annualMedian(year, region); var built builtUpMask(img); var areaM2 built .multiply(ee.Image.pixelArea()) .reduceRegion({ reducer: ee.Reducer.sum(), geometry: region.geometry(), scale: 30, maxPixels: 1e10 }); return ee.Number(areaM2.get(constant)).divide(1e6); }导出环节我用一个嵌套map生成FeatureCollection再用Export.table.toDrive输出CSV。这种方式比在Console里一个一个看结果要靠谱得多因为八年的数据用眼睛看根本看不过来。var provinces ee.FeatureCollection(users/your_world/turkey_provinces); var cities ee.List([Istanbul, Ankara, Izmir, Bursa, Antalya, Konya, Adana]); var years ee.List.sequence(2018, 2025); var expansionFC ee.FeatureCollection( cities.map(function(cityName) { var cityGeom provinces.filter(ee.Filter.eq(province, cityName)); return years.map(function(y) { var area builtUpAreaKm2(y, cityGeom); return ee.Feature(cityGeom.geometry(), { city: cityName, year: y, area_km2: area }); }); }).flatten() ); Export.table.toDrive({ collection: expansionFC, description: turkey_urban_expansion_2018_2025, fileFormat: CSV });需要注意provinces这份数据需要你提前在Assets里上传自己整理的城市边界。没有合适的行政边界的话用GEE公开数据里的任意城市边界也行但一定要保证年份跨度内边界没有明显变化。4. 结果整理与解读面积变化和扩张方向怎么看4.1 年度面积表怎么组织跑完导出后你得到的是一张多行CSV每行对应一个城市某年的建成区面积。下面是我为了展示表格结构放的一组演示数据不是真实计算结果你跑完自己的代码后要把真实数值填进来城市2018面积(km²)2025面积(km²)增幅伊斯坦布尔1520170512.2%安卡拉64076219.1%伊兹密尔50559818.4%布尔萨42051221.9%安塔利亚31040229.7%科尼亚29533212.5%阿达纳27031817.8%这张表直接看总量不太能说明问题因为每个城市的初始建成区规模差异很大。伊斯坦布尔增加185平方千米安塔利亚增加92平方千米绝对量差距明显但安塔利亚的增幅才是最高的。所以做分析时绝对值要配合增幅一块看。4.2 扩张率怎么算才不误导人我一般的做法是算两类指标年度新增面积和年扩张率。年度新增面积 当年面积 - 上年面积年扩张率 (当年面积 - 上年面积) / 上年面积这两类指标放在折线图里会很有信息量。比如安塔利亚如果旅游地产在某个时间段集中开发对应年份的新增面积会出现一个明显的高峰而不是均匀增长。伊斯坦布尔这种基数大的城市面积虽然增加很多但年扩张率可能反而低于中等城市。这类时序图表可以直接在GEE里用ui.Chart.feature.byFeature画也可以把CSV拿到Excel或Python里处理看个人习惯。我建议在GEE里先跑出来看一遍确认没有异常年份再导出到外部做发表级图表。4.3 判断往哪个方向扩张的进阶思路面积统计只能回答“扩张了多少”回答不了“朝哪个方向扩张”。要回答后者还需要做方向扇区分析。做法是以2018年的城市质心为圆心把城市周围按45度分成8个扇区统计每个扇区内新增城市像元的面积。比如安卡拉如果沿北部新区的交通干线扩展N和NW扇区的新增面积会明显高于其他扇区安塔利亚如果是沿海岸线东西延伸那么E和W扇区会占据主导。这个分析在GEE里也能做生成8个扇形多边形然后和新增城市像元做区域统计。不过我当时是在ArcGIS里完成的这一步因为需要结合道路、地形等矢量数据综合判断桌面GIS会更顺手一些。GEE负责提取和统计桌面GIS负责空间分析这个组合效率其实很高。4.4 结果可信度怎么验证城市扩张监测如果只跑完算法就写结论审稿人或合作方大概率会问一句你提取的建成区边界到底准不准我常用的验证方法有三类。第一随机点目视验证。每个城市每年随机抽100个点叠加到当年的高分辨率影像上人工判断每个点是城市还是非城市然后算混淆矩阵。第二跨尺度验证。用Sentinel-2 10米数据对某一年做同样的提取和Landsat 30米结果对比看空间格局是否一致。第三外部数据交叉验证。VIIRS夜间灯光月度合成数据NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG是个不错的选择建成区通常有稳定的夜间灯光信号而裸土和农田没有。这些验证步骤不是可选项尤其是多城市对比研究里每个城市的光谱环境不一样没有统一的绝对真相只有通过多重证据交叉验证面积统计才站得住脚。5. 实操中的坑与调试经验5.1 水体干扰海峡和海湾边的“假城市”我第一次跑伊斯坦布尔的结果时发现博斯普鲁斯海峡沿岸和海湾区出现了一大片新增城市像元明显不对。水体的光谱特征在某些情况下会落到NDBI大于NDVI的区间特别是含沙量高的浅水区或城市近岸的混合像元不遮水体的话这些都会变成“新增城市”。解决办法就是前面代码里的MNDWI掩膜。实际调试时要注意阈值MNDWI大于0掩膜水体是个常见取值但在部分浑浊水域可能要调到0.05甚至0.1。建议在某个城市区域先打印MNDWI的直方图看看水体峰值在哪里再根据直方图谷底定阈值不要一套参数打天下。5.2 裸土干扰科尼亚大平原上的误判科尼亚是安纳托利亚高原上的农业大市周边有大量半干旱农田和裸露土地。裸土在SWIR1波段反射很强在近红外波段也不算低所以NDBI经常会大于NDVI算法会把大片农田错判成建成区。我处理这个问题时试过两种方式。一种是给每个城市单独调阈值例如把判定规则改成NDBI - NDVI 0.05甚至更高。这个偏移量不能拍脑袋定要对着高分辨率影像看效果。另一种是统计每个城市NDBI减NDVI的直方图正常情况下城市和非城市会形成双峰谷底位置就是该城市的合适阈值。如果双峰不明显说明光谱区分度本身就差这时候需要叠加夜间灯光数据做人机交互筛选而非单纯依赖光谱指数。5.3 面积统计时坐标系和尺度必须显式指定很多GEE新手会在reduceRegion里不写scale和crs直接跑结果和预期差了十万八千里。scale不指定的话GEE会使用影像的默认本机分辨率但不同影像、不同波段可能不一样。crs不指定的话会自动用reduce时的默认投影通常不是面积保持投影统计出的面积可能偏大或偏小。我在这套流程里统一写的scale: 30因为Landsat 8的SR波段本征分辨率就是30米。投影方面土耳其覆盖UTM 35、36、37三个区带若想严格精确可以按城市所在区带分别设置crs到对应的UTM投影。不想那么麻烦的话统一用EPSG:32636对面积统计的影响也很小但至少比默认的Web Mercator靠谱得多。5.4 多云年份合成质量下降怎么办2018到2025年间某些年份、某些城市在夏季云量特别大可用影像数量很少中值合成后局部区域会出现明显的NULL值。这些NULL像素不会被算进建成区面积里导致当年统计面积偏小形成时间序列上的假凹陷。排查方法很简单把每年的合成影像先Map.addLayer叠到地图上看一眼确认城市区域没有明显空洞。如果有空洞解决思路有三个。第一是扩展合成窗口把5到10月的数据都纳进来代价是生长季边缘的NDVI值可能会有偏差。第二是合并Landsat 9数据Landsat 9于2021年底起运行波段设置和Landsat 8基本一致可以补强2022年之后的可用样本。第三对于缺失特别严重的年份可以直接用前后两年的结果做插值但在报告里要标注清楚。5.5 千万别忽视波段选择的一致性GEE的波段名在不同产品集合里可能变化。Landsat C2 L2的波段名是SR_B1到SR_B7不是旧版LANDSAT/LC08/C01/T1_SR里的B1、B2。如果你在代码里写的是select(B4)不会直接报错但选出来的是原始TOA反射率波段不是SR产品缩放系数也对不上结果当然是错的。检查方法很直接在算出NDVI后随便选一个城市的一年输出几个点的值和USGS官网上下载的同日期影像手动算一遍对比。数值如果对不上优先检查波段名和缩放系数。其实这套流程跑稳定之后我后来又把它扩展到了巴尔干半岛的几个国家发现真正可复用的部分不是代码本身而是“先明确实验设计、再写算法、最后用多种方式验证”这个工作思路。城市扩张监测里提取得再漂亮没有可靠的面积统计和空间验证支撑都只是地图上的视觉产物而已。如果你正打算用GEE做类似的长时序城市分析可以把上面的代码块当作起点但一定记得针对自己的研究区和年份重新校准阈值。数据是公开的算法也不复杂真正决定结果质量的永远是你对影像和误差的理解。