☰
S1与S2影像在GEE中的多波段合并:时间匹配与踩坑复盘
2026/10/15 11:29:27 网站建设 项目流程

两年前我第一次跑农田长势监测的时候,被一件事磨得没脾气——六月份正是看苗情的关键期,Sentinel-2光学影像却连续三旬被云层糊得干干净净;Sentinel-1雷达成像倒是风雨无阻,一周给我递一双新数据,可纯雷达图上的信号我又不敢轻易解读,因为后向散射对地表状态的响应方式和可见光反射率完全是两套语言。后来我在Google Earth Engine(GEE)里试着把Sentinel-1和Sentinel-2连接起来做分析,把同一时间段内的雷达极化波段和光学多光谱波段合并成一张复合影像,问题才真正解开。这篇文章就把我这套“同一时间多波段合并”的做法、代码和踩过的坑完整梳理一遍,适合正在做农业遥感、洪水监测、湿地分类或者任何需要光电协同的同行参考。

1. 雷达与光学影像的互补逻辑:为什么非要连接分析

1.1 光学看得见“内容”,雷达看得见“状态”

很多人第一次听到“连接分析”这个概念会觉得玄乎,其实说穿了特别简单:Sentinel-2记录的是地表对太阳光的反射,反映的是叶绿素含量、水分吸收、地表材质这些“物质属性”;Sentinel-1记录的是地表对C波段微波的后向散射,反映的是介电常数、表面粗糙度、目标几何结构这些“状态属性”。

我用一个生活化类比来解释:Sentinel-2像是你白天用眼睛看一棵树,能分辨出叶子绿不绿、密不密;Sentinel-1像是你在夜里用手电筒斜着扫同一棵树,看到的是枝干走向、树冠凹凸和叶子湿不湿。前者看“是什么”,后者看“怎么样”,两者本来就不是替代关系,而是正交的互补关系。

在实际遥感任务里,这种互补非常珍贵。比如洪水监测,Sentinel-2的短波红外波段能识别水体,但洪水期往往正好是云层最厚的时期,光学影像经常一景干净的都没有;而Sentinel-1的VH极化对静水表面非常敏感,水面会呈现强烈的前向镜面反射,后向散射信号断崖式下降,几乎不受云层干扰。反过来,纯雷达数据对作物的具体种类区分度不够,因为水稻、玉米、林地在某些物候期的后向散射范围会重叠,但加上Sentinel-2的红边波段和植被指数之后,分类精度会立刻上一个台阶。

1.2 什么场景必须做“光电连接分析”

不是所有任务都非要把S1和S2绑在一起,我自己判断的标准很简单:只要地表状态随时间变化比较快、同时光学影像又无法保证连续覆盖,那“连接分析”几乎就是刚需。

典型的场景有这么几个。第一个是作物分类和长势监测:光学波段提供NDVI、NDWI、红边参数,雷达的VV/VH提供冠层结构信息和土壤水分信号,两者组合以后,对水稻、玉米、小麦这类作物的识别精度能明显提升,尤其在生育期内不同作物的雷达时序差异会逐渐显现。第二个是洪涝灾害快速评估:洪水发生时往往连天阴雨,S2即便能成像也是大片云掩膜,这时候靠S1的VH下降幅度圈定淹没范围,再用淹没前的一帧S2做土地利用底图,直接估算不同地类下的受灾面积。第三个是土壤水分反演:光学数据能剥掉植被覆盖的影响,雷达数据则直接对表层土壤介电常数敏感,两者联合反演比单用其中一个稳健得多。

还有一个经常被忽略的用途是样本标注。做分类任务时,很多人用光学影像人工勾选样本点,但云遮挡区域根本没有光学底图可看;如果用S1的假彩色合成辅助判断地类边界,很多模糊的样本就能定下来。这种“雷达辅助光学解译”的做法,本质上也是一种连接分析。

2. 时间匹配窗口的选法:重访周期与地表变化速度的博弈

2.1 Sentinel-1与Sentinel-2的时间节奏差

连接分析最容易翻车的地方不在波段处理,而在时间匹配。“同一时间”这四个字说起来轻巧,做起来却有讲究。

Sentinel-1单颗卫星的重访周期是12天,但Sentinel-1A和1B双星编队飞行时能达到6天,在赤道附近结合升降轨甚至能拿到更密的观测;Sentinel-2A/B双星则是5天重访。看起来两边节奏都挺快,但别忘了Sentinel-2的5天是“理论重访”,真正有效的光学观测还要看云,所以实际应用中经常出现“S2已经一个月没给干净影像、S1倒是勤勤恳恳拍了七八景”的局面。

