拓冰网站定制 
这些年做生态环境遥感评估我周围很多同行一开始都习惯用NDVI单指标来看一个区域的生态状况但说实话单看绿度太片面了。一个城市周边可能有大量农田NDVI很高但地表温度、建筑裸土占比这些生态压力因素完全没体现出来。后来接触到徐涵秋老师提出的遥感生态指数RSEI模型才真正觉得找到了一个既科学严谨、又能快速出结果的多指标评价框架。这篇文章我就把这套方法从原理到实操完整梳理一遍给正在做生态环境评价、城市规划研究、或者毕业论文需要出生态质量图的朋友做个参考。RSEI的核心思路其实很直接把影响生态环境的四大要素——绿度、湿度、热度、干度——用遥感数据分别量化再用主成分分析把它们客观地压缩成一个综合指数。整个流程在ENVI、ArcGIS或者Google Earth Engine上都能跑通关键是理解每一步为什么要这么做。1. RSEI模型到底在评价什么四大分量的生态语义在动手算之前我们得先把RSEI的四个输入分量吃透。这决定了你后面选数据、调参数的时候心里有没有底。1.1 绿度植被覆盖的健康表绿度指标用的是NDVI即归一化植被指数。它的原理是利用植被在红光波段强吸收、近红外波段强反射的光谱特性把植被覆盖信息从地表信号里分离出来。公式是NDVI (NIR - Red) / (NIR Red)对Landsat 8/9来说就是Band5 - Band4/Band5 Band4。NDVI的取值范围在-1到1之间理论上裸土接近0茂密森林可以到0.8以上水体为负值。我在实际处理中最想提醒的一点是NDVI虽然看起来简单但不同季节的植被状态差异极大。如果是做多时期对比影像的获取月份必须尽量一致比如都用每年8-9月份的影像否则你很难区分生态变化和物候变化。1.2 湿度土壤含水与地表湿润度湿度分量来自缨帽变换Tasseled Cap Transformation的湿度分量WET。缨帽变换是一种针对植被、土壤、水体等典型地物设计的线性变换它能把原始多波段数据压缩成几个有明确生态含义的组分亮度、绿度和湿度。不同传感器的WET系数差别很大这点非常容易踩坑。Landsat 8 OLI地表反射率产品的WET系数是WET 0.1511 * Blue 0.1973 * Green 0.3283 * Red 0.3407 * NIR - 0.7117 * SWIR1 - 0.4559 * SWIR2注意这里的系数是基于地表反射率数据推导的如果你用的是大气表观反射率TOA产品系数必须换一套。另外有些教程里给的是Landsat 5 TM的系数直接套到Landsat 8上结果出来全是乱的这个坑后面我会详细讲。1.3 热度地表温度的生态压力热度分量就是地表温度LST。城市热岛效应是生态质量恶化的最直接表现之一所以在RSEI里热度被当作重要压力指标。LST的反演方法有很多种比如辐射传输方程法、劈窗算法、单窗算法等。实操中我们最常用的是辐射传输方程法大气校正法它的基本思路是先把热红外波段的DN值转换为辐射亮度再转换为亮温然后结合地表比辐射率做发射率校正最终得到真实的地表温度。这个指标的计算链路最长涉及到好几个中间参数很多新手就是在这步被劝退的。别急第3章我会把每一步的推导和代码都给你列出来。1.4 干度建筑与裸土的人工印记干度指标NDBSINormalized Difference Bare Soil and Building Index是裸土指数SI和建筑指数IBI的均值NDBSI (SI IBI) / 2为什么要用两个指数的均值因为单一指数很难同时兼顾自然裸土和人工建筑。比如传统的裸土指数容易把建筑也识别进来建筑指数又容易漏掉裸土。把两者平均能在一定程度上相互弥补将这两类不透水面和裸露地表统一表征为生态干度。这个指标是RSEI里最能反映人类活动干扰强度的分量。一个区域如果NDBSI高通常意味着城市化程度高、地表硬化严重或者水土流失明显这些都属于生态退化的重要信号。理解了这四大分量的生态含义之后你会发现RSEI其实构建了一个压力-状态-响应的简化框架植被状态绿度、湿度代表生态系统的基础状况人类活动干扰热度、干度代表生态压力。综合起来就得到了一个能同时反映自然禀赋和人为影响的生态质量指数。2. 数据源选型与预处理决定RSEI成败的第一道关口方向不对努力白费。数据选型和预处理的质量直接决定了后面所有计算结果的可靠性。2.1 Landsat 8/9还是Sentinel-2怎么选目前做RSEI用得最多的是Landsat系列主要原因是它拥有从1980年代至今的连续存档非常适合做几十年的生态变化分析。具体到传感器我推荐以下选型原则数据源分辨率优势注意点Landsat 8/9 OLITIRS30m热红外100m2013年至今存档连续有成熟的热红外波段重访周期16天受云影响大Landsat 5 TM30m热红外120m1984-2011年历史数据数据较老辐射定标参数需仔细处理Landsat 7 ETM30m1999年至今2003年后有条带需要插值修复Sentinel-2 MSI10/20m分辨率高重访周期5天没有热红外波段LST需要单独处理有一点必须澄清Sentinel-2本身没有热红外传感器直接用Sentinel-2做完整的RSEI是算不出LST的。虽然有些研究用Landsat的LST降尺度配合Sentinel-2做高分辨率RSEI但那超出了本文基础教程范畴。我建议新手先老老实实用Landsat 8/9的Collection 2 Level 2产品。2.2 云掩膜与影像合成云是光学遥感最大的敌人。一片薄云或者云阴影覆盖的区域反射率信号完全失真如果没做掩膜就参与计算那几个像元的NDVI、LST都会变成离谱的异常值直接影响后续PCA的特征向量。在Google Earth Engine里我们可以直接用Landsat Collection 2 Level 2产品自带的QA_PIXEL波段做云掩膜这是目前最省事也最靠谱的做法。我常用的JavaScript代码如下function maskL8sr(image) { // 获取QA波段 var qa image.select(QA_PIXEL); // 设置要掩膜的位标志云、云阴影、冰雪 var cloudBitMask 1 1; // Cloud var shadowBitMask 1 3; // Cloud Shadow var snowBitMask 1 4; // Snow // 创建掩膜 var mask qa.bitwiseAnd(cloudBitMask).eq(0) .and(qa.bitwiseAnd(shadowBitMask).eq(0)) .and(qa.bitwiseAnd(snowBitMask).eq(0)); // 去掉缩放系数保留SR波段和热红外波段 return image.updateMask(mask) .select([SR_B2,SR_B3,SR_B4,SR_B5,SR_B6,SR_B7,ST_B10], [Blue,Green,Red,NIR,SWIR1,SWIR2,Thermal]); }这里要特别注意热红外波段ST_B10的定标单位。Collection 2 Level 2的地表温度产品单位是开尔文乘以10在后续计算时需要先乘以0.00341802并加上149.0来还原成真实的开尔文温度这个细节很坑后面演示代码里我会带上。如果你用的是ENVI/ArcGIS路线那么下载数据时应优先选Landsat Collection 2 Level 2科学产品而不是Level 1原始DN产品。Level 2已经帮你完成了大气校正可以直接拿SR波段做指数计算省掉了最耗时的大气校正环节。2.3 投影、裁剪与统一分辨率RSEI的四个分量中NDVI、WET、NDBSI是反射率比值或线性组合不受投影影响但LST是热红外波段计算来的原始分辨率是100米。如果在ENVI里直接做波段合成需要先把热红外重采样到30米否则最后PCA的输入图层分辨率不一致会产生错位。在GEE中处理非常简单用select和resample即可// 重采样到30米并统一投影 var thermal image.select(Thermal).resample(bilinear); var sortedImage image.addBands(thermal, null, true);统一投影时建议全部转到UTM投影。如果研究区跨多个UTM带可以用Albers等积投影保证面积计算不变形。另外别忘了裁剪到研究区边界。用clip接研究区的FeatureCollection即可var roi ee.FeatureCollection(projects/your-project/assets/roi); var imageClipped sortedImage.clip(roi);这一步能显著减少后续PCA的计算量尤其是研究区只占影像一小部分的时候。3. 四大指标计算实操公式、参数与GEE代码预处理做完就到了最核心的指标计算环节。我直接以Landsat 8/9 Collection 2 Level 2为输入给出完整的GEE实现。3.1 NDVI与WET的计算NDVI前面已经给了公式在GEE里用normalizedDifference一行就能解决var ndvi image.normalizedDifference([NIR, Red]).rename(NDVI);WET的计算需要用到缨帽变换系数。这里我给出Landsat 8 OLI地表反射率对应的完整公式var wet image.expression( 0.1511 * B2 0.1973 * B3 0.3283 * B4 0.3407 * B5 - 0.7117 * B6 - 0.4559 * B7, { B2: image.select(Blue), B3: image.select(Green), B4: image.select(Red), B5: image.select(NIR), B6: image.select(SWIR1), B7: image.select(SWIR2) } ).rename(WET);我把这一段单独列出来是因为系数搞混的概率实在太高了。网上流传的缨帽变换系数有很多版本有基于TOA反射率的、有基于地表反射率的、有TM的、有ETM的用错一个版本湿度分量就可能出现大面积负数后面PCA的结果也会跟着错。3.2 Landsat 8地表温度LST反演的完整链路LST是RSEI四个分量里计算链路最长的我按顺序把每一步拆开讲清楚。第一步读取热红外波段的DN值并定标。GEE里ST_B10已经是Level 2处理过的地表温度产品但单位是开尔文×10需要还原var thermalDN image.select(Thermal); var lstKelvin thermalDN.multiply(0.00341802).add(149.0).rename(LST_K);如果用的是Collection 2 Level 1的原始数据则需要先做辐射定标得到辐射亮度L再用热红外波段的K1、K2常数换算亮温。Landsat 8 Band 10的K1774.8853、K21321.0789这是当前Collection 2版本的官方常数不同版本有差异务必查询影像MTL文件// 如果从L1数据出发DN转辐射亮度 var ML 0.0003342; // 从MTL文件获取 var AL 0.1; // 从MTL文件获取 var radiance image.select(B10).multiply(ML).add(AL); // 辐射亮度转亮温 var brightnessTemp radiance.expression( K2 / log(K1 / rad 1), {rad: radiance, K1: 774.8853, K2: 1321.0789} );第二步计算地表比辐射率。这一步需要用NDVI估算植被覆盖度FVC再根据植被和裸土的比辐射率加权平均// 植被覆盖度 var ndviMin 0.05; var ndviMax 0.75; var fvc ndvi.subtract(ndviMin).divide(ndviMax - ndviMin).clamp(0, 1).rename(FVC); // 地表比辐射率水体取0.99城镇和裸地取0.97自然地表用FVC加权 var emissivity fvc.expression( (FVC 0) (NDVI 0.05) ? 0.995 : ((NDVI 0.05 NDVI 0.75) ? (0.004 * FVC 0.986) : 0.99), {FVC: fvc, NDVI: ndvi} ).rename(EM);第三步用单窗算法反演真实地表温度。经典的简化公式是LST BT / [1 (λ * BT / ρ) * ln(ε)]其中λ是热红外波段中心波长Landsat 8 Band 10取10.895微米ρ h * c / σ ≈ 1.438 * 10^-2 m·Kε为比辐射率。代码实现如下var lambda 10.895; // 微米 var rho 1.438e-2; // 米·开尔文 var lst brightnessTemp .divide(brightnessTemp.multiply(lambda).divide(rho).log().multiply(-1).add(1))等等这个公式我写反了。正确的单窗算法简化式是LST BT / (1 (λ * BT / ρ) * ln(ε))因为对数项是负值ε 1所以分母小于1LST会比亮温略高。GEE里正确写法是var lst brightnessTemp.divide( ee.Image(1).add( brightnessTemp.multiply(lambda).divide(rho).multiply(emissivity.log()) ) ).rename(LST);最后别忘了把开尔文转摄氏度方便制图时标注var lstC lst.subtract(273.15).rename(LST_C);3.3 NDBSI建筑与裸土指数计算NDBSI包含SI和IBI两个指数先分别算再取平均。SI的公式为SI [(SWIR1 Red) - (NIR Blue)] / [(SWIR1 Red) (NIR Blue)]IBI的公式为IBI {2SWIR1/(SWIR1NIR) - [NIR/(NIRRed) Green/(GreenSWIR1)]} / {2SWIR1/(SWIR1NIR) [NIR/(NIRRed) Green/(GreenSWIR1)]}从结构上可以看出IBI利用的是建筑在短波红外波段的反射率明显高于植被的特点以及植被在近红外高反射、建筑在绿光波段反射率相对高的光谱差异。代码实现如下var si image.expression( ((SWIR1 Red) - (NIR Blue)) / ((SWIR1 Red) (NIR Blue)), { SWIR1: image.select(SWIR1), Red: image.select(Red), NIR: image.select(NIR), Blue: image.select(Blue) } ).rename(SI); var ibi image.expression( (2*SWIR1/(SWIR1NIR) - (NIR/(NIRRed) Green/(GreenSWIR1))) / (2*SWIR1/(SWIR1NIR) (NIR/(NIRRed) Green/(GreenSWIR1))), { SWIR1: image.select(SWIR1), NIR: image.select(NIR), Red: image.select(Red), Green: image.select(Green) } ).rename(IBI); var ndbsi si.add(ibi).multiply(0.5).rename(NDBSI);到这里NDVI、WET、LST、NDBSI四个分量影像就全部准备好了。接下来的PCA才是RSEI的灵魂。4. 主成分分析PCA与RSEI合成从四个指标到一个指数为什么要用PCA而不是简单地给四个指标各赋一个权重求平均这是RSEI模型最有技术含量的地方。4.1 为什么先做归一化四个分量的量纲和数值范围差异巨大NDVI大致在-1到1之间LST可能在280到320开尔文WET的范围因传感器而异NDBSI也在-1到1附近。如果不归一化直接做PCA数值范围大的LST会主导主成分的方差贡献导致特征向量严重偏向热度分量其余分量的信息被压制。归一化公式是每个指标按0-1区间拉伸NI (I - I_min) / (I_max - I_min)在GEE里实现需要注意这里的I_min和I_max应该基于研究区影像的实际像元分布取而不是理论值。但直接取影像的最小最大值容易被异常像元干扰。我习惯用2%到98%的分位数来截断归一化这样能有效抑制极少数云残留或水体异常值的影响function normalize(img) { var min img.reduceRegion({ reducer: ee.Reducer.percentile([2]), geometry: roi, scale: 30, maxPixels: 1e10 }).values().get(0); var max img.reduceRegion({ reducer: ee.Reducer.percentile([98]), geometry: roi, scale: 30, maxPixels: 1e10 }).values().get(0); return img.subtract(ee.Image.constant(min)).divide(ee.Image.constant(max).subtract(ee.Image.constant(min))) .clamp(0, 1); }这里get(0)取到的值是数组形式在GEE里要小心处理。更稳妥的做法是用reduceRegion返回的字典按波段名分别取min和max。另外注意maxPixels要设得足够大否则大范围研究区会报错。4.2 PCA的输入与主成分贡献率判断归一化完成后将四个分量合成一个多波段影像然后在GEE里用reduceRegion配合ee.Reducer.principalComponentAnalysis()做PCA。PCA的核心逻辑是求四个波段构成的协方差矩阵的特征值和特征向量按特征值大小排序提取第一主成分PC1。var image4bands ndvi.addBands([wet, lst, ndbsi]).select([NDVI,WET,LST,NDBSI]).float(); var pca ee.Reducer.principalComponentAnalysis(4); var pcImage image4bands.reduceNeighborhood(pca, ee.Kernel.square(1));等等reduceNeighborhood不是我们想要的方案。它是在滑动窗口里做PCA计算量和结果的解释方式都不对。正确的做法应该是先在全影像范围内采样构建协方差矩阵再对影像做线性变换// 先采样得到均值与协方差 var mean image4bands.reduceRegion({ reducer: ee.Reducer.mean(), geometry: roi, scale: 30, maxPixels: 1e10 }); var centered image4bands.select([NDVI,WET,LST,NDBSI]) .subtract(ee.Image.constant(mean.values())); var covar centered.toArray().reduceRegion({ reducer: ee.Reducer.covariance(), geometry: roi, scale: 30, maxPixels: 1e10 }); // 从协方差做PCA var eigen ee.Array(covar.get(array)).eigen(); var eigenVectors eigen.slice(1, 0, 1).slice(0, 0, 4); var pc1 centered.toArray().matrixMultiply(eigenVectors.transpose()).arrayProject([0]) .arrayFlatten([[PC1]]);这一段代码有几个技术细节值得展开eigen()返回的是一个二维数组每一行格式是[特征值, 特征向量分量...]按特征值降序排列。所以slice(1, 0, 1)取出的是第一行后面的部分即最大特征值对应的特征向量。这个特征向量就是PC1对应的权重它的四个分量分别代表NDVI、WET、LST、NDBSI在PC1中的载荷方向。拿到了PC1之后我们要看两个关键指标第一是特征值贡献率。PC1对应的特征值占所有特征值之和的比例一般要求大于80%。如果贡献率偏低说明四个分量的信息没有很好地收敛到一个主轴上模型解释力不够。第二是特征向量的符号。在原始的RSEI框架中我们期望NDVI和WET的载荷为正生态好的指标贡献正向LST和NDBSI的载荷为负生态压力指标贡献负向。如果你算出来的PC1载荷方向完全相反即NDVI为负、LST为正说明第一主成分表达的是生态退化而不是生态质量。这种情况下需要把RSEI反过来用1减去PC1来校正方向否则高值代表恶劣、低值代表良好语义全反了。4.3 RSEI合成、再归一化与方向校正无论特征向量方向如何最终都要做一个关键步骤把PC1再归一化到0-1区间然后判断是否需要方向修正。// 先归一化PC1到0-1 var pc1Norm normalize(pc1); // 判断方向并校正 // 如果NDVI、WET在特征向量中的系数之和为正说明PC1越大生态越好直接用pc1Norm // 如果为负则需要用 1 - pc1Norm var rsei pc1Norm; // 或者 rsei ee.Image(1).subtract(pc1Norm);这里我强烈建议你打印出特征向量来人工检查一下别全指望自动判断。具体做法是print(eigenVectors)在Console里查看四分量各自的载荷值。我处理过不少区域大部分时候NDVI和WET的载荷为正但也有例外——比如在大范围荒漠地区PC1可能完全由LST主导这时候就要特别小心。合成之后的RSEI理论上取值范围是0到1数值越高代表生态质量越好。但实际像元分布往往不均匀可能集中在0.4到0.7之间。这时候我会做一次直方图拉伸或分位数拉伸方便后续分级和可视化。5. 结果分级、变化检测与制图输出RSEI算出来是个连续变量但研究上最常用的呈现方式是分级。分级既方便空间对比也好做面积统计。5.1 五级分类阈值的设定目前学界用得最多的是等间距五级分类法把0-1区间均分为五等份RSEI区间等级生态含义0.0 - 0.2差生态质量极差以高强度开发或严重退化区域为主0.2 - 0.4较差生态质量偏差人为干扰明显0.4 - 0.6中等生态质量一般自然与人为活动交织0.6 - 0.8良生态质量较好植被覆盖较高0.8 - 1.0优生态质量优异基本为茂密植被区等间距分类最大的优势是可对比性。不同年份、不同研究区的RSEI结果都能用同一套标准比较。不过也有学者用自然断点法Jenks分级让每个类别的内部差异最小化。我自己做多时期对比时坚持用等间距法因为自然断点法每年的断点不同变化统计的结果会被分类方法本身干扰很难说清是生态真变了还是分级变了。GEE里的分级可以用expression或者where实现var rseiClassified rsei.expression( (rsei 0 rsei 0.2) ? 1 : (rsei 0.4) ? 2 : (rsei 0.6) ? 3 : (rsei 0.8) ? 4 : 5, {rsei: rsei} );5.2 多时期变化检测RSEI最有价值的应用场景其实是时间序列变化分析。比如评估一个城市十年的生态环境演变或者一个生态修复工程实施前后的效果对比。操作思路很简单分别计算2005年、2010年、2015年、2020年等时间节点的RSEI然后做差值或矩阵转移分析。差值法就是后一期RSEI减前一期RSEI正值代表生态改善负值代表生态退化。矩阵转移分析则需要统计不同等级之间转移的面积。GEE里做差值非常方便var rsei2020 ...; // 2020年RSEI var rsei2010 ...; // 2010年RSEI var change rsei2020.subtract(rsei2010);然后设置阈值判断改善区域和退化区域。但我建议差值分析用连续RSEI值做分类转移矩阵则用分级后的结果做。两条技术路线互补一个评估程度一个评估类别的转变轨迹。5.3 制图输出的规范RSEI制图有几个容易忽略但影响观感的地方第一是配色方案。生态质量从差到优我推荐使用从红色经黄色、绿色到深绿的连续色带。具体来说差用暗红#8B0000较差用橙红#FF4500中等用黄色#FFD700良用浅绿#32CD32优用深绿#006400。这套配色有很直观的生态认知感审稿人也习惯这种表达。千万别用蓝色当优因为蓝色通常被读者理解为水体。第二是统计信息。图上一定要附上各等级面积和百分比统计表。很多人在ArcGIS里割了图就走评委一问中等生态区占比多少就答不上来这就尴尬了。第三是空间位置信息。研究区示意图、行政区边界、比例尺、指北针一个都不能少这是学术制图的基本盘。GEE里用Export.image.toDrive导出整幅栅格再拿到ArcGIS或QGIS里做正式制图。导出时注意scale设置成30米对应Landsat分辨率crs设置成UTM投影maxPixels适当调大。6. 我踩过的坑和给你的建议RSEI的坑说多不多但每一个都能让结果偏得离谱。我把这几年实际处理中踩过、以及看别人踩过的坑集中说一下。6.1 坑一归一化时被0除这是最常见也最隐蔽的问题。归一化的分母是max - min如果某个波段的max和min相等比如研究区内全是均匀地物或者掩膜后有效像元极少分母为0结果全面变成NoData或者Infinity。解决办法有两个一是归一化之前先做掩膜统计确保有效像元数量充足二是在公式里加一个极小值偏移量防止除零// 加一个极小值防止除零 var range max.subtract(min).add(0.0001); var normalized img.subtract(ee.Image.constant(min)).divide(ee.Image.constant(range));6.2 坑二PCA主成分被单一指标主导有次我帮人复查数据发现他的PC1特征向量里LST载荷达到0.9以上其他三个分量几乎可以忽略。查了一圈原因一是他直接用原始LST没有归一化把整个分析带偏了二是他的研究区有大量火烧迹地LST异常高而且没做掩膜。后来我把归一化修好、火烧迹地掩掉PC1的贡献结构就正常了。所以每次算出特征向量后我都建议顺手把载荷打印出来看一眼四个分量载荷量级应该相差不多如果某一个分量载荷超过0.8基本可以断定预处理出了问题。6.3 坑三大水体干扰严重LST在大水面上会异常偏低NDVI在水体上是负值WET在水体上反而极高。如果不处理整个PCA会被水体带节奏尤其是南方湖泊密集的研究区这个问题特别突出。常规做法是先用水体指数如MNDWI把水体掩膜掉再算RSEI最后制图时把水体单独叠加回来。这样既避免了水体对统计的干扰图面上也能看到水体的空间分布。公式是var mndwi image.normalizedDifference([Green, SWIR1]).rename(MNDWI); var waterMask mndwi.gt(0.2).not(); // MNDWI 0.2 视为水体掩膜掉6.4 经验和流程建议最后总结一套我在多个项目里验证过的稳定流程你可以直接照抄明确研究区和时间范围确定影像年份和季节夏季晴天优先。在GEE里批量筛选目标年份、云量低于10%的Landsat 8/9影像做云掩膜后取中值合成。计算NDVI、WET、LST、NDBSI四个分量并做水体掩膜。每个分量分别做2%-98%分位数截断归一化。合成后做PCA检查PC1贡献率需大于80%和特征向量方向。对PC1归一化得RSEI按需做方向校正。五级分类、统计各等级面积占比、导出栅格。在ArcGIS/QGIS里出图配上等级统计表和比例尺。做RSEI这个模型公式本身不复杂难点全在预处理和参数细节上。如果你能在自己的研究区把上面这八步完整跑通一遍你会发现它不仅是一个评价工具更是一套理解自然-人类耦合系统的思维方式。后续如果你想扩展还可以把RSEI和人口密度、GDP等社会经济数据结合做相关性分析或者用它来筛选生态修复的优先区域这都是一旦用熟了RSEI之后很自然的延伸方向。 遥感生态指数RSEI模型原理与GEE实现全解析
这些年做生态环境遥感评估我周围很多同行一开始都习惯用NDVI单指标来看一个区域的生态状况但说实话单看绿度太片面了。一个城市周边可能有大量农田NDVI很高但地表温度、建筑裸土占比这些生态压力因素完全没体现出来。后来接触到徐涵秋老师提出的遥感生态指数RSEI模型才真正觉得找到了一个既科学严谨、又能快速出结果的多指标评价框架。这篇文章我就把这套方法从原理到实操完整梳理一遍给正在做生态环境评价、城市规划研究、或者毕业论文需要出生态质量图的朋友做个参考。RSEI的核心思路其实很直接把影响生态环境的四大要素——绿度、湿度、热度、干度——用遥感数据分别量化再用主成分分析把它们客观地压缩成一个综合指数。整个流程在ENVI、ArcGIS或者Google Earth Engine上都能跑通关键是理解每一步为什么要这么做。1. RSEI模型到底在评价什么四大分量的生态语义在动手算之前我们得先把RSEI的四个输入分量吃透。这决定了你后面选数据、调参数的时候心里有没有底。1.1 绿度植被覆盖的健康表绿度指标用的是NDVI即归一化植被指数。它的原理是利用植被在红光波段强吸收、近红外波段强反射的光谱特性把植被覆盖信息从地表信号里分离出来。公式是NDVI (NIR - Red) / (NIR Red)对Landsat 8/9来说就是Band5 - Band4/Band5 Band4。NDVI的取值范围在-1到1之间理论上裸土接近0茂密森林可以到0.8以上水体为负值。我在实际处理中最想提醒的一点是NDVI虽然看起来简单但不同季节的植被状态差异极大。如果是做多时期对比影像的获取月份必须尽量一致比如都用每年8-9月份的影像否则你很难区分生态变化和物候变化。1.2 湿度土壤含水与地表湿润度湿度分量来自缨帽变换Tasseled Cap Transformation的湿度分量WET。缨帽变换是一种针对植被、土壤、水体等典型地物设计的线性变换它能把原始多波段数据压缩成几个有明确生态含义的组分亮度、绿度和湿度。不同传感器的WET系数差别很大这点非常容易踩坑。Landsat 8 OLI地表反射率产品的WET系数是WET 0.1511 * Blue 0.1973 * Green 0.3283 * Red 0.3407 * NIR - 0.7117 * SWIR1 - 0.4559 * SWIR2注意这里的系数是基于地表反射率数据推导的如果你用的是大气表观反射率TOA产品系数必须换一套。另外有些教程里给的是Landsat 5 TM的系数直接套到Landsat 8上结果出来全是乱的这个坑后面我会详细讲。1.3 热度地表温度的生态压力热度分量就是地表温度LST。城市热岛效应是生态质量恶化的最直接表现之一所以在RSEI里热度被当作重要压力指标。LST的反演方法有很多种比如辐射传输方程法、劈窗算法、单窗算法等。实操中我们最常用的是辐射传输方程法大气校正法它的基本思路是先把热红外波段的DN值转换为辐射亮度再转换为亮温然后结合地表比辐射率做发射率校正最终得到真实的地表温度。这个指标的计算链路最长涉及到好几个中间参数很多新手就是在这步被劝退的。别急第3章我会把每一步的推导和代码都给你列出来。1.4 干度建筑与裸土的人工印记干度指标NDBSINormalized Difference Bare Soil and Building Index是裸土指数SI和建筑指数IBI的均值NDBSI (SI IBI) / 2为什么要用两个指数的均值因为单一指数很难同时兼顾自然裸土和人工建筑。比如传统的裸土指数容易把建筑也识别进来建筑指数又容易漏掉裸土。把两者平均能在一定程度上相互弥补将这两类不透水面和裸露地表统一表征为生态干度。这个指标是RSEI里最能反映人类活动干扰强度的分量。一个区域如果NDBSI高通常意味着城市化程度高、地表硬化严重或者水土流失明显这些都属于生态退化的重要信号。理解了这四大分量的生态含义之后你会发现RSEI其实构建了一个压力-状态-响应的简化框架植被状态绿度、湿度代表生态系统的基础状况人类活动干扰热度、干度代表生态压力。综合起来就得到了一个能同时反映自然禀赋和人为影响的生态质量指数。2. 数据源选型与预处理决定RSEI成败的第一道关口方向不对努力白费。数据选型和预处理的质量直接决定了后面所有计算结果的可靠性。2.1 Landsat 8/9还是Sentinel-2怎么选目前做RSEI用得最多的是Landsat系列主要原因是它拥有从1980年代至今的连续存档非常适合做几十年的生态变化分析。具体到传感器我推荐以下选型原则数据源分辨率优势注意点Landsat 8/9 OLITIRS30m热红外100m2013年至今存档连续有成熟的热红外波段重访周期16天受云影响大Landsat 5 TM30m热红外120m1984-2011年历史数据数据较老辐射定标参数需仔细处理Landsat 7 ETM30m1999年至今2003年后有条带需要插值修复Sentinel-2 MSI10/20m分辨率高重访周期5天没有热红外波段LST需要单独处理有一点必须澄清Sentinel-2本身没有热红外传感器直接用Sentinel-2做完整的RSEI是算不出LST的。虽然有些研究用Landsat的LST降尺度配合Sentinel-2做高分辨率RSEI但那超出了本文基础教程范畴。我建议新手先老老实实用Landsat 8/9的Collection 2 Level 2产品。2.2 云掩膜与影像合成云是光学遥感最大的敌人。一片薄云或者云阴影覆盖的区域反射率信号完全失真如果没做掩膜就参与计算那几个像元的NDVI、LST都会变成离谱的异常值直接影响后续PCA的特征向量。在Google Earth Engine里我们可以直接用Landsat Collection 2 Level 2产品自带的QA_PIXEL波段做云掩膜这是目前最省事也最靠谱的做法。我常用的JavaScript代码如下function maskL8sr(image) { // 获取QA波段 var qa image.select(QA_PIXEL); // 设置要掩膜的位标志云、云阴影、冰雪 var cloudBitMask 1 1; // Cloud var shadowBitMask 1 3; // Cloud Shadow var snowBitMask 1 4; // Snow // 创建掩膜 var mask qa.bitwiseAnd(cloudBitMask).eq(0) .and(qa.bitwiseAnd(shadowBitMask).eq(0)) .and(qa.bitwiseAnd(snowBitMask).eq(0)); // 去掉缩放系数保留SR波段和热红外波段 return image.updateMask(mask) .select([SR_B2,SR_B3,SR_B4,SR_B5,SR_B6,SR_B7,ST_B10], [Blue,Green,Red,NIR,SWIR1,SWIR2,Thermal]); }这里要特别注意热红外波段ST_B10的定标单位。Collection 2 Level 2的地表温度产品单位是开尔文乘以10在后续计算时需要先乘以0.00341802并加上149.0来还原成真实的开尔文温度这个细节很坑后面演示代码里我会带上。如果你用的是ENVI/ArcGIS路线那么下载数据时应优先选Landsat Collection 2 Level 2科学产品而不是Level 1原始DN产品。Level 2已经帮你完成了大气校正可以直接拿SR波段做指数计算省掉了最耗时的大气校正环节。2.3 投影、裁剪与统一分辨率RSEI的四个分量中NDVI、WET、NDBSI是反射率比值或线性组合不受投影影响但LST是热红外波段计算来的原始分辨率是100米。如果在ENVI里直接做波段合成需要先把热红外重采样到30米否则最后PCA的输入图层分辨率不一致会产生错位。在GEE中处理非常简单用select和resample即可// 重采样到30米并统一投影 var thermal image.select(Thermal).resample(bilinear); var sortedImage image.addBands(thermal, null, true);统一投影时建议全部转到UTM投影。如果研究区跨多个UTM带可以用Albers等积投影保证面积计算不变形。另外别忘了裁剪到研究区边界。用clip接研究区的FeatureCollection即可var roi ee.FeatureCollection(projects/your-project/assets/roi); var imageClipped sortedImage.clip(roi);这一步能显著减少后续PCA的计算量尤其是研究区只占影像一小部分的时候。3. 四大指标计算实操公式、参数与GEE代码预处理做完就到了最核心的指标计算环节。我直接以Landsat 8/9 Collection 2 Level 2为输入给出完整的GEE实现。3.1 NDVI与WET的计算NDVI前面已经给了公式在GEE里用normalizedDifference一行就能解决var ndvi image.normalizedDifference([NIR, Red]).rename(NDVI);WET的计算需要用到缨帽变换系数。这里我给出Landsat 8 OLI地表反射率对应的完整公式var wet image.expression( 0.1511 * B2 0.1973 * B3 0.3283 * B4 0.3407 * B5 - 0.7117 * B6 - 0.4559 * B7, { B2: image.select(Blue), B3: image.select(Green), B4: image.select(Red), B5: image.select(NIR), B6: image.select(SWIR1), B7: image.select(SWIR2) } ).rename(WET);我把这一段单独列出来是因为系数搞混的概率实在太高了。网上流传的缨帽变换系数有很多版本有基于TOA反射率的、有基于地表反射率的、有TM的、有ETM的用错一个版本湿度分量就可能出现大面积负数后面PCA的结果也会跟着错。3.2 Landsat 8地表温度LST反演的完整链路LST是RSEI四个分量里计算链路最长的我按顺序把每一步拆开讲清楚。第一步读取热红外波段的DN值并定标。GEE里ST_B10已经是Level 2处理过的地表温度产品但单位是开尔文×10需要还原var thermalDN image.select(Thermal); var lstKelvin thermalDN.multiply(0.00341802).add(149.0).rename(LST_K);如果用的是Collection 2 Level 1的原始数据则需要先做辐射定标得到辐射亮度L再用热红外波段的K1、K2常数换算亮温。Landsat 8 Band 10的K1774.8853、K21321.0789这是当前Collection 2版本的官方常数不同版本有差异务必查询影像MTL文件// 如果从L1数据出发DN转辐射亮度 var ML 0.0003342; // 从MTL文件获取 var AL 0.1; // 从MTL文件获取 var radiance image.select(B10).multiply(ML).add(AL); // 辐射亮度转亮温 var brightnessTemp radiance.expression( K2 / log(K1 / rad 1), {rad: radiance, K1: 774.8853, K2: 1321.0789} );第二步计算地表比辐射率。这一步需要用NDVI估算植被覆盖度FVC再根据植被和裸土的比辐射率加权平均// 植被覆盖度 var ndviMin 0.05; var ndviMax 0.75; var fvc ndvi.subtract(ndviMin).divide(ndviMax - ndviMin).clamp(0, 1).rename(FVC); // 地表比辐射率水体取0.99城镇和裸地取0.97自然地表用FVC加权 var emissivity fvc.expression( (FVC 0) (NDVI 0.05) ? 0.995 : ((NDVI 0.05 NDVI 0.75) ? (0.004 * FVC 0.986) : 0.99), {FVC: fvc, NDVI: ndvi} ).rename(EM);第三步用单窗算法反演真实地表温度。经典的简化公式是LST BT / [1 (λ * BT / ρ) * ln(ε)]其中λ是热红外波段中心波长Landsat 8 Band 10取10.895微米ρ h * c / σ ≈ 1.438 * 10^-2 m·Kε为比辐射率。代码实现如下var lambda 10.895; // 微米 var rho 1.438e-2; // 米·开尔文 var lst brightnessTemp .divide(brightnessTemp.multiply(lambda).divide(rho).log().multiply(-1).add(1))等等这个公式我写反了。正确的单窗算法简化式是LST BT / (1 (λ * BT / ρ) * ln(ε))因为对数项是负值ε 1所以分母小于1LST会比亮温略高。GEE里正确写法是var lst brightnessTemp.divide( ee.Image(1).add( brightnessTemp.multiply(lambda).divide(rho).multiply(emissivity.log()) ) ).rename(LST);最后别忘了把开尔文转摄氏度方便制图时标注var lstC lst.subtract(273.15).rename(LST_C);3.3 NDBSI建筑与裸土指数计算NDBSI包含SI和IBI两个指数先分别算再取平均。SI的公式为SI [(SWIR1 Red) - (NIR Blue)] / [(SWIR1 Red) (NIR Blue)]IBI的公式为IBI {2SWIR1/(SWIR1NIR) - [NIR/(NIRRed) Green/(GreenSWIR1)]} / {2SWIR1/(SWIR1NIR) [NIR/(NIRRed) Green/(GreenSWIR1)]}从结构上可以看出IBI利用的是建筑在短波红外波段的反射率明显高于植被的特点以及植被在近红外高反射、建筑在绿光波段反射率相对高的光谱差异。代码实现如下var si image.expression( ((SWIR1 Red) - (NIR Blue)) / ((SWIR1 Red) (NIR Blue)), { SWIR1: image.select(SWIR1), Red: image.select(Red), NIR: image.select(NIR), Blue: image.select(Blue) } ).rename(SI); var ibi image.expression( (2*SWIR1/(SWIR1NIR) - (NIR/(NIRRed) Green/(GreenSWIR1))) / (2*SWIR1/(SWIR1NIR) (NIR/(NIRRed) Green/(GreenSWIR1))), { SWIR1: image.select(SWIR1), NIR: image.select(NIR), Red: image.select(Red), Green: image.select(Green) } ).rename(IBI); var ndbsi si.add(ibi).multiply(0.5).rename(NDBSI);到这里NDVI、WET、LST、NDBSI四个分量影像就全部准备好了。接下来的PCA才是RSEI的灵魂。4. 主成分分析PCA与RSEI合成从四个指标到一个指数为什么要用PCA而不是简单地给四个指标各赋一个权重求平均这是RSEI模型最有技术含量的地方。4.1 为什么先做归一化四个分量的量纲和数值范围差异巨大NDVI大致在-1到1之间LST可能在280到320开尔文WET的范围因传感器而异NDBSI也在-1到1附近。如果不归一化直接做PCA数值范围大的LST会主导主成分的方差贡献导致特征向量严重偏向热度分量其余分量的信息被压制。归一化公式是每个指标按0-1区间拉伸NI (I - I_min) / (I_max - I_min)在GEE里实现需要注意这里的I_min和I_max应该基于研究区影像的实际像元分布取而不是理论值。但直接取影像的最小最大值容易被异常像元干扰。我习惯用2%到98%的分位数来截断归一化这样能有效抑制极少数云残留或水体异常值的影响function normalize(img) { var min img.reduceRegion({ reducer: ee.Reducer.percentile([2]), geometry: roi, scale: 30, maxPixels: 1e10 }).values().get(0); var max img.reduceRegion({ reducer: ee.Reducer.percentile([98]), geometry: roi, scale: 30, maxPixels: 1e10 }).values().get(0); return img.subtract(ee.Image.constant(min)).divide(ee.Image.constant(max).subtract(ee.Image.constant(min))) .clamp(0, 1); }这里get(0)取到的值是数组形式在GEE里要小心处理。更稳妥的做法是用reduceRegion返回的字典按波段名分别取min和max。另外注意maxPixels要设得足够大否则大范围研究区会报错。4.2 PCA的输入与主成分贡献率判断归一化完成后将四个分量合成一个多波段影像然后在GEE里用reduceRegion配合ee.Reducer.principalComponentAnalysis()做PCA。PCA的核心逻辑是求四个波段构成的协方差矩阵的特征值和特征向量按特征值大小排序提取第一主成分PC1。var image4bands ndvi.addBands([wet, lst, ndbsi]).select([NDVI,WET,LST,NDBSI]).float(); var pca ee.Reducer.principalComponentAnalysis(4); var pcImage image4bands.reduceNeighborhood(pca, ee.Kernel.square(1));等等reduceNeighborhood不是我们想要的方案。它是在滑动窗口里做PCA计算量和结果的解释方式都不对。正确的做法应该是先在全影像范围内采样构建协方差矩阵再对影像做线性变换// 先采样得到均值与协方差 var mean image4bands.reduceRegion({ reducer: ee.Reducer.mean(), geometry: roi, scale: 30, maxPixels: 1e10 }); var centered image4bands.select([NDVI,WET,LST,NDBSI]) .subtract(ee.Image.constant(mean.values())); var covar centered.toArray().reduceRegion({ reducer: ee.Reducer.covariance(), geometry: roi, scale: 30, maxPixels: 1e10 }); // 从协方差做PCA var eigen ee.Array(covar.get(array)).eigen(); var eigenVectors eigen.slice(1, 0, 1).slice(0, 0, 4); var pc1 centered.toArray().matrixMultiply(eigenVectors.transpose()).arrayProject([0]) .arrayFlatten([[PC1]]);这一段代码有几个技术细节值得展开eigen()返回的是一个二维数组每一行格式是[特征值, 特征向量分量...]按特征值降序排列。所以slice(1, 0, 1)取出的是第一行后面的部分即最大特征值对应的特征向量。这个特征向量就是PC1对应的权重它的四个分量分别代表NDVI、WET、LST、NDBSI在PC1中的载荷方向。拿到了PC1之后我们要看两个关键指标第一是特征值贡献率。PC1对应的特征值占所有特征值之和的比例一般要求大于80%。如果贡献率偏低说明四个分量的信息没有很好地收敛到一个主轴上模型解释力不够。第二是特征向量的符号。在原始的RSEI框架中我们期望NDVI和WET的载荷为正生态好的指标贡献正向LST和NDBSI的载荷为负生态压力指标贡献负向。如果你算出来的PC1载荷方向完全相反即NDVI为负、LST为正说明第一主成分表达的是生态退化而不是生态质量。这种情况下需要把RSEI反过来用1减去PC1来校正方向否则高值代表恶劣、低值代表良好语义全反了。4.3 RSEI合成、再归一化与方向校正无论特征向量方向如何最终都要做一个关键步骤把PC1再归一化到0-1区间然后判断是否需要方向修正。// 先归一化PC1到0-1 var pc1Norm normalize(pc1); // 判断方向并校正 // 如果NDVI、WET在特征向量中的系数之和为正说明PC1越大生态越好直接用pc1Norm // 如果为负则需要用 1 - pc1Norm var rsei pc1Norm; // 或者 rsei ee.Image(1).subtract(pc1Norm);这里我强烈建议你打印出特征向量来人工检查一下别全指望自动判断。具体做法是print(eigenVectors)在Console里查看四分量各自的载荷值。我处理过不少区域大部分时候NDVI和WET的载荷为正但也有例外——比如在大范围荒漠地区PC1可能完全由LST主导这时候就要特别小心。合成之后的RSEI理论上取值范围是0到1数值越高代表生态质量越好。但实际像元分布往往不均匀可能集中在0.4到0.7之间。这时候我会做一次直方图拉伸或分位数拉伸方便后续分级和可视化。5. 结果分级、变化检测与制图输出RSEI算出来是个连续变量但研究上最常用的呈现方式是分级。分级既方便空间对比也好做面积统计。5.1 五级分类阈值的设定目前学界用得最多的是等间距五级分类法把0-1区间均分为五等份RSEI区间等级生态含义0.0 - 0.2差生态质量极差以高强度开发或严重退化区域为主0.2 - 0.4较差生态质量偏差人为干扰明显0.4 - 0.6中等生态质量一般自然与人为活动交织0.6 - 0.8良生态质量较好植被覆盖较高0.8 - 1.0优生态质量优异基本为茂密植被区等间距分类最大的优势是可对比性。不同年份、不同研究区的RSEI结果都能用同一套标准比较。不过也有学者用自然断点法Jenks分级让每个类别的内部差异最小化。我自己做多时期对比时坚持用等间距法因为自然断点法每年的断点不同变化统计的结果会被分类方法本身干扰很难说清是生态真变了还是分级变了。GEE里的分级可以用expression或者where实现var rseiClassified rsei.expression( (rsei 0 rsei 0.2) ? 1 : (rsei 0.4) ? 2 : (rsei 0.6) ? 3 : (rsei 0.8) ? 4 : 5, {rsei: rsei} );5.2 多时期变化检测RSEI最有价值的应用场景其实是时间序列变化分析。比如评估一个城市十年的生态环境演变或者一个生态修复工程实施前后的效果对比。操作思路很简单分别计算2005年、2010年、2015年、2020年等时间节点的RSEI然后做差值或矩阵转移分析。差值法就是后一期RSEI减前一期RSEI正值代表生态改善负值代表生态退化。矩阵转移分析则需要统计不同等级之间转移的面积。GEE里做差值非常方便var rsei2020 ...; // 2020年RSEI var rsei2010 ...; // 2010年RSEI var change rsei2020.subtract(rsei2010);然后设置阈值判断改善区域和退化区域。但我建议差值分析用连续RSEI值做分类转移矩阵则用分级后的结果做。两条技术路线互补一个评估程度一个评估类别的转变轨迹。5.3 制图输出的规范RSEI制图有几个容易忽略但影响观感的地方第一是配色方案。生态质量从差到优我推荐使用从红色经黄色、绿色到深绿的连续色带。具体来说差用暗红#8B0000较差用橙红#FF4500中等用黄色#FFD700良用浅绿#32CD32优用深绿#006400。这套配色有很直观的生态认知感审稿人也习惯这种表达。千万别用蓝色当优因为蓝色通常被读者理解为水体。第二是统计信息。图上一定要附上各等级面积和百分比统计表。很多人在ArcGIS里割了图就走评委一问中等生态区占比多少就答不上来这就尴尬了。第三是空间位置信息。研究区示意图、行政区边界、比例尺、指北针一个都不能少这是学术制图的基本盘。GEE里用Export.image.toDrive导出整幅栅格再拿到ArcGIS或QGIS里做正式制图。导出时注意scale设置成30米对应Landsat分辨率crs设置成UTM投影maxPixels适当调大。6. 我踩过的坑和给你的建议RSEI的坑说多不多但每一个都能让结果偏得离谱。我把这几年实际处理中踩过、以及看别人踩过的坑集中说一下。6.1 坑一归一化时被0除这是最常见也最隐蔽的问题。归一化的分母是max - min如果某个波段的max和min相等比如研究区内全是均匀地物或者掩膜后有效像元极少分母为0结果全面变成NoData或者Infinity。解决办法有两个一是归一化之前先做掩膜统计确保有效像元数量充足二是在公式里加一个极小值偏移量防止除零// 加一个极小值防止除零 var range max.subtract(min).add(0.0001); var normalized img.subtract(ee.Image.constant(min)).divide(ee.Image.constant(range));6.2 坑二PCA主成分被单一指标主导有次我帮人复查数据发现他的PC1特征向量里LST载荷达到0.9以上其他三个分量几乎可以忽略。查了一圈原因一是他直接用原始LST没有归一化把整个分析带偏了二是他的研究区有大量火烧迹地LST异常高而且没做掩膜。后来我把归一化修好、火烧迹地掩掉PC1的贡献结构就正常了。所以每次算出特征向量后我都建议顺手把载荷打印出来看一眼四个分量载荷量级应该相差不多如果某一个分量载荷超过0.8基本可以断定预处理出了问题。6.3 坑三大水体干扰严重LST在大水面上会异常偏低NDVI在水体上是负值WET在水体上反而极高。如果不处理整个PCA会被水体带节奏尤其是南方湖泊密集的研究区这个问题特别突出。常规做法是先用水体指数如MNDWI把水体掩膜掉再算RSEI最后制图时把水体单独叠加回来。这样既避免了水体对统计的干扰图面上也能看到水体的空间分布。公式是var mndwi image.normalizedDifference([Green, SWIR1]).rename(MNDWI); var waterMask mndwi.gt(0.2).not(); // MNDWI 0.2 视为水体掩膜掉6.4 经验和流程建议最后总结一套我在多个项目里验证过的稳定流程你可以直接照抄明确研究区和时间范围确定影像年份和季节夏季晴天优先。在GEE里批量筛选目标年份、云量低于10%的Landsat 8/9影像做云掩膜后取中值合成。计算NDVI、WET、LST、NDBSI四个分量并做水体掩膜。每个分量分别做2%-98%分位数截断归一化。合成后做PCA检查PC1贡献率需大于80%和特征向量方向。对PC1归一化得RSEI按需做方向校正。五级分类、统计各等级面积占比、导出栅格。在ArcGIS/QGIS里出图配上等级统计表和比例尺。做RSEI这个模型公式本身不复杂难点全在预处理和参数细节上。如果你能在自己的研究区把上面这八步完整跑通一遍你会发现它不仅是一个评价工具更是一套理解自然-人类耦合系统的思维方式。后续如果你想扩展还可以把RSEI和人口密度、GDP等社会经济数据结合做相关性分析或者用它来筛选生态修复的优先区域这都是一旦用熟了RSEI之后很自然的延伸方向。