1. 辐射定标这件事到底在定什么Landsat 8 的 Level-1 数据产品不管是 Collection 1 还是现在的 Collection 2打开之后你会发现它的像元值DN值是一个 0 到 65535 之间的整数。这个数字本身没有物理意义它只是传感器记录下来的一个相对响应强度。你要拿它做植被指数、地表温度反演、水体提取第一步就必须把这个无量纲的整数转换成有物理量纲的辐射亮度或者反射率。这个过程就是辐射定标。很多人第一次接触遥感影像处理看到“辐射定标”四个字觉得特别高深其实用一句话就能说清楚把传感器记录的原始数字换算成真实世界里的物理量。就像你用手机拍了一张照片照片里某个像素的亮度值是 230这个 230 是相机内部电路给你的一个编码值它不代表那个物体实际反射了多少光。辐射定标就是找到这个编码值和实际物理量之间的换算关系。Landsat 8 上面搭载了两个主要传感器OLIOperational Land Imager陆地成像仪和 TIRSThermal Infrared Sensor热红外传感器。OLI 有 9 个波段覆盖可见光到短波红外TIRS 有 2 个热红外波段。这两个传感器的定标参数是分开的处理的时候要区别对待。OLI 的定标通常输出的是大气表观反射率TOA ReflectanceTIRS 的定标输出的是辐射亮度Radiance或者亮温Brightness Temperature。在 Google Earth Engine后面简称 GEE出现之前做这一步要么用 ENVI 手动跑要么用 IDL 写脚本批处理要么用 Python 的 rasterio 加 numpy 自己算。每种方式都有各自的麻烦ENVI 操作繁琐且不好批量IDL 学习曲线陡峭Python 方案虽然灵活但你要自己管理元数据文件。GEE 的好处在于它已经把 Landsat 8 的定标系数作为影像属性存好了你只需要调用几个内置函数就能在云端完成整个定标流程不用下载原始数据不用装任何软件。这篇文章面向的是已经接触过 GEE 基础操作、但对辐射定标的具体实现还不太有把握的人。我会把 OLI 和 TIRS 的定标逻辑分别拆开讲给出可以直接复制运行的代码并且把我在实际项目中踩过的坑和总结的经验一并分享出来。不管你是做土地利用分类、水体监测还是城市热岛分析这套定标流程都是绕不过去的基础环节。2. 定标之前必须搞清楚的几个概念2.1 DN值、辐射亮度、反射率三者的关系Landsat 8 的原始影像像元值叫 DNDigital Number是一个整数。这个整数经过第一步换算变成辐射亮度Radiance单位是 W/(m²·sr·μm)。辐射亮度再经过第二步换算变成大气表观反射率TOA Reflectance这是一个 0 到 1 之间的无量纲比值。为什么要分两步因为辐射亮度是传感器入瞳处的实际辐射能量它受到太阳高度角、日地距离、传感器增益等多个因素影响。而反射率是把这些影响因素归一化之后的结果不同时间、不同地点拍的影像可以直接对比。做时间序列分析的时候你几乎一定要用反射率而不是辐射亮度。对于热红外波段TIRS 的 B10 和 B11情况稍微不同。热红外波段的物理意义是地表发射的热辐射通常转换成亮温Brightness Temperature单位是开尔文。亮温还不是真正的地表温度它还需要做大气校正和比辐射率校正但那是下一步的事了。辐射定标阶段你至少要把 DN 值转成亮温。2.2 GEE 中 Landsat 8 数据的存储结构GEE 里的 Landsat 8 影像集Collection 2 的命名规则是LANDSAT/LC08/C02/T1_L2和LANDSAT/LC08/C02/T1_TOA。注意这里的区别T1_L2是 Level-2 产品已经做了大气校正直接就是地表反射率T1_TOA是 Level-1 产品只做了辐射定标和几何校正输出的是大气表观反射率。如果你要自己走一遍辐射定标流程应该选T1_TOA这个数据集。虽然它名字里带 TOA但它的 DN 值仍然是原始的量化值你需要自己用定标系数去转换。GEE 已经把定标系数RADIANCE_MULT_BAND_x、RADIANCE_ADD_BAND_x、REFLECTANCE_MULT_BAND_x、REFLECTANCE_ADD_BAND_x、K1_CONSTANT_BAND_x、K2_CONSTANT_BAND_x作为影像属性存好了你直接读取就行。这里有一个容易混淆的点Collection 1 和 Collection 2 的定标系数存储方式略有不同。Collection 1 的系数是直接写在影像属性里的Collection 2 也是但属性名称的格式有细微差别。如果你在网上找到的代码跑不通很可能就是 Collection 版本不匹配导致的。2.3 太阳高度角和日地距离修正做反射率定标的时候有两个修正因子必须考虑太阳高度角Sun Elevation和日地距离Earth-Sun Distance。太阳高度角影响的是太阳光到达地表的入射角度。如果太阳高度角是 60 度那太阳光不是垂直照射地面的实际到达地表的能量要除以 sin(60°)。GEE 里每景影像都有一个SUN_ELEVATION属性单位是度。日地距离是因为地球绕太阳的轨道是椭圆形的不同季节日地距离不一样到达大气层顶的太阳辐射强度也不一样。这个修正因子通常用ESUN大气层顶太阳辐照度和d日地距离天文单位来计算。GEE 的TOA数据集里已经内置了这些参数你调用ee.Algorithms.Landsat.TOA()的时候它会自动处理。但如果你要手动一步步算就得自己把这些因子乘进去。手动算的好处是你能完全掌控每一步出了问题也知道去哪里排查。自动算的好处是省事但万一结果不对你很难定位是哪个环节出了偏差。3. OLI 传感器辐射定标全流程拆解3.1 从 DN 值到辐射亮度OLI 的辐射亮度计算公式是L M * DN A其中L是辐射亮度M是RADIANCE_MULT_BAND_xA是RADIANCE_ADD_BAND_x。这两个系数在 GEE 的影像属性里可以直接读到。在 GEE 里实现这一步代码大概长这样var radiance function(image) { var mult image.get(RADIANCE_MULT_BAND_4); var add image.get(RADIANCE_ADD_BAND_4); var rad image.select(B4).multiply(mult).add(add); return rad.rename(B4_radiance); };这段代码看起来简单但有几个细节要注意。第一get()返回的是ee.Number类型直接和影像做乘法运算时 GEE 会自动广播但如果你要做更复杂的运算最好先ee.Number()显式转换一下。第二不同波段的系数是不一样的你不能用一个系数套所有波段。第三Collection 2 的系数值可能和 Collection 1 不同因为传感器的响应会随时间衰减定标系数是定期更新的。我实测下来B4红光波段的RADIANCE_MULT_BAND_4大约是 0.012RADIANCE_ADD_BAND_4大约是 -63.6。这个负的加系数是因为 DN 值 0 对应的辐射亮度不是 0而是经过偏移的。Collection 2 引入了这个偏移量目的是让 DN 值的动态范围更大避免高亮区域饱和。3.2 从辐射亮度到大气表观反射率反射率的计算公式是ρ (π * L * d²) / (ESUN * cos(θ))其中ρ是反射率L是上一步算出来的辐射亮度d是日地距离天文单位ESUN是大气层顶太阳辐照度θ是太阳天顶角90° 减去太阳高度角。GEE 提供了一个更简便的公式直接用 DN 值算反射率ρ (M * DN A) / sin(θ)其中M是REFLECTANCE_MULT_BAND_xA是REFLECTANCE_ADD_BAND_xθ是太阳高度角。这个公式把日地距离修正和 ESUN 都合并到系数里了用起来更方便。代码实现var reflectance function(image) { var mult image.get(REFLECTANCE_MULT_BAND_4); var add image.get(REFLECTANCE_ADD_BAND_4); var sunElev image.get(SUN_ELEVATION); var ref image.select(B4).multiply(mult).add(add) .divide(Math.sin(sunElev * Math.PI / 180)); return ref.rename(B4_reflectance); };这里最容易出错的地方是太阳高度角的单位。GEE 里SUN_ELEVATION是度数但 JavaScript 的Math.sin()函数接受的是弧度。你必须先转换成弧度否则算出来的反射率会完全不对。我见过不止一个人在这里栽跟头算出来的反射率要么是负数要么大得离谱。还有一个坑是反射率算出来之后可能会有一些像元的值小于 0 或者大于 1。小于 0 通常是因为大气散射或者传感器噪声大于 1 可能是因为云层或者水面镜面反射。这些异常值在后续分析中需要处理但定标阶段不用管保留原值就行。3.3 批量处理整景影像的所有波段实际项目中你不会只处理一个波段而是要把 OLI 的 9 个波段全部定标。手动一个波段一个波段写太累了用循环或者map()函数批量处理更高效。var oliBands [B1,B2,B3,B4,B5,B6,B7,B8,B9]; var calibrateOLI function(image) { var calibrated oliBands.map(function(band) { var mult image.get(REFLECTANCE_MULT_ band); var add image.get(REFLECTANCE_ADD_ band); var sunElev image.get(SUN_ELEVATION); var ref image.select(band).multiply(mult).add(add) .divide(Math.sin(sunElev * Math.PI / 180)); return ref.rename(band _ref); }); return image.addBands(calibrated, null, true); };这段代码里addBands(calibrated, null, true)的第三个参数true表示覆盖同名的原始波段。如果你不想覆盖改成false就会保留原始波段新增定标后的波段。我个人的习惯是覆盖因为原始 DN 值在后续分析中基本用不到留着反而占内存。批量处理的时候要注意波段名称的对应关系。REFLECTANCE_MULT_B1对应的是 B1 波段REFLECTANCE_MULT_B2对应的是 B2 波段以此类推。但 B8 是 panchromatic 波段它的定标系数和别的波段不一样分辨率也不同15米 vs 30米。如果你要把 B8 和其他波段一起处理后面做波段运算的时候需要先重采样到统一分辨率。4. TIRS 传感器辐射定标与亮温计算4.1 热红外波段的特殊性TIRS 的两个波段 B1010.60-11.19 μm和 B1111.50-12.51 μm和 OLI 的波段有本质区别。OLI 测的是反射的太阳辐射TIRS 测的是地表发射的热辐射。所以 TIRS 的定标流程不是算反射率而是算亮温。亮温的计算分两步先把 DN 值转成辐射亮度再把辐射亮度转成亮温。辐射亮度的公式和 OLI 一样L M * DN A但亮温的公式是普朗克函数的反函数T K2 / ln(K1 / L 1)其中K1和K2是定标常数存在影像属性K1_CONSTANT_BAND_10和K2_CONSTANT_BAND_10里。L是上一步算出来的辐射亮度。算出来的T单位是开尔文。4.2 GEE 中亮温计算的代码实现var brightnessTemp function(image) { var radMult image.get(RADIANCE_MULT_BAND_10); var radAdd image.get(RADIANCE_ADD_BAND_10); var k1 image.get(K1_CONSTANT_BAND_10); var k2 image.get(K2_CONSTANT_BAND_10); var radiance image.select(B10).multiply(radMult).add(radAdd); var temp radiance.multiply(0).add(k2).divide( radiance.divide(k1).add(1).log() ); return temp.rename(B10_temp); };这段代码里radiance.multiply(0).add(k2)是一个小技巧目的是创建一个和radiance同样大小的影像每个像元的值都是k2。GEE 里影像和数字做运算时数字会自动广播到所有像元但如果你直接写k2.divide(...)k2是一个ee.Number它不能直接和影像做除法。所以要先把它变成影像。log()是自然对数GEE 里还有log10()是常用对数别搞混了。普朗克反函数用的是自然对数。算出来的亮温单位是开尔文如果你想要摄氏度减去 273.15 就行。但做地表温度反演的时候通常保留开尔文因为后续的比辐射率校正和大气校正公式都是用开尔文推导的。4.3 B10 和 B11 的差异处理Landsat 8 的 B11 波段从 2013 年发射以来就一直有问题它的定标精度不如 B10而且存在杂散光stray light影响。USGS 官方建议在做地表温度反演的时候只用 B10B11 作为辅助或者干脆不用。在 GEE 里处理的时候如果你要做分裂窗算法split-window algorithm需要同时用 B10 和 B11。但如果你只是做单通道算法用 B10 就够了。我个人的经验是如果研究区不是特别大直接用 B10 做单通道反演结果更稳定也省去了处理 B11 杂散光的麻烦。5. 用 GEE 内置函数一键定标5.1 ee.Algorithms.Landsat.TOA() 的用法和局限GEE 提供了一个内置函数ee.Algorithms.Landsat.TOA()可以直接把 Landsat 8 的 DN 值转成 TOA 反射率或亮温。用法很简单var toa ee.Algorithms.Landsat.TOA(image);这个函数会自动判断哪些波段是反射波段哪些是热红外波段然后分别处理。反射波段输出反射率热红外波段输出亮温。你不需要手动读定标系数也不需要管太阳高度角修正。但它有几个局限。第一它只适用于 Collection 1 的 TOA 数据Collection 2 的数据用这个函数可能会报错或者结果不对。第二它不输出辐射亮度如果你需要辐射亮度作为中间结果得自己算。第三它不处理 B8 全色波段因为 B8 的定标逻辑和别的波段不一样。我实测下来这个函数在处理 Collection 1 数据时很稳但 Collection 2 数据最好还是手动定标。USGS 从 2022 年开始全面推行 Collection 2新项目建议直接用手动方案。5.2 手动定标 vs 自动定标的取舍手动定标的好处是透明、可控、可调试。每一步你都知道自己在算什么出了问题能定位到具体环节。缺点是代码量大容易写错。自动定标的好处是省事一行代码搞定。缺点是黑箱你不知道它内部做了什么处理万一结果不对很难排查。我的建议是如果你只是做快速预览或者教学演示用自动定标没问题。但如果是正式的研究项目或者生产环境一定要用手动定标。因为手动定标的过程中你会对数据的物理意义有更深的理解而且可以根据具体需求调整参数。还有一个实际考虑GEE 的免费计算资源是有限的。自动定标函数内部做了很多额外的检查和处理计算量比手动定标大。如果你要处理几百景影像手动定标能省不少计算时间。6. 定标结果验证与常见问题排查6.1 如何判断定标结果是否合理定标做完之后怎么知道结果对不对我通常用以下几个方法交叉验证。第一看数值范围。OLI 的 TOA 反射率正常范围是 0 到 1大部分地物的反射率在 0.05 到 0.6 之间。植被在近红外波段B5反射率可以到 0.5 以上水体在近红外波段反射率接近 0。如果你算出来的反射率大面积超过 1 或者全是负数那肯定有问题。第二看光谱曲线。选一个植被像元画出它的光谱曲线应该是红光波段B4吸收、近红外波段B5反射的典型植被光谱。如果曲线形状不对说明定标系数用错了。第三和 Level-2 产品对比。GEE 里有LANDSAT/LC08/C02/T1_L2数据集这是官方做好的地表反射率产品。你可以选同一景影像把你的 TOA 反射率和 L2 的地表反射率对比。两者不会完全一样因为 L2 做了大气校正但趋势应该一致。如果差异特别大说明你的定标流程有问题。6.2 常见报错与解决方法问题现象可能原因解决方法反射率全为负数太阳高度角未转弧度检查Math.sin()参数是否乘以 π/180反射率大于 10定标系数读取错误确认波段名称和系数属性名匹配亮温计算报错K1/K2 常数未读取到检查影像是否有K1_CONSTANT_BAND_10属性波段名称冲突addBands 覆盖了原始波段设置第三个参数为 false 或重命名计算结果全为 0影像被过滤掉了检查 filter 条件是否过严6.3 几个我踩过的坑第一个坑Collection 2 的 DN 值有偏移。Collection 1 的 DN 值范围是 0 到 255对于 8 位数据或者 0 到 65535对于 16 位数据但 Collection 2 引入了 -0.2 的反射率偏移和 -63.6 的辐射亮度偏移。如果你用 Collection 1 的公式去算 Collection 2 的数据结果会偏。一定要确认你用的系数和数据类型匹配。第二个坑太阳高度角在极地地区可能很小。如果你研究的是高纬度地区冬季的太阳高度角可能只有几度sin(θ)接近 0算出来的反射率会非常大。这种情况下要么换季节要么对反射率做截断处理。第三个坑B8 全色波段的分辨率是 15 米其他波段是 30 米。如果你把 B8 和其他波段一起做定标然后叠加GEE 会自动重采样但重采样方法默认是最近邻可能会引入误差。如果需要精确的波段运算最好先统一分辨率再定标。第四个坑GEE 的get()函数返回的是ee.Number但如果你在map()函数里用get()每次迭代都会发起一次服务器请求速度很慢。优化方法是先把所有需要的系数一次性读出来存成字典然后在map()里查字典。7. 定标之后的下一步从 TOA 到地表参数辐射定标只是预处理的第一步。定标完成之后你得到的是大气表观反射率或者亮温这些还不是最终的地表参数。如果你要做植被覆盖度分析下一步是算 NDVI然后做像元二分模型。如果你要做地表温度反演下一步是大气校正和比辐射率校正。NDVI 的计算很简单var ndvi image.normalizedDifference([B5_ref, B4_ref]).rename(NDVI);但要注意NDVI 对大气效应很敏感TOA 反射率算出来的 NDVI 和地表反射率算出来的 NDVI 会有差异。如果研究区大气比较干净比如高原地区差异不大如果大气污染严重或者水汽含量高差异可能超过 0.1。这种情况下最好用 Level-2 产品或者自己做大气校正。地表温度反演的单通道算法需要三个输入亮温、比辐射率、大气水汽含量。亮温你已经有了比辐射率可以从 NDVI 估算大气水汽含量可以从再分析数据获取。GEE 里有 MODIS 的大气产品可以拿来用。但这是另一个话题了展开讲又是一篇文章。8. 一些提高效率的实操技巧8.1 用 ee.Dictionary 缓存定标系数前面提到过在map()里反复调用get()会很慢。优化方法是在循环外面把所有系数读出来存成一个字典var coeffs ee.Dictionary({ B1: {mult: image.get(REFLECTANCE_MULT_B1), add: image.get(REFLECTANCE_ADD_B1)}, B2: {mult: image.get(REFLECTANCE_MULT_B2), add: image.get(REFLECTANCE_ADD_B2)}, // ... });然后在map()里用coeffs.get(B1)来取。这样只需要一次服务器请求速度能快好几倍。8.2 用 export 批量导出定标结果GEE 的交互式地图适合预览但不适合大批量处理。如果你要处理几十景甚至上百景影像最好用Export.image.toDrive()批量导出到云盘然后在本地做后续分析。导出的时候要注意设置scale参数。Landsat 8 的 OLI 波段是 30 米分辨率TIRS 是 100 米重采样到 30 米。如果你导出的时候不指定scaleGEE 会用默认值可能导致分辨率不一致。Export.image.toDrive({ image: calibratedImage, description: L8_calibrated, scale: 30, region: studyArea, maxPixels: 1e13 });maxPixels要设大一点否则大区域导出会报错。1e13是我常用的值一般够用。8.3 用 QA 波段做云掩膜定标完成之后下一步通常是去云。Landsat 8 的 Collection 2 产品自带 QA_PIXEL 波段里面用比特位标记了云、云阴影、水体等信息。你可以用这个波段做云掩膜var cloudMask function(image) { var qa image.select(QA_PIXEL); var cloud qa.bitwiseAnd(1 3).eq(0); var shadow qa.bitwiseAnd(1 4).eq(0); return image.updateMask(cloud.and(shadow)); };1 3是比特位运算表示第 3 位是云标记。Collection 2 的 QA 波段比特位定义和 Collection 1 不同用之前一定要查官方文档确认。8.4 定标系数的时效性Landsat 8 的定标系数不是一成不变的。USGS 会定期根据传感器状态更新系数通常每年更新一次。GEE 里的数据会自动同步最新的系数所以你不用手动更新。但如果你在本地用 Python 处理下载下来的原始数据就要注意检查系数文件的版本。我遇到过一种情况用两年前下载的系数处理今年的影像结果反射率整体偏高了 5% 左右。虽然差异不大但对于做定量分析的项目来说5% 的偏差可能影响结论。所以本地处理的时候一定要用和影像同期发布的系数文件。9. 不同场景下的定标策略选择9.1 时间序列分析做时间序列分析的时候最重要的是不同时相影像之间的可比性。这时候你必须用反射率而不是辐射亮度因为辐射亮度受太阳高度角和日地距离影响不同季节的影像没法直接对比。另外时间序列分析通常需要多景影像合成。合成之前要确保所有影像都经过了相同的定标流程用了相同的系数版本。如果中间换了系数版本合成结果会出现突变。9.2 单景高精度分析如果只处理一景影像而且对精度要求很高可以考虑用辐射亮度而不是反射率。因为反射率的计算引入了太阳高度角修正如果太阳高度角数据不准确反而会引入误差。辐射亮度是传感器直接测量的物理量不依赖外部参数。但辐射亮度的缺点是不同影像之间没法直接对比。所以如果你的分析不涉及多时相对比用辐射亮度也没问题。9.3 热红外应用热红外应用比如城市热岛、火山监测、水温反演必须用亮温或者地表温度。定标阶段算出来的亮温已经可以用于很多定性分析了比如找热异常区域。但如果要做定量的温度反演还需要做大气校正和比辐射率校正。大气校正需要大气水汽含量、气温、气压等参数。这些参数可以从气象再分析数据获取也可以用 MODIS 的大气产品。GEE 里有 NCEP/NCAR 再分析数据可以拿来用。但大气校正的精度很大程度上取决于大气参数的精度如果研究区地形复杂或者天气多变校正结果可能不理想。10. 关于定标精度的一点个人体会辐射定标的精度直接影响后续所有分析的可靠性。我见过很多项目前面分类、反演做得花里胡哨但定标这一步没做好结果整个结论都站不住脚。定标精度主要受三个因素影响定标系数本身的精度、太阳高度角数据的精度、以及计算过程中的数值精度。定标系数是 USGS 提供的你没法改变太阳高度角是影像自带的也没法改变你能控制的只有计算过程。计算过程中最容易损失精度的地方是数据类型转换。GEE 默认用浮点数存储影像但如果你在计算过程中不小心把结果转成了整数小数部分就被截断了。比如image.select(B4).multiply(mult)如果mult是整数结果就是整数。所以一定要确保所有参与运算的系数都是浮点数。另一个精度问题是log()函数的计算。亮温公式里有ln(K1/L 1)当L很小时K1/L很大ln的结果对L的微小变化很敏感。如果L的精度不够亮温的误差会被放大。所以热红外波段的定标一定要用双精度浮点数。最后说一个我自己的习惯每次定标完成之后我都会随机选几个像元手动用计算器验算一遍。虽然麻烦但能确保代码逻辑没问题。这个习惯帮我发现过好几次系数读取错误的问题。