这就引出第一个原则:时间匹配窗口不能一刀切,必须根据地表变化速度和任务需求来定。变化越快的场景,比如洪水演进、台风过后、作物收割期,窗口必须收紧;变化相对缓慢的场景,比如作物营养状况监测、土壤水分旬报,窗口可以放宽到半月甚至一个月,用中值合成来换数据完整性。

2.2 严格同步与宽松同步:如何取舍

我一般把时间匹配策略分成两档。

“严格同步”指的是围绕某个关键时间点,S1和S2的获取日期相差不超过3到5天。这种做法适用于灾害应急,比如某天发生溃坝,我就以那一天的S1影像为锚点,向前找一景10天内的S2做地类底图,向后找一景暴雨后的S2做对比。代价是窗口内可能只有一两景影像,噪声压制不够,局部会有斑点残留。

“宽松同步”指的是按自然月或自定义的30天窗口,把该时段内的所有S1和S2影像分别做中值合成,再拼接。这种做法适用于周期性监测,比如物候分析、月度土地覆盖制图。中值合成能有效压制S1的斑点噪声和S2的云残影,影像看起来干净很多,但代价是“同一时间”变成了“同一个月”,地表如果在一个月内发生过剧变(比如中间来了一场强降雨),合成的S1和S2描述的实际是不同阶段的状态。

我自己的经验是:如果任务目标是分类制图,宽松同步足够;如果任务目标是变化检测或参数反演,至少要以S1的获取日期为准做逐景配对,不要简单拿两个月度合成硬拼。

2.3 以Sentinel-1日期为准的配对技巧

这里分享一个我后来一直在用的技巧:不要试图让S1去迁就S2的日期,而是反过来,以每一景S1的获取时间为中心,向前后各取一个时间窗(推荐10天),在这个窗口内找S2影像做中值或最邻近合成,然后得到一景“雷达时刻”的复合影像。

道理很简单:雷达不受云影响,它的时间序列是连续的,可以作为时间轴的主干;光学影像断断续续,应该被“挂载”到雷达时间节点上。比如一景S1是6月8日拍的,我就取5月29日到6月18日之间的所有合格S2影像,合成这一帧光学底图,与6月8日的S1捆绑成一景复合数据。这样生成的时间序列,每一个时间节点的雷达信号都是真实的、独立观测的,不会被月合成平均掉。

3. 合并前的数据预处理:坐标系、分辨率、去云和去噪

3.1 集合ID的选择:一个容易被忽略的坑

GEE里Sentinel数据有好几个集合,选错会让你后面所有的波段数值都出问题。Sentinel-2这边,我建议直接用大气校正后的地表反射率产品,也就是COPERNICUS/S2_SR;不要用COPERNICUS/S2那个L1C顶级产品,因为它给的是大气顶部反射率,做定量分析还得自己校正一遍,非常麻烦。如果遇到新的处理基线导致数值不连续的问题,可以考虑COPERNICUS/S2_SR_HARMONIZED,这个集合人为统一了不同处理基线的尺度,时间序列分析更稳。

Sentinel-1那边常用的是COPERNICUS/S1_GRD。这个产品已经做过轨道校正、热噪声去除、辐射定标和地形校正,拿到手可以直接用。但要特别注意,它存储的数值不是标准后向散射系数,官方说明是“需要乘以0.01才是线性σ°”。如果你不想手动处理这一步,可以改用较新的COPERNICUS/S1_GRD_FLOAT集合,那个数据直接就是线性σ°,省心很多。

3.2 尺度与投影:统一坐标系和分辨率

合并多波段之前最要命的一件事就是“网格对齐”。Sentinel-2的10米波段(B2、B3、B4、B8)和20米波段(B11、B12)本身分辨率就不一样,Sentinel-1 GRD产品虽然标称10米,但它的投影网格和S2的网格不是天然对齐的。如果直接把两个影像拼在一起而不做任何重采样,后续按像素提取样本时会出现很细微但真实存在的错位。

处理思路是选一个统一目标网格,一般以S1的投影为准,把S2的波段重采样到10米。具体做法是在GEE里对S2合成影像调用resample('bilinear'),然后reproject到S1网格的CRS。这里我要提醒一句,reproject本质上是在目标网格上重新采样,它会改变像素值,所以不要反反复复做,确定一次即可。

3.3 Sentinel-2去云与波段选择

