1. 项目概述与核心需求解析1.1 这个项目解决了什么问题做遥感的人绝大多数时间都在跟NDVI、地表温度、水体提取这些常规指标打交道。但真正落到山区、峡谷、高纬度地区时山体阴影这个看似简单的问题往往会变成一道很难绕过去的坎。光伏电站选址需要避开阴影区域生态学里评估林下光照条件需要知道阴影分布农田遥感反演时阴影区的地表反射率会明显偏低如果没做阴影识别就直接拿去算植被指数结果必然被系统性拉低。你可能会想MODIS那么粗的分辨率动不动就是250米到1公里的像元算山体阴影有什么意义恰恰相反MODIS的应用场景正是大面积、长时序的宏观分析。比如整个青藏高原的太阳能资源评估、横断山区的植被生产力修正、甚至全国尺度的地形遮蔽效应统计——这些工作在Landsat 30米尺度下数据量是天文数字但在MODIS尺度下几个Tile一拼就能快速出结果。这个项目要做的就是在Google Earth Engine题目标注为Open Earth Engine library实际指的是GEE这个开放的地球观测云计算平台里用MODIS数据配合DEM高程模型计算出研究区域的山体阴影分布并精确统计出被阴影覆盖的区域面积。先说结论整个方案不需要本地下载任何影像所有计算全部在云端完成跑一个中等省份范围的分幅统计一般几十秒就能出结果。这个效率你让ENVI或ArcGIS处理同样的数据量光导出就要等到怀疑人生。1.2 适合谁使用这套方案我把话放前面这个方案对以下三类人群价值最大做新能源选址的技术人员。光伏电站、风电场的日照评估需要对阴影分布做定量分析但很多人手里只有CAD和气象站数据没有处理遥感影像的完整流程。GEE这套方案可以直接一键出面积统计表。高校和科研院所的生态遥感团队。做高寒草甸、林线过渡带、冰川变化研究的都需要区分地形阴影与真实地表变化这个流程可以作为数据预处理的标准模块。刚接触GEE的初学者。因为这套方案涉及影像筛选、波段计算、区域统计、面积换算、结果导出这些最基本的操作练完这一整套你对GEE的理解会上一个台阶。我在下面几节里会把数据选型逻辑、山体阴影计算的原理、代码实现、以及我踩过的坑全部展开讲清楚。2. MODIS数据选型与GEE环境准备2.1 为什么选择MODIS而不是Landsat这是很多新手最先问的问题。Landsat的空间分辨率更高做小范围阴影识别确实更精细但一旦研究区扩大到省级以上Landsat的条带拼接和时间合成会带来非常麻烦的问题。MODIS的优势有三个幅宽达到2330公里一天就能覆盖全球绝大部分区域研究区跨多个省也不需要手动拼接。有成熟的大气校正产品和逐像元质量评估波段State QA处理阴影问题时可用的波段信息更完整。GEE里MODIS Collection 6.1数据已经做了标准化处理反射率产品自带缩放系数直接使用非常省心。当然MODIS的短板也很明显几何分辨率粗250m/500m/1000m你无法指望它识别出单棵树的阴影。但它做的是区域尺度的宏观统计这个尺度下地形阴影的主导因素是山体本身的坡度坡向而不是局部地物。2.2 MODIS产品选择MOD09GA还是MOD09A1做山体阴影分析我们需要的是地表反射率产品。GEE里常见的MODIS产品有MOD09GA每日表面反射率和MOD09A18天合成表面反射率。如果你要研究的是特定时刻或特定过境日期的阴影分布用MOD09GA如果你要的是某个月份或季度内具有代表性的阴影状态用MOD09A1。从数据质量稳定性来说MOD09A1通过8天合成能有效去除云污染和异常值对于面积统计类的项目我更推荐它。但需要注意的是MODIS本身携带的角度信息里可以直接获取太阳高度角和方位角而阴影计算需要这些参数。MOD09GA的SR_Date和时间属性里就包含了太阳位置信息这个对我们的项目非常关键。2.3 GEE环境初始化与数据导入方法如果你还没用过GEE先花两分钟注册一个账号并访问GEE代码编辑器Code Editor。所有处理都在浏览器内完成不需要安装任何本地包。如果你习惯Python也可以选择用Python API但下面我统一用JavaScript版本的GEE代码编辑器来讲因为它在数据可视化方面更直接。初始化时建议先加载研究区边界、DEM数据和MODIS影像集。代码框架如下// 设置研究区以某个山区县为例这里用几何坐标近似表示 var roi ee.Geometry.Polygon([ [[101.5, 27.0], [102.5, 27.0], [102.5, 28.0], [101.5, 28.0]] ]); // 加载SRTM 30米DEM数据 var dem ee.Image(USGS/SRTMGL1_003).clip(roi); // 加载MODIS MOD09GA每日地表反射率产品 var modis ee.ImageCollection(MODIS/061/MOD09GA) .filterBounds(roi) .filterDate(2023-11-01, 2023-11-30) .first(); print(modis);这里有几个容易踩的坑第一ROI范围不要太大否则后续shadow计算时云计算耗时很长甚至超时第二MODIS数据一定要指定日期草地和裸地的反射率差异会干扰阈值的选取第三DEM与MODIS的坐标系是不同源数据自动融合的GEE会自动重投影但我们在面积统计时要格外注意投影设置后面有详细说明。3. 山体阴影计算的原理与参数选择3.1 山体阴影的本质是什么山体阴影从遥感物理本质来说就是太阳入射光线被地形起伏遮挡后形成的暗区。这个暗区不一定等于“夜晚区”而是太阳直接辐射照不到、只有天空散射辐射到达的区域。计算山体阴影的标准方法是用山体阴影函数Hillshade。核心输入是三个变量太阳高度角Solar Elevation太阳相对地平线的高度角。太阳方位角Solar Azimuth太阳光线从正北方向顺时针转过的角度。地形参数坡度Slope和坡向Aspect。计算公式如下hillshade cos(90°-太阳高度角) * cos(坡度) sin(90°-太阳高度角) * sin(坡度) * cos(太阳方位角 - 坡向)当hillshade值小于某个阈值比如0或负数时说明该像元接收不到直接的太阳辐射即为阴影区。这里我提醒一个关键点很多人直接把遥感影像按亮度阈值切割来提取阴影这在山体地区经常出错因为水体、暗色裸岩、云影都会干扰。正确做法是先算地形阴影再用它来指导遥感影像的阴影识别。3.2 太阳位置的获取方式MODIS数据在采集时卫星过境时间基本固定太阳位置也就相对固定。我们要从MOD09GA的影像属性中读取太阳高度角与方位角。具体来说MOD09GA的metadata里包含SOLAR_ZENITH_ANGLE和SOLAR_AZIMUTH_ANGLE这类字段GEE里可以直接通过get方法获取。如果你的研究区很大覆盖多个MODIS条带不同条带的过境时间略有差异太阳位置也会有差别。这时候建议对每个影像分别计算hillshade再取均值避免统一参数带来的系统误差。但是在绝大多数项目里研究区内太阳位置的差异对面积统计的结果影响极小——真正影响大的是DEM的精度。3.3 DEM的选择SRTM与ASTER的区别GEE里最容易获得的DEM数据是SRTM30米和ASTER GDEM30米。两者对比SRTM在平坦地区噪声小ASTER在部分山区有更丰富的地形细节但也伴随更多空洞和异常值。我建议优先用SRTM。原因很简单山体阴影计算对地形高程的细微误差非常敏感SRTM的数据整体平顺性好不像ASTER在高海拔山区容易出现条带与碎点。如果你有更高分辨率的本地DEM也可以通过ee.Image.load上传但SRTM对于绝大多数场景已经够用。3.4 阈值的选择到底怎么定阴影和非阴影之间的切割阈值没有万能值。靠单一hillshade的0界线来分割在实际山地中往往不够准确因为清晨和傍晚时太阳高度角很低散射光占主导即便在受光面地表反射率也很低。我常采用的做法是组合使用两个条件阴影像元 (hillshade 0.15) 且原始影像红波段反射率 某数值这个“某数值”要根据研究区地表覆盖类型来定。冬季草地一般在0.05-0.1之间裸岩和沙地在0.12-0.2之间。你可以先输出一个阴影二值图叠加在影像上目视检查再调整阈值。4. 实操代码MODIS影像预处理与山体阴影面积计算4.1 MODIS影像的缩放与质量过滤MODIS MOD09GA产品的波段值实际上是原始整型数据乘以0.0001换算成反射率。计算前必须先缩放否则你拿到的反射率是0-32767这样的原始DN值没法用。同时需要进行云掩膜把云和云影的像元剔除否则云影会被错误统计成山体阴影。这里使用MOD09GA自带的StateQA波段// 读取MODIS影像波段 var surfRef modis.select(sur_refl_b01).multiply(0.0001).rename(RED); var stateQA modis.select(StateQA); // 云掩膜函数 function cloudMask(image) { var qa image.select(StateQA); var cloudBit qa.bitwiseAnd(2).eq(0); // bit 1 云像元标记 return image.updateMask(cloudBit); } var modisMasked modis .map(cloudMask) .select(sur_refl_b01) .multiply(0.0001) .rename(RED) .first();注意这里演示的是逐日影像配合云掩膜的做法。如果使用MOD09A1合成产品云掩膜逻辑略有不同但有一样的QA波段可参考。4.2 基于SRTM的坡度、坡向与hillshade计算在GEE中ee.Terrain.slope和ee.Terrain.aspect可以直接从DEM提取坡度坡向。然后我们根据MODIS metadata里获取的太阳高度角和方位角调用ee.Terrain.hillshade生成山体阴影栅格。// 提取DEM的坡度和坡向 var slope ee.Terrain.slope(dem); var aspect ee.Terrain.aspect(dem); // 获取MODIS过境时的太阳高度角和方位角 // 若使用单景影像直接从属性中读取 var elevation ee.Number(modis.get(SOLAR_ZENITH_ANGLE)).multiply(-1).add(90); var azimuth ee.Number(modis.get(SOLAR_AZIMUTH_ANGLE));这里有一个非常重要的换算问题MODIS数据里直接提供的是太阳天顶角Solar Zenith Angle而hillshade函数需要的是太阳高度角Solar Elevation Angle两者关系为太阳高度角 90° - 太阳天顶角千万别直接拿天顶角去算算出来的阴影全反过来了。4.3 hillshade二值化与阴影面积统计有了太阳高度角和方位角就可以计算山体阴影并生成二值栅格。二值栅格中值为1代表阴影区值为0代表非阴影区。然后用ee.Image.pixelArea()计算每个像元的真实地表面积考虑了地球曲率和投影面积变形把阴影区的像元面积累加就得到了山体阴影的总面积。这一步特别关键——如果直接用像元尺寸相乘在高纬度地区误差会非常大。// 计算山体阴影 var hillshade ee.Terrain.hillshade(dem, azimuth, elevation); // 生成阴影二值图像 var shadow hillshade.lt(0.15).rename(shadow); // 获取研究区内的像元面积平方米 var areaImage ee.Image.pixelArea().updateMask(shadow).clip(roi); // 统计阴影面积 var shadowArea areaImage.reduceRegion({ reducer: ee.Reducer.sum(), geometry: roi, scale: 250, maxPixels: 1e13, bestEffort: true }); print(阴影面积平方米:, shadowArea.get(area)); // 换算为平方千米 print(阴影面积平方千米:, shadowArea.get(area).divide(1e6));上面用的scale是250米这是为了匹配MODIS的波段分辨率。有些教程会把scale设为30米然后告诉你投影到SRTM的分辨率上更精细这是错误的。MODIS本身是250米分辨率你硬用30米的网格去插值统计看起来精细了实际是把原始数据重复采样并不会增加任何真实信息。4.4 对整个研究区逐块统计面积如果你的ROI是省级甚至全国范围直接对整个区域reduceRegion可能计算内存不足或超时。应对办法是把ROI划分成渔网格Fishnet逐个格网计算后汇总。// 创建渔网格每个格网大小为0.5度 var grid ee.FeatureCollection( ee.Geometry.Polygon(roi.coordinates()).coveringGrid(0.5, null, false) ); // 对每个格网计算阴影面积 var shadowStats grid.map(function(feature) { var area areaImage.reduceRegion({ reducer: ee.Reducer.sum(), geometry: feature.geometry(), scale: 250, maxPixels: 1e13, bestEffort: true }); return feature.set(shadow_area, area.get(area)); }); print(shadowStats);这种方法比一次性统计更稳定也方便后续做空间分布地图。5. 常见问题与排查技巧实录5.1 为什么算出来的阴影面积为0或极大我遇到最多的问题是阴影面积直接为0。排查思路如下先检查太阳高度角是不是负数。MODIS某些影像metadata里的太阳天顶角在早晚轨道时大于90度换算出来的高度角是负的此时全区域都是理论阴影面积就会极大。再检查hillshade阈值。常见反映率在0-255的产品没有乘以0.0001会导致后续判断条件永远不满足。最简单的排查方法是把hillshade直接加到地图上可视化看它是否呈现明显的明暗起伏而不是整片纯黑或纯白。5.2 MODIS云云影仍然混入怎么办如果你的研究区在热带或雨季云影问题非常头疼。山体阴影是用DEM算出来的地形阴影云影是云层遮挡造成的两者不是一个来源。解决云影的正确姿势是优先使用MOD09A1 8天合成产品。它已经做了时序合成能有效去云。如果你坚持每日数据建议用MOD09GA的StateQA把所有有云像元做掩膜同时在最终面积统计时将云掩膜后的有效像元作为分母计算出“有效观测范围内阴影面积占比”这样结果才有可比性。5.3 不同投影和scale设置对面积的影响有人习惯把scale设为30米认为更准确。其实MODIS原始分辨率只有250米你用30米去计算只是做了一个双线性/最近邻重采样得到的边缘像元会变得平滑但本质上没有新增真实信息。我在实际项目中测过同一ROI在250米和30米不同scale下的阴影面积结果差异在2%-5%之间。这个差异不算大但如果你统计的是省级大范围误差就会以百平方公里计。为了可重复性建议固定250米别折腾。还有一个坑是pixelArea()本身返回的是SRTM像元实际地表面积包含了坡度引起的面积放大效应。这个功能非常实用因为山区斜坡上的真实地表面积一定大于水平投影面积如果你直接拿水平面积去估算阴影覆盖结果就偏小。5.4 太阳高度角较低时的误判处理冬季太阳高度角低时大片坡面都处在阴影中hillshade阈值设为0.15会切割出过多的阴影。此时建议分两步走第一步分别统计坡向朝北与朝南的部分再做对比。这可以通过重分类坡向来实现var northAspect aspect.gt(90).and(aspect.lt(270));第二步对阴影栅格加载显示肉眼对比典型地物如河流峡谷、山脊线看边界是否符合地形特征。如果阴影边界明显越过了山脊线说明阈值过低适当调高。这里没有标准答案每个区域的太阳角度和地形条件都不同项目的核心价值就是让你能灵活调整并快速评估结果而不是机械套公式。5.5 常见问题速查表问题现象可能原因排查与解决办法阴影面积为0太阳高度角为负值或缩放系数未应用检查metadata中的天顶角换算确认反射率已乘0.0001阴影面积异常偏大云影被误判为山体阴影使用MOD09A1合成产品或加强QA云掩膜两次运行结果不一致使用了不同日期的MODIS影像固定日期或对多日期结果取平均大面积计算超时直接对整个ROI做reduceRegion改用coveringGrid分格网统计边界处阴影锯齿严重分辨率不匹配将scale统一设为250米避免产生过度锯齿6. 结果可视化与扩展应用思路6.1 在地图上直接叠加显示GEE最让人舒服的地方是能直接在地图上预览结果。把阴影二值影像叠加在MODIS假彩色影像上可以看出阴影区域是否与山系走向一致。Map.centerObject(roi, 10); Map.addLayer(modisMasked, {min: 0, max: 0.3}, MODIS RED); Map.addLayer(shadow.updateMask(shadow), {palette: [blue]}, Shadow);蓝色区域就是统计出的山体阴影。我每次都会这样先目视检查一遍再去看面积数字。6.2 多日期对比与阴影变化的时序分析单一时期的阴影面积只能给出一个静态结果但如果你做光伏选址需要知道从冬至到夏至阴影面积的季节变化。利用GEE的时序能力可以循环逐月计算阴影面积输出一条季节变化曲线。这个功能在环境评估里尤其有用判断山谷里的某块规划用地到底是全年光照充足还是只有夏季几个月能照到太阳。实操上只需把上面的计算封装成一个函数传入不同月份的MODIS影像集用map循环计算即可。不需要额外写复杂代码但对结果的解读要有清醒认识MODIS影像每天过境时刻是固定的这个时刻的阴影面积并不完全等于全天平均阴影面积。6.3 扩展到Landsat/Sentinel-2精细识别如果你觉得250米分辨率太粗也可以用同样的方法在Landsat或Sentinel-2上跑原理完全一致只是把影像集合和scale替换掉。但反过来我劝你意识到一个问题精细尺度下的阴影识别难点不再是面积统计而是云影、建筑物阴影、树木阴影的分离。那时候单靠hillshade地形阴影就不够了建议引入光谱形状分析或太阳-传感器入射余弦模型配合红蓝波段比值做多条件判定。7. 实操心得与几点忠告最后再分享一些我自己做这类项目时沉淀下来的经验。有人会纠结于MODIS的BRDF校正波段其实对于山体阴影识别来说BRDF调整不是必须的因为阴影识别的核心是地形和太阳角度而不是反射率绝对值。但另一个细节常常被忽略。MODIS的像元在地形起伏大的区域其实际地面覆盖范围与水平投影差异很大这种几何扭曲在陡峭峡谷里会造成边界位置偏移。如果你研究区的坡度超过30度直接套用hillshade阈值并不严谨需要结合太阳直射与散射比例进一步做辐射传输模型校正。学术级研究建议参考文献里的COS(z)模型但在工程预算有限的条件下用hillshade统计面积已经是性价比最高的方案。最后给一个技术建议。GEE的ee.Terrain.hillshade函数其实有两种调用参数一种是传入太阳方位角和高度角另一种是提取影像自带的太阳几何。两者结果略有差异。我的经验是用影像自带太阳几何因为山地环境下的云层和大气状况会影响有效照明角度用实际过境时的角度比用理论值更贴近地表真实光照状态。写到这里项目主体内容已经完整了。你拿到这套方案之后可以先用小范围研究区跑通流程再逐步放大到省域。每一步都能通过可视化图层验证面积统计结果也能导出成CSV或Shapefile用于后续制图。这套流程我重复用了不下二十次稳得很。