简介面向遥感初学者与SNAP用户的Sentinel-1/2数据预处理图文指南聚焦SAR与光学影像的完整处理链路。文档以SNAP工具为主线详细讲解Sentinel-1的辐射定标、几何校正、斑点滤波、多视处理与投影输出同时针对Sentinel-2介绍了基于Sen2Cor 2.4.0的辐射定标、大气校正及重采样方法并附有Sen2Cor的下载安装说明。全流程按实际操作顺序组织从数据导入到结果导出环环相扣每个步骤均提供界面截图与参数设置示意尤其对辐射定标中的极化选择、地形校正中的DEM设置、滤波算法选取等易错环节做了明确演示可有效降低上手门槛。资源为单个docx文档约24.27MB内容组织清晰便于对照操作与后续查阅截至当前已有6429人学习下载适合需要系统掌握SNAP预处理流程、开展SAR或光学遥感应用研究的入门及进阶用户。1. 哨兵1号2号预处理S1的雷达信号和S2的光学反射率为什么都被SNAP收了第一次接触SARP的人最容易问一个问题S1是发射微波脉冲、靠回波成像的主动雷达S2是晒太阳、记录地表反射的阳光成像的光学传感器两者物理机理完全不一样为什么数据处理教程最后都指向同一款软件答案很简单处理链路不一样但工作流骨架一样——都需要从原始产品走到一个几何精校正、辐射一致、能直接算地学指数的栅格。SNAP的优势在于把S1的轨道校正、热噪声去除、辐射定标、地形校正和S2的大气校正、重采样、波段运算统一在同一个工程里中间结果不用在多个工具之间倒来倒去。这篇文章适合刚拿到S1 GRD或S2 L1C产品、准备做时序分析或光学与SAR联合的从业者我会按环境准备、两条预处理链路、常见坑和批量验证的顺序把这套流程讲透。2. 先把环境搭对SNAP版本、JVM内存与snappy初始化2.1 桌面版还是命令行版按使用场景选SNAP的图形界面适合做三件事单景数据参数调试、中间结果目视检查、拿着鼠标一步步理解每个算子到底改了图像的什么。比如第一次处理S1 GRD我建议你在图形界面里把Read、Apply-Orbit-File、Calibration几个算子逐个加进去每执行一步就看一眼图像直方图这样你对“热噪声去掉之后暗区抬升了多少”“定标之后数值从幅度变成后向散射系数”会有直观记忆。但图形界面的问题也很明显处理大景SLC或长时间序列时渲染每个中间结果都要消耗内存和显存而且参数不可复现换一台机器就得重新点一遍。命令行版叫gpt是SNAP安装目录下的一个可执行脚本。它读取一个Graph文件就是你在图形界面里搭好的算子链可以导出成XML再通过命令行参数注入输入输出路径。我的习惯是先在一景小数据上用桌面版把Graph调通导成XML之后所有批量任务都走gpt。这样处理和验证完全分离同一套Graph在不同日期、不同订户的产品上反复跑结果质量是一致的。2.2 JVM内存与Python接口snappy-conf一步到位SNAP本身是Java应用处理S1全分辨率GRD或S2的10m波段时默认JVM堆内存经常不够。很多人在图形界面里点“Run”之后发现进度条卡在70%然后整个程序无响应十有八九是内存没放开。改法很简单找到SNAP安装目录下的配置文件搜Xmx字符串默认一般是2g改成8g或者更大。# 先找到配置文件位置再修改内存参数 grep -r Xmx /opt/snap/etc/ # 把 -Xmx2g 改成 -Xmx8g然后重启SNAP生效配置完之后重点是把Python接口接上。SNAP的Python库叫snappy它本身不通过pip直接安装而是用SNAP自带的snappy-conf脚本生成绑定。这样做的目的是让snappy和你机器上的Python解释器精确对应避免系统里有多个Python版本时导入到错误的库。# 把snappy绑定到当前需要用的Python环境 /opt/snap/bin/snappy-conf /usr/bin/python3 # 验证是否绑定成功 python3 -c import snappy; print(snappy.__file__)注意snappy-conf的第一个参数必须是完整的Python解释器路径。如果你用虚拟环境就传虚拟环境里的python路径比如/home/user/venv/bin/python。绑定完成之后snappy就可以通过GPF.load读取Graph文件再用ProductIO读写产品这在后面做批量处理时会用到。如果这一步只报错不通过最常见的原因是Python版本太新或太旧SNAP每个大版本支持的Python范围不同先换成系统自带的Python再绑一次往往就好了。2.3 数据准备与目录组织SLC/GRD、L1C/L2A怎么放动手处理之前先搞清手里是什么产品。S1的Level-1产品分为SLC和GRDSLC保留了相位信息单视复数数据适合做干涉、形变、极化分解GRD是多视强度数据已经做了距离向和方位向的多视平均丢掉相位但噪声更小适合做后向散射分类、水体识别、时序振幅分析。如果你做的方向不是干涉老老实实用GRD处理速度快得多文件也小一个量级。S2这边L1C是大气顶反射率包含辐射校正和几何校正但没有做大气校正L2A是地表反射率由L1C经过大气校正得到同时附带云和云影掩膜。SNAP能直接读L1C和L2A但L1C必须先通过Sen2Cor或其他大气校正工具转成L2A否则算出的NDVI会被大气散射污染尤其蓝波段和近红外波段偏差明显。数据下载之后建议统一归档不要散落在各个下载目录里。我一般按卫星和级别分三层存放mkdir -p /data/rs_raw/s1_grd mkdir -p /data/rs_raw/s2_l1c mkdir -p /data/rs_l2a mkdir -p /data/rs_outS1的产品是一个zip压缩包SNAP可以直接读zip不需要解压但zip内的结构不要手动改动S2的L1C产品是一个SAFE目录目录里是各种xml和dat文件SNAP靠manifest.safe文件识别产品里面任何一级文件的改名都会导致读取失败。这一点在高版本SNAP上尤其严格很多人下载完顺手改了个带中文的目录名结果Read算子报错实际就是文件名破坏了。产品级别内容适合场景处理重点S1 GRD多视强度无相位水体/植被/地表分类、时序振幅热噪声、定标、地形校正S1 SLC单视复数含相位干涉、形变、极化配准、去平地、相位滤波S2 L1C大气顶反射率需要自己控制大气校正参数Sen2Cor或ATCORS2 L2A地表反射率云掩膜植被指数、土地利用、时序重采样、波段运算另外提醒一句处理前先看一眼数据的轨道方向Ascending还是Descending。S1的地形校正和S2的投影都依赖轨道和姿态信息如果你的研究区刚好在陡峭地形不同轨道方向叠出来的阴影方向完全不同这会影响后续分析的可靠性。3. Sentinel-1 GRD进SNAP轨道校正、热噪声、定标与地形校正的命令链3.1 算子处理顺序为什么一个都不能省S1 GRD的预处理在SNAP里是一条固定算子链顺序是Read → Apply-Orbit-File → ThermalNoiseRemoval → Calibration → Terrain-Correction → 可选的dB转换和裁剪。每步都有明确的物理目的顺序不能乱。轨道校正把产品自带的粗略轨道替换为精轨。轨道信息直接影响后续地形校正的几何精度尤其在山区轨道差几十米就能让叠出来的SAR图像错位好几个像元。热噪声去除针对GRD产品Sensor内部的噪声在暗区如平静水面会产生横向条带不除掉的话后向散射系数会被系统性抬高。定标把数字像元值换算为物理意义的后向散射系数选择输出Sigma0。地形校正解决的是SAR侧视成像带来的几何形变和辐射畸变把图像投影到标准地图坐标系同时用DEM做辐射归一化。最后一步是把线性功率值转成dB方便与其他数据的统计区间对比。3.2 从Read到Terrain-Correction一个可直接套用的Graph文件下面这个Graph是我在GRD IW模式下常用的完整链路双极化VV和VH同时输出中间用热噪声去除和定标最后地形校正输出10m分辨率GeoTIFF。?xml version1.0 encodingUTF-8? graph idS1_GRD_Preprocess version1.0/version node idRead operatorRead/operator sources/ parameters file${source}/file /parameters /node node idOrbit operatorApply-Orbit-File/operator sources sourceProduct refRead/ /sources parameters orbitTypeSentinel-1/orbitType polyDegree3/polyDegree /parameters /node node idThermalNoise operatorThermalNoiseRemoval/operator sources sourceProduct refOrbit/ /sources parameters/ /node node idCalibration operatorCalibration/operator sources sourceProduct refThermalNoise/ /sources parameters outputSigmaBandtrue/outputSigmaBand outputImageScaleInDbfalse/outputImageScaleInDb selectedPolarisations polarisationsVV/polarisations polarisationsVH/polarisations /selectedPolarisations /parameters /node node idTerrainCorrection operatorTerrain-Correction/operator sources sourceProduct refCalibration/ /sources parameters demNameSRTM 3Sec/demName imageResamplingMethodBILINEAR_INTERPOLATION/imageResamplingMethod pixelSizeInMeter10/pixelSizeInMeter mapProjectionUTM / WGS84 (Automatic)/mapProjection /parameters /node node iddBConv operatorBandMaths/operator sources sourceProduct refTerrainCorrection/ /sources parameters targetBands targetBand nameSigma0_VV_db/name expression10 * log10(Sigma0_VV 1e-6)/expression descriptionVV in dB/description /targetBand targetBand nameSigma0_VH_db/name expression10 * log10(Sigma0_VH 1e-6)/expression descriptionVH in dB/description /targetBand /targetBands /parameters /node node idWrite operatorWrite/operator sources sourceProduct refdBConv/ /sources parameters formatNameGeoTIFF/formatName file${outputDir}/S1_GRD_preprocessed.tif/file /parameters /node /graph保存为s1_grd_preprocess.xml后用gpt执行gpt s1_grd_preprocess.xml \ -Ssource/data/rs_raw/s1_grd/S1A_IW_GRDH_1SDV_*.zip \ -PoutputDir/data/rs_out整个Graph的逻辑是读取zip产品先用精轨替换原始轨道然后做热噪声去除这一步对GRD来说是必须的对SLC不是接着定标输出Sigma0线性功率再做地形校正把图像重采样到10m网格并投影到UTM最后用BandMaths生成两个dB波段。之所以在dB转换里给Sigma0_VV加了1e-6的极小值是为了防止暗目标回波接近0时计算出-inf这个值会在后续统计和可视化里污染整景影像的范围。3.3 关键参数表定标输出盯Sigma0地形校正盯分辨率Calibration算子的核心参数是输出类型和是否直接转dB。输出类型有Beta0、Sigma0、Gamma0三种Beta0是斜距方向的雷达亮度Sigma0是地面投影的后向散射系数Gamma0是按椭球面归一化的值。绝大多数地表参数反演和分类任务用Sigma0Gamma0在森林和粗糙度研究里用得多一些。我建议在线性域做定标dB转换放到后面这样遇到多视处理或滤波算子时数值平稳不会因为dB域的非线性放大噪声。算子参数推荐值说明Apply-Orbit-FileorbitTypeSentinel-1用精密轨道若下载失败自动用恢复轨道CalibrationoutputSigmaBandtrue输出Sigma0不输出Beta0/Gamma0CalibrationoutputImageScaleInDbfalse先保留线性功率dB放在最后Terrain-CorrectiondemNameSRTM 3Sec默认DEM平坦地区够用Terrain-CorrectionpixelSizeInMeter10与GRD原始分辨率匹配别填100Terrain-CorrectionimageResamplingMethodBILINEAR_INTERPOLATION比Cubic快对噪声数据更稳地形校正里的pixelSizeInMeter是我最常被问到的一个参数。GRD的标称分辨率是20m×22m多视处理之后实际像元地面尺寸约10m所以重采样到10m是正确的选择。有人为了“快一点”填50或者100结果图像像被磨了一层水体边界糊成一片后面对比S2 10m数据时完全对不上。另外如果研究区在高纬度地区UTM自动分带一般没问题但如果你的范围跨了两个投影带建议手动指定统一的EPSG代码否则多景拼接时会出现错位。4. Sentinel-2 L1C到L2A大气校正、重采样与波段运算的落地路径4.1 Sen2Cor的角色大气校正到底改了什么S2的L1C产品已经做了辐射定标和正射校正但像元值是大气顶反射率光线在进出大气层的过程中受到分子散射、气溶胶散射和水汽吸收的影响蓝波段被大气抬升得最明显近红外波段在植被区域被散射后数值偏低。直接拿L1C算NDVI植被区域会系统性偏低水体反射率也会失真。Sen2Cor是官方的大气校正处理器输入L1C的SAFE目录输出L2A地表反射率同时生成云、云影、雪和水的概率掩膜。L2A的反射率才真正适合做时序分析和跨时间对比。Sen2Cor有两种使用方式一种是SNAP图形界面里的集成入口另一种是独立安装命令行版本。我推荐命令行独立版本因为输出路径可控而且可以和gpt的Graph无缝衔接。4.2 从L1C到L2ASen2Cor命令与参数登录服务器后进入Sen2Cor的bin目录执行L2A处理cd /opt/sen2cor/bin # 把L1C转成L2A输出10m分辨率地表反射率 ./L2A_Process --resolution 10 /data/rs_raw/s2_l1c/S2A_MSIL1C_*.SAFE这条命令里的--resolution有三个可选值10、20、60。选10表示输出波段统一为10m网格其中原生20m和60m波段也会被重采样到10m选20则只保留原生20m和60m波段重采样到20m选60速度最快但红边波段等细节会丢失。我一般分任务取舍做NDVI或者水体提取用10做土地覆盖分类用20快速预览用60。执行过程中会在终端输出气溶胶反演、水汽反演、场景分类等阶段中间如果报“ERROR”不要慌先看是气溶胶厚度的错误还是DEM缺失的错误。Sen2Cor处理完成后在L1C目录的同级会生成一个以S2A_MSIL2A_开头的SAFE目录这就是后续SNAP读取的输入。4.3 重采样与NDVI计算一个可直接复用的GraphL2A产品里波段分辨率不统一B2、B3、B4、B8是10mB5、B6、B7、B8A、B11、B12是20mB1、B9是60m。直接做波段运算会因网格不一致报错所以第一步用Resample算子统一到B8的10m网格。下面这个Graph读入L2A重采样后计算NDVI并输出GeoTIFF?xml version1.0 encodingUTF-8? graph idS2_L2A_NDVI version1.0/version node idRead operatorRead/operator sources/ parameters file${source}/file /parameters /node node idResample operatorResample/operator sources sourceProduct refRead/ /sources parameters referenceBandB8/referenceBand /parameters /node node idNDVI operatorBandMaths/operator sources sourceProduct refResample/ /sources parameters targetBands targetBand nameNDVI/name expression(B8 - B4) / (B8 B4)/expression descriptionNormalized Difference Vegetation Index/description /targetBand /targetBands /parameters /node node idWrite operatorWrite/operator sources sourceProduct refNDVI/ /sources parameters formatNameGeoTIFF/formatName file${outputDir}/S2_NDVI.tif/file /parameters /node /graph保存执行gpt s2_ndvi.xml \ -Ssource/data/rs_l2a/S2A_MSIL2A_*.SAFE \ -PoutputDir/data/rs_out这里Resample只用referenceBand指定参考波段B8所有波段自动重采样到B8的10m网格省去手动填像元尺寸的麻烦。BandMaths表达式里的B8对应近红外波段B4对应红波段如果SNAP报“B8 not found”多半是表达式里的波段名带引号或数据集命名不同在图形界面里打开BandMaths从可用波段列表里双击插入波段名就不会错。NDVI输出后建议顺手在S2的L2A产品里找一下场景分类图层SCL它包含了云、云影、水的掩膜可以用来过滤NDVI时序中受云影响的像元。4.4 裁剪与统一投影多时相叠图的习惯做法S2的L2A产品覆盖范围约290km宽但很多研究区只有一小块全图处理既慢又浪费磁盘空间。S1地形校正之后的输出坐标是UTMS2原生也接近UTM但两者自动投影可能不在同一个投影带。做时序分析或S1/S2联合时我会先确定一个基准投影通常选择研究区中心所在的UTM带然后所有产品的Graph尾部加Spatial-Subset算子裁剪到统一范围。裁剪有两种方式按经纬度范围直接裁剪或者用一个矢量shp文件裁剪。按范围裁剪时我会用S2 L1C产品在SNAP里打开目视确定研究区范围后记下坐标再写进Graph的Subset参数。按矢量裁剪更精确先Import-Vector导入shp再在Subset里选择“使用矢量裁剪”这样输出产品与研究区边界完全一致。多时相叠图最怕的是每景的投影和范围不统一宁可处理慢一点也要把这一步在流程里固定下来。5. 预处理避坑指南从现象到原因的5条排查记录5.1 热噪声去除顺序错乱Sigma0整体偏高且暗区有条纹现象处理完的GRD图像在水体等暗区看到明显的横向条带后向散射系数的直方图整体偏右水面的平均Sigma0比文献值高了2~3dB。原因热噪声去除被放到了定标之后而定标算子把所有像元统一乘以定标因子暗区的系统噪声也被当作有效回波放大了后续地形校正无法消除这种条带。解决严格按Read → Apply-Orbit-File → ThermalNoiseRemoval → Calibration的顺序执行热噪声去除一定要在定标之前因为噪声分量是加性的先在原始强度域里减掉噪声再定标才是正确的物理顺序。5.2 JVM内存不足gpt跑一半直接卡死现象gpt处理S1 GRD全分辨率时日志停在Terrain-Correction节点CPU占用不高但内存接近占满然后进程挂起kill都费劲。原因SRTM DEM和S1图像都在内存里做重采样默认JVM堆内存不够GC频繁导致假死。解决先在配置文件里把-Xmx改到8g甚至更高再执行处理如果单景范围特别大可以先把S1产品用Subset裁到研究区再进处理链路内存占用可以降一半以上。5.3 Sen2Cor在沿海场景报错或输出条带现象L2A处理到气溶胶反演阶段直接报错或者处理完成后10m波段的沿海区域出现东西向条带。原因Sen2Cor的暗像元算法在均匀水体地区反演不出足够的空间变异性沿海场景常出现气溶胶厚度无解或解不稳定。解决先用--sc_only参数只跑场景分类生成云掩膜后跳过气溶胶依赖的波段或者干脆改用SNAP里自带的替代大气校正工具如果只是做NDVI这种比值指数L1C的大气误差在比值里会大幅抵消应急时也可以直接用L1C算但要说明局限性。5.4 地形校正分辨率乱填处理时间暴涨且结果模糊现象地形校正设置了pixelSizeInMeter100处理速度确实快但输出的地物边缘像被磨砂过和同一时期的S2影像叠加时明显对不齐反过来填成2m处理时间从10分钟涨到2小时图像却并没有变清晰。原因pixelSizeInMeter是输出像素尺寸不是处理精度改得比原始分辨率更细只是插值加密不会增加真实信息填得过粗则丢失真实纹理。解决查一下你手里GRD产品的标称地面分辨率通常IW模式填10m或15m别只看处理速度。5.5 S1和S2投影不一致叠加显示错位现象S1处理完的sigma0和S2的RGB图像在SNAP里叠加河流边界错位10~20个像元但单独看每张图都在正确位置。原因S1的地形校正选了自动UTM带S2产品可能用了另一个投影带或S2的L2A在Sen2Cor处理时坐标基准与SNAP默认不一致。解决两个Graph里都显式固定同一个EPSG代码比如研究区在UTM 49N就写死32649不要用Automatic最后在SNAP里把两景图同时打开用Geocoding信息检查像元坐标是否一致如果仍有偏移优先检查DEM和轨道文件是否正常更新。6. 批量跑S1/S2预处理的脚本写法与结果验证技巧6.1 用bash循环把gpt串起来单景跑通过之后批量处理的关键是保持Graph不动只替换输入输出路径。Graph里的${source}和${outputDir}就是为这个设计的。下面这段脚本会把/data/rs_raw/s1_grd下所有S1 zip逐景处理输出目录按输入文件名命名方便后续追溯mkdir -p /data/rs_out for f in /data/rs_raw/s1_grd/*.zip; do name$(basename $f .zip) echo Processing $name gpt s1_grd_preprocess.xml \ -Ssource$f \ -PoutputDir/data/rs_out/$name done如果你的机器内存足够可以把for循环改成并行执行比如一次跑两个gpt进程。但我不建议并行超过两个因为SNAP的JVM内存是独立分配的两个进程各占8g的话16g内存的机器会立刻吃紧。串行处理虽然慢但稳定日志也容易排查。S2的批量处理思路一样先写一个循环跑Sen2Cor把L1C全转成L2A再用gpt的NDVI Graph批量产出指数图。6.2 结果怎么算“处理成功”三个区间一张掩膜预处理是不是成功不是看gpt有没有返回0而是看输出产品的统计量合不合理。我习惯打开SNAP用查看元数据和直方图来验证三件事。第一件事是S1的Sigma0_VV_dB值域。典型地物的区间大体是平静水体在-20dB到-10dB之间城市建成区在-5dB到0dB之间裸地和低矮植被在-10dB到-5dB之间。如果整个产品的最小值低于-30dB先检查是不是dB转换时出现了-inf那就是前面说的加极小值那步没做对如果城市区域的值高于0dB定标可能选成了Beta0而不是Sigma0。第二件事是S2的NDVI区间。水体NDVI一般在0以下裸地在0到0.2茂密植被在0.6以上。如果整景NDVI最大值不超过0.4优先怀疑大气校正没成功或者Sen2Cor输出的L2A波段名被BandMaths表达式引错了。另外一个常用验证是打开真彩色合成B4、B3、B2看植被的绿色是否自然如果图像发蓝或发灰大概率是大气校正失败。第三件事是几何一致性。把S1的sigma0和S2的NDVI在SNAP里用同一个EPSG叠加误差超过一个像元就要查轨道和DEM山区尤其明显地形校正后山谷和山脊的错位会直接暴露你的DEM分辨率是否够用。我还有一个土办法取研究区里一根笔直的公路或河堤分别在两景图像上量同一段长度差值超过5%就说明几何处理还有问题。最后分景输出一个检查记录把每景的统计值和当时用的DEM、轨道文件写进去后面做时序分析出异常时能快速定位是哪一景的数据质量引起的而不是怀疑算法。这一点是我踩过最深的坑以前某项目里我一口气批量处理了30多景S1没做逐景验证后面做洪水淹没范围提取时发现其中几景的水体面积系统性偏大查了一天才发现是这几景当时轨道文件没更新到精轨地形校正后坡度区域的值全被抬高了。从那以后我不管批量跑多少景每景都保留一个统计量清单。数据处理这事儿图省事后面往往要花几倍时间还回去。希望帮到你。本文还有配套的精品资源点击获取