S2影像合并前要做两类预处理:云掩膜和尺度缩放。云检测使用QA60波段里的云与卷云位掩膜,把云和卷云的像素置为无效。这里很多人会忽略卷云波段,但夏季影像里的薄卷云是最难用肉眼识别的,它会让反射率整体抬升,造成NDVI计算偏差,所以QA60的bit 10和bit 11必须同时处理。

波段选择上,做通用分类和参数反演时我用得最多的是B2(蓝)、B3(绿)、B4(红)、B8(近红外)、B11(短波红外1)、B12(短波红外2)这六个波段。前四个提供植被和土地覆盖的核心信息,后两个对植被含水量和土壤湿度很敏感,而且它们是区分许多地类的关键波段。如果对红边感兴趣,可以再加B5或B8A,但20米分辨率的红边波段参与分类前建议先重采样到10米。

3.4 Sentinel-1的去噪与极化波段处理

S1这边主要处理两个极化通道:VV和VH。选择IW成像模式的数据,同时必须确认影像同时包含这两个极化,很多GRD影像只有单一极化,过滤时要用listContains检查极化列表。

预处理的第一步是把原始值乘以0.01换算成线性σ°,然后做中值合成。噪声压制方面,单景雷达的斑点噪声很严重,表现为同质地物上的随机颗粒感;中值合成能压掉大部分,因为斑点噪声近似服从乘性分布,中位数估计比均值稳健得多。如果中值合成后还有明显的椒盐感,可以对VH波段做一次3×3的平滑,但注意不要反复平滑,否则会牺牲边缘信息。

还要留意一个容易被忽略的点:不同轨道方向(升轨/降轨)拍摄的同一地物,后向散射会有1到2dB的差异,因为入射角和方位角不同。做长时间序列分析时,尽量只保留同一轨道方向的数据,或者在合并前按入射角波段做一个简单的归一化。

4. 核心实操:在GEE里面完成S1与S2多波段合并

4.1 完整代码:同一个月S1+S2复合影像

下面这套代码是我日常用的基础版,逻辑很清晰:先分别构建S2和S1的月度中值合成,统一网格,然后拼接波段。复制到GEE编辑器里,把roi和日期换掉就能跑。

// 1. 确定研究区和时间范围 var roi = ee.Geometry.Point([110.2, 35.3]).buffer(5000); var start = '2023-10-01'; var end = '2023-11-01'; // 2. Sentinel-2地表反射率:去云、缩放、选波段 function maskS2clouds(image) { var qa = image.select('QA60'); var cloudBitMask = 1 << 10; var cirrusBitMask = 1 << 11; var mask = qa.bitwiseAnd(cloudBitMask).eq(0) .and(qa.bitwiseAnd(cirrusBitMask).eq(0)); return image.updateMask(mask).divide(10000); } var s2 = ee.ImageCollection('COPERNICUS/S2_SR') .filterBounds(roi) .filterDate(start, end) .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20)) .map(maskS2clouds) .select(['B2','B3','B4','B8','B11','B12']); var s2median = s2.median(); // 3. Sentinel-1 GRD:转σ°、取VV/VH、中值合成 var s1 = ee.ImageCollection('COPERNICUS/S1_GRD') .filterBounds(roi) .filterDate(start, end) .filter(ee.Filter.eq('instrumentMode', 'IW')) .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV')) .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH')) .map(function(image) { var vv = image.select('VV').multiply(0.01).rename('VV'); var vh = image.select('VH').multiply(0.01).rename('VH'); return ee.Image.cat([vv, vh]).copyProperties(image, ['system:time_start']); }); var s1vv = s1.select('VV').median().rename('VV'); var s1vh = s1.select('VH').median().rename('VH'); // 4. 统一网格:S2重采样到S1的10米网格 var s2resampled = s2median .select(['B2','B3','B4','B8','B11','B12']) .resample('bilinear') .reproject({crs: s1vv.projection(), scale: 10}); // 5. 多波段合并 var combo = s2resampled.addBands(s1vv).addBands(s1vh); // 6. 快速检查 Map.centerObject(roi, 12); Map.addLayer(combo, {bands: ['B4','B3','B2'], min: 0, max: 0.3}, 'S2 RGB'); Map.addLayer(combo.select('VV'), {min: 0, max: 0.25}, 'S1 VV'); Map.addLayer(combo.select('VH'), {min: 0, max: 0.15}, 'S1 VH'); print('Composite bands:', combo.bandNames());

这段代码里S2的合成用了median(),它会把多景影像的每一个波段取中位数,云掩膜后的无效像素会被剔除,因此合成结果对云残影有很好的抵抗力。S1同样用中值合成,既压斑点又保留边缘。最后addBands只是把两个影像的波段在通道维上拼起来,不会改变像素尺寸。

4.2 逐景配对版:保持时间同步的时间序列复合

如果你要做的是变化检测或者物候分析,月合成会掩盖细节,这时候就需要逐景配对。下面这段代码以每个S1影像的日期为锚点,前后各取10天的S2窗口,生成一组S1时刻的复合影像序列。

var s1List = s1.toList(100); var pairedList = s1List.map(function(item) { var img = ee.Image(item); var date = img.date(); var s2Window = s2.filterDate(date.advance(-10, 'day'), date.advance(10, 'day')); var s2Ref = s2Window.median().select(['B2','B3','B4','B8','B11','B12']); var radar = img.select(['VV','VH']); return ee.Image.cat([radar, s2Ref]) .set('system:time_start', img.date()) .set('date', date.format('YYYY-MM-dd')); }); var pairedCol = ee.ImageCollection(pairedList);

这段代码的运行效率没有月合成高,因为每个S1节点都要临时查一次S2集合,但它在时间精度上要严谨得多。要注意的是,当S2窗口内没有合格影像时,s2Ref会是空波段,生成的那一景复合影像就没有光学信息,后续应直接丢弃或做标记。

4.3 合并后的波段梳理与快速目视检查

合并完成后,建议第一时间打印波段名确认顺序。我这套基础版产出的复合影像一共有8个波段:B2、B3、B4、B8、B11、B12、VV、VH。波段顺序看起来无关紧要,但在后续做样本采集或分类输入时,波段顺序直接决定特征矩阵的列名,一旦错位很难排查。

目视检查也有技巧。先单独看S2的RGB组合,确认没有大片云残影;再单独加载S1的VV或VH,看地物边界是否清晰;最后把S1的VH波段作为红色通道、S2的NDVI作为绿色通道叠加显示,能快速发现两者之间有没有错位或异常。这个叠加画面里,植被区域应该是绿色主导,水体区域应该是红色偏暗,如果出现大量奇怪的橙色,多半是配准出了问题。

5. 实测踩坑记录:坑的成因、排查链路与对策

5.1 偏移一个像素:坐标系不一致引发的对齐问题

有一次我在一个山地试验区做复合影像验证,把S2的RGB和S1的VH叠在底图上,怎么都对不齐,沟谷边缘总是错开一个像素左右。我第一反应是雷达侧视造成的几何形变,但单独对比S1和参考底图,发现它本身并没有那么大的错位;再对比S2和参考底图,也只在陡坡处有轻微偏移。问题最后锁定在两者网格不一致上:S1 GRD使用它自带的投影网格,S2_SR使用的是另一个UTM分带的网格,GEE在动态渲染时会对两者分别重投影,而重投影参数不一致就产生了亚像元到像元级的错位。

排查链路是这样的:先分别显示两个单波段影像,确认各自的几何形态;然后在同一个放大级别下叠加,统计边缘错位方向;最后查看两个影像的projection()信息,确认CRS不同。解决方法是把S2重采样到S1的投影网格上,也就是前面代码里reproject那一步。这里再补一个建议:reproject的scale参数必须和S1的一致,否则你统一了坐标系却没统一分辨率,像素边界依然对不齐。

5.2 雷达波段权重失控:量纲混用导致分类结果一边倒

这是我早期分类任务中翻车最狠的一次。当时图省事,直接用原始GRD影像的DN值去做波段合并,S1的VV波段数值通常在几百到上千,S2反射率则在0到1之间,合并后训练随机森林,模型分裂点的选择几乎被S1两个波段主导,光学波段形同虚设。精度看似还行,但画出来的分类图呈现奇怪的条带纹理,完全不符合地类分布。

原因就是量纲不一致。后向散射和反射率的物理单位不同,数值尺度差了几个数量级,分类器计算特征重要性时方差大的波段天然占便宜。解决办法分两步:第一步把S1原始DN乘以0.01变成线性σ°,让量级回到0到0.3左右;第二步对合并后的所有波段做一个z-score标准化,也就是减去均值除以标准差。标准化这一步对距离度量类的算法尤其重要,随机森林虽然对单特征缩放不敏感,但多特征混用后你无法确定树分裂时的基尼增益会不会被某一个大尺度特征带偏,所以做总比不做好。

5.3 “假同步”:同月不同天的复合影像失真

我刚接触S1和S2连接分析时犯过一个认知错误:看到“同一个月”就放心了,直接用S1月度中值和S2月度中值拼成复合影像。有一次做作物收割监测,发现一块田的S1信号显示地表已经变得非常粗糙(像是收割后的茬地),而S2影像里的NDVI却还是绿色高值,两类信号明显冲突。

排查之后发现,那一年的S1月度合成只包含了月初的几景,而S2月度合成主要是月末的影像,两者之间隔着一次收割事件,地表状态早已完全不同。这就是典型的“假同步”问题。从那以后我定了一条规矩:民用产品做监测周期可以按月,但做变化检测必须用4.2节的逐景配对逻辑,以S1的成像时间为准来窗选S2,并且每次合成后都要检查system:time_start,确认S1和S2的日期差在允许范围内。

5.4 斑点噪声与椒盐状斑块

S1单景影像的斑点噪声很重,尤其是VH极化,数值动态范围里夹杂着大量颗粒噪声。如果不做中值合成,合并出来的复合影像在分类时会出现密密麻麻的椒盐状误分斑块,地物边缘也被噪声撕成锯齿。我踩过这个坑后的经验是:能中值合成就不单景硬拼;必须用单景时,至少对VV和VH做一次3×3的focalMean平滑。

还要注意一个问题,focalMean会轻微模糊边缘,如果后续要做地物边界提取,这一步要谨慎。我自己通常只对VH做平滑,VV保留原始信息,因为VH信号弱、噪声占比高,平滑收益更大。另外,在地形起伏大的区域,雷达阴影和叠掩区域会形成异常低值或高值,这些像素建议用坡度数据直接掩膜掉,否则它们会在合并影像里以假地物的身份参与分类。

5.5 轨道方向不同带来的信号断层

长时序复合影像还会遇到一个隐蔽问题:升降轨数据混用后,同一个地物的VV/VH信号在地块内部出现断层,像是被分成了两个不同的类别。这不是地表变了,而是雷达观测几何变了。我的处理办法是在数据过滤阶段就限定轨道方向,优先使用降轨数据,因为降轨覆盖范围通常更适合中低纬度区域;如果某段时间降轨数据不足,再补充升轨,但补充进来的数据我会在后续分析中单独标注,不让它进入同一批训练样本。

6. 复合影像的下游应用与导出细节

6.1 复合影像的下游分析方向

拼好这张8波段复合影像之后,能做的事就非常多了。做土地覆盖分类时,可以直接拿波段名称作为特征,加一个随机森林分类器训练;做作物识别时,可以在复合影像的基础上派生雷达植被指数RVI(4*VH/(VV+VH)),再配合S2的NDVI、NDWI一起进模型,效果比单用光学指数好不少;做洪水监测时,用淹没前后的VH差值影像叠加S2水体指数,能快速区分“积水区”和“过水区”。

这里有一个通用建议:不要把所有原始波段一股脑全喂给模型。S2的六个波段之间相关性很高,S1的VV和VH也密切相关,8个波段直接进模型会出现信息冗余。我一般会做一个特征筛选,保留B4、B8、B11、B12、VV、VH这六个主要特征,必要时再添NDVI、NDWI、RVI三个派生指数,一共九个特征,模型效果和计算效率最平衡。

6.2 导出参数与本地检查要点

复合影像最终要导出到本地做进一步分析,导出时最容易踩的坑是maxPixels限制。默认导出任务一次最多处理约1亿像素,10米分辨率的区域稍大一点就会报错。我习惯直接把maxPixels设置成1e13,并且指定scale: 10,同时手动指定一个合适的UTM投影CRS,避免GEE默认投影带来的分带混乱。

Export.image.toDrive({ image: combo, description: 'S1_S2_composite_202310', folder: 'GEE_exports', region: roi, scale: 10, crs: 'EPSG:32650', maxPixels: 1e13, fileFormat: 'GeoTIFF' });

导出的GeoTIFF在本地打开后,我建议做两件事。第一,检查波段顺序和数值范围,尤其确认S1波段不是0到几千的原始DN值,而是0到0.3左右的线性σ°;第二,叠加一张当地高分辨率影像做目视抽查,看S1的VH边缘是否和光学底图吻合。很多人在GEE里看着没问题,导出来才发现坐标系错了,这一步能用两分钟避免后面一整天的返工。

最后再分享一个我每次合并完都坚持做的检查动作:随机在图上选十个点,分别查看它们在每个波段上的采样值,并打印S1和S2各自的采集日期。十个点不多,但足以暴露“假同步”、值域异常和网格错位这三类最常见的问题。做连接分析最怕的不是处理步骤复杂,而是数据已经错了还在傻傻往下算,所以这个检查动作我强烈建议你也养成习惯。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询