1. 先把 Sentinel-3 的四套载荷摸清楚接手这个项目的第一天我做的第一件事不是写代码而是把产品手册打印出来贴在工位上。原因很简单Sentinel-3 不是一颗“拍照片”的卫星它是一个多载荷、多产品级别的对地观测系统同一个轨道数据在命名、结构、坐标系上都有完全不同的样子。你要是按处理普通卫星影像的直觉去理解它大概率会在第三天崩溃。这套系统的核心价值在于业务化的延续性观测单星重复周期约 27 天双星编队后赤道地区重访缩短到一天左右过境时间被锁定在上午光照条件稳定。这两点决定了它和“看到就拍”的商业卫星完全不同——它追求的是同一地点、同一时刻、同一角度、连续多年可比。做过长期时间序列的人都知道“可比”两个字有多值钱它意味着你不需要每年重新标定一次算法阈值。1.1 一颗星上的四种“眼睛”分别负责什么把四套主力载荷拆开看每个都有自己的脾气。OLCI海陆色仪是使用频率最高的一个21 个波段覆盖 400–1020 nm空间分辨率 300 米幅宽 1270 公里推扫式成像。它的波段排布不是随便定的Oa08约 665 nm和 Oa11约 708 nm这对组合是水色反演的主力前者对应叶绿素吸收谷后者对应反射峰Oa11 还夹在荧光基线位置上所以它同时具备一定的藻类荧光探测能力。更关键的是 Oa13约 761 nm和 Oa14约 767 nm这两个氧 A 吸收带波段专门用来做大气校正中的云和气溶胶检测——这是设计时就把“校正”需求前置了不是事后补的。SLSTR海陆表面温度辐射计才是这个系统里最被低估的部分。它最特别的地方是双视角观测同一目标在天底方向和约 55° 倾斜方向各看一次两个路径的大气厚度不同联立起来就能同时反演地表温度和大气水汽。这个设计的直接收益是海表温度精度目标做到 0.3 K 量级且在潮湿热带海域不容易系统性偏移。它有 9 个通道可见/短波红外 500 米热红外 1 公里。做高温目标监测比如火点时短波红外的饱和特性要提前评估不然强火点直接削顶。SRAL合成孔径雷达高度计用 Ku 和 C 两个频段同时工作Ku 频段测距精度可以做到厘米级。它输出的核心量是海面高度、有效波高、后向散射系数。这个载荷和成像类载荷的思维方式完全不一样——它是一条“沿轨剖面”每秒 20 个测量点每个点是一个波形。你要处理的是时间序列不是图像。很多做遥感出身的人第一次碰高度计会不适应因为它没有二维结构只有一条线。MWR微波辐射计只有两个频段约 18.7 和 23.8 GHz任务单一测大气水汽含量给高度计的湿对流层延迟做校正。它不起眼但没有它海面高度的精度会掉一个数量级。提示这四套载荷的时间同步性很好但空间覆盖差别极大。OLCI 幅宽 1270 kmSLSTR 天底幅宽 1420 km、倾斜视角只有 740 kmSRAL 是星下点窄带。做多源融合之前先算一遍覆盖交集别默认“同一景数据覆盖同一区域”。1.2 产品分级与打包格式SEN3 目录到底装了什么下载下来的是一个以.SEN3结尾的目录不是单个文件。这个细节坑过很多人——用文件选择对话框直接选它会失败必须整目录拖走或者用递归传输。目录里通常包含一个 XML 清单文件记录元数据、处理基线、每个测量文件的校验值、一堆.nc格式的测量文件每个波段或每个地球物理量单独一个、坐标文件、时间文件、质量控制标志文件以及仪器辅助数据。这种“一个量一个文件”的设计对按需读取很友好——你只要叶绿素就不用把 21 个波段的反射率全读进内存。产品级别上能直接用的主要是两级1 级是经过辐射定标和几何定位的辐亮度/亮温物理意义明确但没有地球物理含义2 级是经过大气校正后的地球物理量比如叶绿素浓度、海表温度、气溶胶光学厚度。中间还有一个 1A 级主要是原始计数加部分校准日常应用基本不碰。真正需要记住的是后缀。同一颗星同一个载荷不同的两个字母代表完全不同的产品后缀含义分辨率典型用途EFR全分辨率辐亮度300 m自定义算法、高分辨率水体ERR降分辨率辐亮度1.2 km大范围快速浏览、云检测WFR水体全分辨率地球物理量300 m叶绿素、悬浮物、水体透明度WRR水体降分辨率1.2 km全球尺度长时间序列LFR陆地全分辨率300 m陆地反射率、植被指数RBT辐射亮温含双视角500 m / 1 km地表温度、火点WST海表温度1 km海洋温度场LST陆表温度1 km城市热环境、农业FRP火点辐射功率1 km火情监测WAT测高海洋产品沿轨 20 Hz海面高度、波高我见过最典型的误用是拿 WFR 去做内陆浑浊水体的叶绿素反演然后抱怨结果全是零。原因后面会讲这不是数据的问题是算法假设的前提被打破了。1.3 它和哨兵2号、哨兵1号的分工边界经常有人问“我有哨兵2号了还要不要用这个”。答案取决于你要颜色还是要光谱细节。哨兵2号是 10–60 米分辨率、13 个波段、5 天重访擅长陆地细节但它没有专门的氧吸收带波段也没有热红外通道水色反演的波段配置也不是为水体优化的。哨兵3号是 300 米看起来粗但它的波段是为水体和大气专门排布的21 个窄波段能支撑更精细的大气校正和生物光学反演同时还有热红外和海面高度。用一句话概括分工看陆地细节找 2 号算水体参数、测温度和测高找 3 号。真正做海岸带或大型湖泊研究的团队通常是 3 号做长时间序列的“面”2 号做关键时段的“点”两个拼起来用。2. 从需求反推数据选型你要的到底是哪个产品数据选型这一步做错后面所有的代码都是白写。我自己的经验是先回答三个问题我要的物理量是什么、我的目标尺度是多大、我要多久的时间跨度。这三个答案基本就能锁定唯一的载荷和产品级别。2.1 五类典型需求与最小可用产品清单把常见的需求列成表比翻手册快得多需求场景首选载荷最小可用产品需要额外准备开阔海域叶绿素/初级生产力OLCI2 级 WFR 或 WRR无直接读内陆湖泊浑浊度、藻华OLCI1 级 EFR自己算自定义指数 现场验证近岸海表温度SLSTR2 级 WST质量标志筛选陆地地表温度、热岛SLSTR2 级 LST比辐射率假设、地形校正大尺度海面高度与波高SRAL2 级 WAT轨道筛选、重跟踪注意“内陆湖泊”那一行的处理方式这是新手最容易搞错的地方。2 级水体产品默认是按一类水体大洋型设计的它的气溶胶模型和生物光学算法都假设水体光学性质由浮游植物主导。内陆水体里悬浮泥沙和有色溶解有机物占主导这套假设直接失效结果通常表现为浓度被严重低估、异常值大片出现。稳妥的做法是回到 1 级反射率自己做大气校正或者用专门的二类水体处理器再套自定义指数。2.2 文件名就是元数据一次读懂命名规则这个系统的文件名长得离谱但每一个下划线分隔的字段都有意义。举个典型结构S3A_OL_2_WFR____20240115T024000_20240115T024300_20240116T071200_0179_108_123______MAR_O_NT_002.SEN3拆开读S3A平台标识A 星或 B 星。做长时间序列必须区分两颗星的定标差异虽然很小但在趋势分析里必须记录在案。OL载荷代码OL 是色仪SL 是温度辐射计SR 是高度计SY 是协同产品。2产品级别。WFR产品类型。两个时间观测起止时间。这个是筛选数据的关键别拿下载时间当观测时间。第三个时间处理生成时间。0179相对轨道号同一个重访轨道周期内序号固定。108绝对轨道号连续递增。123圈内计数。MAR处理中心标识。O产品类别业务化/测试等。NT时效性等级。002处理基线版本。这一项极其重要不同基线版本的算法、定标系数、标志位定义都可能变混用会造成时间序列跳变。我自己踩过的坑把不同处理基线的产品直接拼成十年序列画出来的趋势线在版本切换那一年有一个明显的台阶看起来像生态突变实际上只是算法升级。注意如果你的分析跨越多年一定要在元数据里保留处理基线字段并在做趋势拟合前检查同一条时间线上是否混了版本。2.3 获取路径的取舍滚动存档、开放数据桶与本地缓存获取渠道大致三类各有适用场景。滚动存档适合“临时取几景”的需求界面友好、按时间区域检索但大批量下载会受单 IP 并发限制且没有断点续传保证。我的做法是只在探索阶段用它确认数据能覆盖研究区之后就切走。云端开放数据桶适合批量和重复计算。它的优势是免下载可以在数据所在的地方直接算把结果几 MB 的 GeoTIFF拉回本地而不是把几个 TB 的原始包搬一遍。做全球尺度分析时这一条几乎决定项目成败。代价是需要写一点对接代码并且要习惯对象存储的“按请求计费、内网访问便宜”的逻辑。本地镜像适合长期固定的研究区。把目标区域三五年的数据一次性拉到本地阵列后续所有实验都是本地 IO迭代速度最快。存储估算要提前做一景全分辨率色仪数据碰巧接近 1 GB一年的双星数据是数千景量级几个 TB 的阵列是起步价。三种方式不是互斥的。我目前的组合是用云端做全球筛查和快速验证用本地镜像做研究区的反复调试用滚动存档补零星缺失的个别时相。3. 读取与预处理netCDF 结构的拆解实操这一章是全文的技术核心。前面都是选型这里开始动手。3.1 xarray 打开 SEN3 目录的几个坑先说最基础的一件事这个目录里的每个.nc文件都是一个独立的 NetCDF 数据集xarray.open_dataset一次只能打开一个。所以你不能直接对.SEN3目录调用它必须精确指定内部文件。import xarray as xr import numpy as np base S3A_OL_2_WFR____20240115T024000_20240115T024300_20240116T071200_0179_108_123______MAR_O_NT_002.SEN3 # 坐标是每个像素都有的必须读 geo xr.open_dataset(f{base}/geo_coordinates.nc) lat geo[latitude].values lon geo[longitude].values # 叶绿素注意变量名和文件名要以实际产品为准先 ls 一遍 chl xr.open_dataset(f{base}/chl_oc4me.nc) print(chl.data_vars) # 先看有什么别猜三个细节容易被忽略第一坐标是逐像素数组不是一维轴。latitude和longitude的维度都是(rows, columns)形如 4000×4865。这意味着它不是一个规则网格而是一个“不规则采样点集合”。后面所有的重投影都要围绕这个事实展开。第二有效值掩膜。产品里通常有一个_FillValue或missing_value属性噪声像元会被填成一个极大值比如 9.99e8 之类。如果不先处理就参与统计你的均值会被拉到天上。稳妥做法是读进来立刻替换为np.nan。第三行列顺序和图像直觉相反。数组的第一个维度是沿轨方向的扫描行第二个是跨轨方向的列。有些产品在存储时做了转置用shape和总像素数交叉验证一下别直接假设。3.2 非规则网格为什么不能直接当 GeoTIFF 用这是新手最容易翻车的地方。你把chl数组和lat/lon数组一起丢给常规栅格工具它会报错或者生成一张完全错位的图因为它没有仿射变换参数。成像类卫星的常规产品带GeoTransform告诉你左上角坐标、像素大小、旋转项所以任意像素的坐标可以用一个线性公式算出来。这个产品没有这个东西它给的是每个像素的真实经纬度而且因为推扫式成像、地球自转、轨道弯曲等因素行与行之间的几何关系并不完全均匀。为什么这么设计因为它保真。推扫式传感器在扫描带两端会有明显的几何拉伸用仿射模型近似会引入系统性位置误差。逐像素给坐标虽然占存储两个 float32 数组4000×4865 就是 300 MB 出头但保证了定位精度。对应的处理逻辑就变了必须做重采样重投影把不规则采样点映射到你指定的规则网格上。这一步不做后面没法做多时相叠加也没法和别的数据源对齐。3.3 三种重投影方案与选型对比实际可用的路子有三条我按踩坑次数从多到少排方案工具优点缺点适用场景图形化处理器重投影桌面处理软件 命令行批处理内置支持好参数少启动慢批处理脚本调试麻烦一次性处理、可视化探索Python 库重采样开源重采样库灵活、可编程、能嵌进管道大数据量要自己管分块生产管线、云端计算先转 VRT 再转栅格栅格工具链生态成熟后续处理方便中间文件大IO 压力大需要与已有栅格工作流对接我目前的默认选择是第二条。原因很实在分块和并行可以自己控制而且能直接接上对象存储的流式读取避免中间文件。图形化处理器那套其实也很好用尤其是命令行模式一条命令就能出结果gpt reproject.xml -Ssourceinput.SEN3 -Poutputout_wgs84.tif -c 4G -q 8其中reproject.xml是操作图定义核心参数就三个目标坐标系、重采样方法、输出像元大小。像元大小怎么定给你一个可用的换算方法目标投影下 1 度纬度约 111 公里要把 300 米分辨率大致保住像元取 0.003 度左右约 333 米。如果做全球长时间序列很多人会降到 0.01 度约 1.1 公里来换存储和速度这时候用降分辨率产品反而更合适不要用全分辨率产品硬降。重采样方法上分类量用最近邻连续量用双线性。叶绿素浓度这种跨数量级的量双线性插值在浓度梯度大的锋面附近会产生虚假的中间值我个人更倾向于最近邻海表温度这种空间连续的场双线性更平滑标志位数组绝对只能用最近邻插值出来的标志位是没有物理意义的。3.4 质量控制标志位与云掩膜所有地球物理量产品都带一个质量标志数组通常是uint32的位掩码bitmask每一位代表一种状况。用法是按位与flags xr.open_dataset(f{base}/wqsf.nc)[WQSF].values # 示意具体位定义必须查对应产品手册的版本说明 CLOUD_BIT 1 LAND_BIT 0 INVALID_BIT 25 cloud_mask (flags (1 CLOUD_BIT)) ! 0 land_mask (flags (1 LAND_BIT)) ! 0 invalid (flags (1 INVALID_BIT)) ! 0 valid ~(cloud_mask | land_mask | invalid)这里必须强调位定义会随处理基线变化。我在不同年份的产品里见过同一含义的位被挪到不同位置的情况。所以你的代码里不应该写死位号而应该写一个“版本到位号”的映射表或者至少在管道里加一步断言检查标志位数组的取值范围和分布是否异常。经验值上一景开阔海域数据里云像元比例通常在 20%–60% 之间如果算出来云占比 99%那八成是位读错了不是你运气差。还有一个常被忽略的点云的“膨胀”。位掩膜给出的是被判为云的像元但云的边缘往往有半透明的薄云和阴影这些像元可能没被标为云数值却已经被污染。稳妥做法是对云掩膜做几次形态学膨胀腐蚀掩膜扩大云的判定范围代价是损失一些有效像元收益是减少异常值。我在做内陆水体时一般膨胀 2–3 个像元开阔海域膨胀 1 个就够。4. 典型应用实操海色、海温与海面高度前面铺垫完了这章讲三个最主流的用法每个都给出可直接复现的思路。4.1 水体叶绿素与浑浊度从标准算法到自定义指数标准产品里给的叶绿素用的是最大波段比值法简而言之就是取几个波段反射率的比值组合再套一个多项式。它的设计前提是一类水体光学性质主要由浮游植物及其降解产物决定其他成分可以近似为常数。在内陆和近岸这个前提碎了。悬浮泥沙会让整体反射率抬升有色溶解有机物会强烈吸收短波蓝光两者都会把比值算歪。表现出来就是明明肉眼可见浓绿色藻华叶绿素反演值却只有个位数。我的处理路径是回到 1 级反射率自己做两个东西第一先做归一化水体指数把水陆分开。用近红外波段的低反射特性做阈值分割这个步骤简单但必须做不然陆地上的植被会被算成高叶绿素。第二用双波段比值指数代替标准算法。比如用约 665 nm 和约 708 nm 构造归一化差值指数def ndci(r665, r708): num r708 - r665 den r708 r665 with np.errstate(divideignore, invalidignore): out num / den return np.where(np.abs(den) 1e-6, np.nan, out)这个指数在富营养化水体里和叶绿素的相关性通常明显好于标准算法代价是它没有物理单位只能做相对分级不能直接报浓度。要用浓度就得有现场采样做回归标定通常准备 20–40 个同步采样点覆盖低中高浓度梯度拟合一条经验曲线。别指望别人的系数能直接搬过来——不同湖泊的悬浮物和溶解有机物组成差别很大系数搬过来误差能到几倍。第三做基线插值指数识别表层藻华。取约 665 nm、708 nm、753 nm 三个波段用前后两个波段在 708 nm 处做线性插值得到一条“无藻华时的基线反射率”实际反射率高出基线的部分就是藻华信号。这个方法的优势是对薄云和观测几何不敏感适合大范围筛查。提示所有基于反射率的指数都必须先完成大气校正。1 级数据是大气顶的辐亮度直接算比值意义不大。如果暂时没有好的大气校正方案退而求其次的做法是用两个邻近波段的比值因为气溶胶影响在两个邻近波段上高度相关比值能部分抵消——但这是妥协不是正解。4.2 海表温度双视角的价值在哪海表温度产品的核心难点从来不是测温度是剥离大气。水汽在热红外波段的吸收很强而海面水汽含量变化又大单视角观测无法同时确定“地表温度”和“大气水汽”两个未知数。双视角就是为这个设计的。同一点被看两次两次穿过的大气路径长度不同吸收量不同两个方程解两个未知数。理论上的收益是水汽含量高的热带海域单视角产品的偏差能到 1 K 以上双视角可以压到 0.4 K 以内。实际使用时的几个注意事项云检测宁严勿宽。热红外对亚像元云极其敏感一个 10% 覆盖的薄云就能让温度偏低好几度。产品的质量标志里通常有置信度分级我一般只保留最高等级。近岸要打问号。1 公里像元在近岸会混入陆地信号尤其是潮间带和浅水区。这不是算法能解决的只能靠掩膜剔除或者用更高分辨率数据补。白天和夜间产品要分开看。白天的太阳反射会引入偏差很多应用场景比如日较差分析必须按地方时筛选观测。验证不要用陆地气象站。这是一个很基础的错误但经常出现拿岸边气象站的气温去验证海表温度两者物理意义完全不同一个在陆地 2 米高度一个是水面皮层温度。合理的做法是用浮标数据且要匹配观测时间窗口一般取 ±2 小时内。4.3 测高波形重跟踪与轨道级时间序列高度计的处理思路和成像完全不同我单独拿出来讲因为这块的坑最隐蔽。基本几何关系是卫星到海面的距离由雷达往返时间测得卫星轨道高度已知两者相减得到海面相对参考椭球的高度。这个高度里混合了大量信号大地水准面起伏、潮汐、海况偏差、对流层和电离层延迟。2 级产品里这些校正项基本都提供好了你要做的是选择性地应用并检查每一项的合理性。主要校正项我自己列了一个核对顺序干对流层由大气压算、湿对流层由微波辐射计算这是那个双频辐射计的作用、电离层双频雷达可以自校正、逆气压效应、固体潮、海洋潮、极潮、海况偏差。每一项都是厘米到分米量级海面高度异常信号本身也就十几厘米所以少一项误差就等于把信号淹了。重跟踪是绕不过去的一步。在开阔大洋标准的波形拟合算法足够好但在近岸、海冰边缘、内陆水域回波波形会被陆地或冰面污染标准算法给出的距离偏差可达几十厘米。解决办法是用专门的重跟踪算法重新估计波形前沿。目前的实践是用现成的重跟踪器参数上用窄的跟踪窗口和专门针对近岸的拟合模型。时间序列有一个关键技巧按轨道号筛选不要按经纬度框筛选。高度计的采样是沿轨稀疏的如果你用一个地理范围去裁不同重访周期裁出来的采样点位置会漂移导致时间序列里混进了空间差异。正确做法是固定相对轨道号这样每次过境走的是几乎完全相同的星下点轨迹空间位置一致差异才真的是时间差异。4.4 一个能跑的最小流程把上面的东西串起来给一段可以直接改参数用的骨架代码。我自己写管线时习惯用它起步验证数据读取和掩膜是否正常再往里塞业务逻辑。import numpy as np import xarray as xr import rioxarray from pyproj import CRS def load_olci_l2(base_dir, var_filechl_oc4me.nc, var_nameCHL_OC4ME): 读取 OLCI 2 级产品返回经纬度和目标变量缺失值统一成 nan geo xr.open_dataset(f{base_dir}/geo_coordinates.nc) lat geo[latitude].values lon geo[longitude].values ds xr.open_dataset(f{base_dir}/{var_file}) raw ds[var_name].values # 关键一步把填充值换成 nan fill ds[var_name].attrs.get(_FillValue, 9.99e8) data np.where(raw fill * 0.9, np.nan, raw) return lat, lon, data def make_valid_mask(base_dir, exclude_bits(0, 1, 25)): 按位排除陆地、云、无效像元bit 定义需按产品版本核对 flags xr.open_dataset(f{base_dir}/wqsf.nc)[WQSF].values bad np.zeros_like(flags, dtypebool) for b in exclude_bits: bad | (flags (1 b)) ! 0 return ~bad def build_point_table(lat, lon, data, valid, step4): 降采样抽稀成点位表交给后续 GIS 或数据库处理 sub (slice(None, None, step), slice(None, None, step)) lat, lon, data, valid lat[sub], lon[sub], data[sub], valid[sub] ok valid np.isfinite(data) (lat -89.9) (lat 89.9) return lon[ok], lat[ok], data[ok]step4这个抽稀参数不是随手定的。300 米分辨率的全分辨率产品一景有效像元可能有几百万个直接做点密度图会把内存吃穿。降采样 4 倍后像元间距约 1.2 公里对做区域态势图完全够用内存压力降一个数量级。如果你要做逐像元的精细统计那就别抽稀改成按块流式写入文件。5. 批处理工程化与性能优化单景处理跑通之后真正的难题才开始你有几百上千景数据怎么让它在可接受的时间里跑完。5.1 内存、分块与并行一景全分辨率色仪的反射率21 个波段全部读进来按 4000×4865 的尺寸算是 800 MB 量级加上坐标、标志位、中间结果内存峰值轻松破 3 GB。看着不大但如果你用多进程跑 8 个并发就是 25 GB普通工作站直接被打爆。三个可用的手段分块读取。用支持懒加载的数组结构把数据切成 1024×1024 的块处理完一块写一块内存占用就固定在一个块的大小上。这个改动通常能把内存峰值降到几百兆。限制并发数而不是死磕单进程速度。我的经验值是按“单景峰值内存 × 并发数 ≤ 物理内存的 70%”倒推并发数剩下的 30% 留给系统缓存和垃圾回收。宁可少开两个进程也不要让系统开始换页——一旦开始换页整体速度会慢十倍。IO 与计算分离。从对象存储拉数据是网络瓶颈重投影是 CPU 瓶颈。最有效的办法是先把需要的变量、按空间窗口裁一遍落到本地临时盘再启动计算任务。这一步多花的时间在计算阶段能十倍地省回来。5.2 存储格式与中间结果管理中间结果用什么格式存长期看影响很大。GeoTIFF通用性最好任何工具都能开适合最终交付。缺点是带压缩之后随机读取性能一般而且一个文件只能存一套网格。NetCDF / Zarr适合多维数据支持分块和压缩能直接对接并行读取。做多时相数据立方体时我几乎只用这个因为它能把“时间×纬度×经度”当成一个整体来切片读某一天的数据不需要打开所有文件。中间结果的命名必须带上下文。我见过最惨的一次事故是管线里有两个版本的重投影参数像元大小不同输出文件名一样跑到一半被覆盖最后一批结果和其他批次的网格差了 3 倍图形上看起来“分辨率突然变粗”排查了两天才定位到问题。从那之后我的输出文件名强制包含数据版本、处理基线、投影代码、像元大小、处理时间戳。长是长了点但省下的调试时间远超命名成本。6. 常见问题与排查技巧实录这一节讲我实际遇到过的问题都是文档里不太会写、但一定会碰上的。6.1 六个高频故障现场故障一读出来的图整体错位海岸线对不上。八成是把逐像素经纬度当成了规则网格的一维坐标轴。表现是图形整体平移几十公里或者出现斜向条纹。排查方法很简单把lat数组的第一行和最后一行打印出来如果同一行内的数值不是单调的、或者行内变化量明显不均匀就说明它是不规则网格必须重投影。故障二云掩膜之后几乎没有有效像元。先检查位号是否正确。快速验证的办法是统计标志位数组的取值分布看是不是只有少数几个值出现。如果所有像元都命中了你排除的某一位那肯定是位定义搞错了。另一个可能是你把“水体”位和“云”位搞反了这种情况在跨基线版本混用数据时特别常见。故障三时间序列出现台阶状跳变。三个可能处理基线混用、定标系数更新、传感器状态变化。用元数据按基线分组画箱线图如果跳变点和分组边界完全重合就确认是版本问题。处理方式要么追溯同基线的完整序列数据量可能不够要么在跳变点前后做相对量分析而不是绝对量分析。故障四批处理跑到一半内存爆掉。除了并发数还要检查是不是有某几景数据尺寸异常。少量产品因为扫描带边缘的填充会比其他景大很多。我的做法是先跑一遍“只读元数据”的预扫描记录每个文件的行列数和文件大小把异常大的挑出来单独限流处理。故障五极区附近图形被严重拉伸。等经纬投影在高纬度会极度变形这是投影本身的特性不是数据问题。解法是换投影比如极射赤面或者兰勃特或者干脆按纬度设置分析范围。做全球产品时很多团队会做两套中低纬用等经纬高纬单独用极区投影最后在分析阶段再拼。故障六多源数据叠加后位置差一个像元。像素坐标的定义方式不同有的产品把坐标定义为像元左上角有的定义为像元中心。差半个像元在 300 米分辨率下就是 150 米的系统性偏移做岸线变化分析时会非常致命。叠加之前先在一个已知地物比如一个形状规整的水库边界上目视核对一次花五分钟省两天。6.2 常见问题速查表现象最可能原因验证方法处理方式图形错位、斜条纹不规则网格未重投影打印经纬度行内变化做重采样后使用有效像元几乎为零标志位定义读错统计标志位取值分布按产品版本映射位号数值全是填充值缺失值未替换检查最大值是否接近 9.99e8读入后立刻置 nan叶绿素系统性偏低二类水体用了标准算法与现场采样比对改自定义指数 本地标定温度整体偏高 1–3 K云污染未剔除干净目视叠加原图掩膜膨胀 只保留高置信度时间序列有台阶处理基线混用按元数据分组画图分组分析或统一基线批处理内存溢出并发过高或单景异常记录峰值内存限流 分块多源叠加偏移像素坐标定义差异已知地物目视核对统一到像元中心定义6.3 几条踩坑经验第一条也是最重要的一条先做一景做到位再谈批处理。我见过太多人拿到数据就直接写循环处理几千景跑到第 800 景发现前面的掩膜逻辑错了全部重来。一景做对的标准不是“程序没报错”是“在一景上把结果和已知事实核对过”——比如找一片公认的清水海域看叶绿素是否落在合理区间找一个已知的高温区看温度场分布是否符合预期。第二条保留中间产物至少保留一个批次。管线调优时你会反复怀疑“到底是输入错了还是算法错了”有中间产物就能一层层往回查。我现在习惯在关键节点各落一份原始读取、掩膜之后、重投影之后、最终结果。占空间但省命。第三条不要迷信产品自带的算法但也不要轻易否定它。我一开始对标准叶绿素算法很不信任因为在研究区表现很差。后来发现是我的研究区本身不适合它换到开阔海域它的表现相当不错而且跨年份一致性好。判断顺序应该是先确认我的应用场景是否符合算法假设再决定要不要自己写。第四条关于验证数据。现场采样的成本很高很多人会拿“别人的论文里的数值”当参考。这个做法在我这里是不通过的——不同时间、不同方法、不同深度的采样值根本不可比。要么自己做同步采样要么用公开的浮标或浮台长时间序列数据做统计比对至少保证时间匹配和量纲一致。最后分享一个我自己常用的小技巧在正式分析前先做一个“裸数据体检”。写一段几十行的脚本把一景数据的行列数、经纬度范围、各变量的有效像元比例、数值分位数、标志位分布全部打印出来。这个动作只花几分钟但能提前发现八成以上的数据质量问题。我现在拿到任何新产品第一件事都是跑这个体检脚本已经成了肌肉记忆。这个内容后续还可以往两个方向扩展一是多源融合和更高分辨率的光学数据做时空融合二是长时序趋势分析配合同一区域的多年代理数据做变化检测——后者对处理基线一致性的要求比我上面讲的还要苛刻动工前务必先把元数据梳理清楚。
1. 先把 Sentinel-3 的四套载荷摸清楚接手这个项目的第一天我做的第一件事不是写代码而是把产品手册打印出来贴在工位上。原因很简单Sentinel-3 不是一颗“拍照片”的卫星它是一个多载荷、多产品级别的对地观测系统同一个轨道数据在命名、结构、坐标系上都有完全不同的样子。你要是按处理普通卫星影像的直觉去理解它大概率会在第三天崩溃。这套系统的核心价值在于业务化的延续性观测单星重复周期约 27 天双星编队后赤道地区重访缩短到一天左右过境时间被锁定在上午光照条件稳定。这两点决定了它和“看到就拍”的商业卫星完全不同——它追求的是同一地点、同一时刻、同一角度、连续多年可比。做过长期时间序列的人都知道“可比”两个字有多值钱它意味着你不需要每年重新标定一次算法阈值。1.1 一颗星上的四种“眼睛”分别负责什么把四套主力载荷拆开看每个都有自己的脾气。OLCI海陆色仪是使用频率最高的一个21 个波段覆盖 400–1020 nm空间分辨率 300 米幅宽 1270 公里推扫式成像。它的波段排布不是随便定的Oa08约 665 nm和 Oa11约 708 nm这对组合是水色反演的主力前者对应叶绿素吸收谷后者对应反射峰Oa11 还夹在荧光基线位置上所以它同时具备一定的藻类荧光探测能力。更关键的是 Oa13约 761 nm和 Oa14约 767 nm这两个氧 A 吸收带波段专门用来做大气校正中的云和气溶胶检测——这是设计时就把“校正”需求前置了不是事后补的。SLSTR海陆表面温度辐射计才是这个系统里最被低估的部分。它最特别的地方是双视角观测同一目标在天底方向和约 55° 倾斜方向各看一次两个路径的大气厚度不同联立起来就能同时反演地表温度和大气水汽。这个设计的直接收益是海表温度精度目标做到 0.3 K 量级且在潮湿热带海域不容易系统性偏移。它有 9 个通道可见/短波红外 500 米热红外 1 公里。做高温目标监测比如火点时短波红外的饱和特性要提前评估不然强火点直接削顶。SRAL合成孔径雷达高度计用 Ku 和 C 两个频段同时工作Ku 频段测距精度可以做到厘米级。它输出的核心量是海面高度、有效波高、后向散射系数。这个载荷和成像类载荷的思维方式完全不一样——它是一条“沿轨剖面”每秒 20 个测量点每个点是一个波形。你要处理的是时间序列不是图像。很多做遥感出身的人第一次碰高度计会不适应因为它没有二维结构只有一条线。MWR微波辐射计只有两个频段约 18.7 和 23.8 GHz任务单一测大气水汽含量给高度计的湿对流层延迟做校正。它不起眼但没有它海面高度的精度会掉一个数量级。提示这四套载荷的时间同步性很好但空间覆盖差别极大。OLCI 幅宽 1270 kmSLSTR 天底幅宽 1420 km、倾斜视角只有 740 kmSRAL 是星下点窄带。做多源融合之前先算一遍覆盖交集别默认“同一景数据覆盖同一区域”。1.2 产品分级与打包格式SEN3 目录到底装了什么下载下来的是一个以.SEN3结尾的目录不是单个文件。这个细节坑过很多人——用文件选择对话框直接选它会失败必须整目录拖走或者用递归传输。目录里通常包含一个 XML 清单文件记录元数据、处理基线、每个测量文件的校验值、一堆.nc格式的测量文件每个波段或每个地球物理量单独一个、坐标文件、时间文件、质量控制标志文件以及仪器辅助数据。这种“一个量一个文件”的设计对按需读取很友好——你只要叶绿素就不用把 21 个波段的反射率全读进内存。产品级别上能直接用的主要是两级1 级是经过辐射定标和几何定位的辐亮度/亮温物理意义明确但没有地球物理含义2 级是经过大气校正后的地球物理量比如叶绿素浓度、海表温度、气溶胶光学厚度。中间还有一个 1A 级主要是原始计数加部分校准日常应用基本不碰。真正需要记住的是后缀。同一颗星同一个载荷不同的两个字母代表完全不同的产品后缀含义分辨率典型用途EFR全分辨率辐亮度300 m自定义算法、高分辨率水体ERR降分辨率辐亮度1.2 km大范围快速浏览、云检测WFR水体全分辨率地球物理量300 m叶绿素、悬浮物、水体透明度WRR水体降分辨率1.2 km全球尺度长时间序列LFR陆地全分辨率300 m陆地反射率、植被指数RBT辐射亮温含双视角500 m / 1 km地表温度、火点WST海表温度1 km海洋温度场LST陆表温度1 km城市热环境、农业FRP火点辐射功率1 km火情监测WAT测高海洋产品沿轨 20 Hz海面高度、波高我见过最典型的误用是拿 WFR 去做内陆浑浊水体的叶绿素反演然后抱怨结果全是零。原因后面会讲这不是数据的问题是算法假设的前提被打破了。1.3 它和哨兵2号、哨兵1号的分工边界经常有人问“我有哨兵2号了还要不要用这个”。答案取决于你要颜色还是要光谱细节。哨兵2号是 10–60 米分辨率、13 个波段、5 天重访擅长陆地细节但它没有专门的氧吸收带波段也没有热红外通道水色反演的波段配置也不是为水体优化的。哨兵3号是 300 米看起来粗但它的波段是为水体和大气专门排布的21 个窄波段能支撑更精细的大气校正和生物光学反演同时还有热红外和海面高度。用一句话概括分工看陆地细节找 2 号算水体参数、测温度和测高找 3 号。真正做海岸带或大型湖泊研究的团队通常是 3 号做长时间序列的“面”2 号做关键时段的“点”两个拼起来用。2. 从需求反推数据选型你要的到底是哪个产品数据选型这一步做错后面所有的代码都是白写。我自己的经验是先回答三个问题我要的物理量是什么、我的目标尺度是多大、我要多久的时间跨度。这三个答案基本就能锁定唯一的载荷和产品级别。2.1 五类典型需求与最小可用产品清单把常见的需求列成表比翻手册快得多需求场景首选载荷最小可用产品需要额外准备开阔海域叶绿素/初级生产力OLCI2 级 WFR 或 WRR无直接读内陆湖泊浑浊度、藻华OLCI1 级 EFR自己算自定义指数 现场验证近岸海表温度SLSTR2 级 WST质量标志筛选陆地地表温度、热岛SLSTR2 级 LST比辐射率假设、地形校正大尺度海面高度与波高SRAL2 级 WAT轨道筛选、重跟踪注意“内陆湖泊”那一行的处理方式这是新手最容易搞错的地方。2 级水体产品默认是按一类水体大洋型设计的它的气溶胶模型和生物光学算法都假设水体光学性质由浮游植物主导。内陆水体里悬浮泥沙和有色溶解有机物占主导这套假设直接失效结果通常表现为浓度被严重低估、异常值大片出现。稳妥的做法是回到 1 级反射率自己做大气校正或者用专门的二类水体处理器再套自定义指数。2.2 文件名就是元数据一次读懂命名规则这个系统的文件名长得离谱但每一个下划线分隔的字段都有意义。举个典型结构S3A_OL_2_WFR____20240115T024000_20240115T024300_20240116T071200_0179_108_123______MAR_O_NT_002.SEN3拆开读S3A平台标识A 星或 B 星。做长时间序列必须区分两颗星的定标差异虽然很小但在趋势分析里必须记录在案。OL载荷代码OL 是色仪SL 是温度辐射计SR 是高度计SY 是协同产品。2产品级别。WFR产品类型。两个时间观测起止时间。这个是筛选数据的关键别拿下载时间当观测时间。第三个时间处理生成时间。0179相对轨道号同一个重访轨道周期内序号固定。108绝对轨道号连续递增。123圈内计数。MAR处理中心标识。O产品类别业务化/测试等。NT时效性等级。002处理基线版本。这一项极其重要不同基线版本的算法、定标系数、标志位定义都可能变混用会造成时间序列跳变。我自己踩过的坑把不同处理基线的产品直接拼成十年序列画出来的趋势线在版本切换那一年有一个明显的台阶看起来像生态突变实际上只是算法升级。注意如果你的分析跨越多年一定要在元数据里保留处理基线字段并在做趋势拟合前检查同一条时间线上是否混了版本。2.3 获取路径的取舍滚动存档、开放数据桶与本地缓存获取渠道大致三类各有适用场景。滚动存档适合“临时取几景”的需求界面友好、按时间区域检索但大批量下载会受单 IP 并发限制且没有断点续传保证。我的做法是只在探索阶段用它确认数据能覆盖研究区之后就切走。云端开放数据桶适合批量和重复计算。它的优势是免下载可以在数据所在的地方直接算把结果几 MB 的 GeoTIFF拉回本地而不是把几个 TB 的原始包搬一遍。做全球尺度分析时这一条几乎决定项目成败。代价是需要写一点对接代码并且要习惯对象存储的“按请求计费、内网访问便宜”的逻辑。本地镜像适合长期固定的研究区。把目标区域三五年的数据一次性拉到本地阵列后续所有实验都是本地 IO迭代速度最快。存储估算要提前做一景全分辨率色仪数据碰巧接近 1 GB一年的双星数据是数千景量级几个 TB 的阵列是起步价。三种方式不是互斥的。我目前的组合是用云端做全球筛查和快速验证用本地镜像做研究区的反复调试用滚动存档补零星缺失的个别时相。3. 读取与预处理netCDF 结构的拆解实操这一章是全文的技术核心。前面都是选型这里开始动手。3.1 xarray 打开 SEN3 目录的几个坑先说最基础的一件事这个目录里的每个.nc文件都是一个独立的 NetCDF 数据集xarray.open_dataset一次只能打开一个。所以你不能直接对.SEN3目录调用它必须精确指定内部文件。import xarray as xr import numpy as np base S3A_OL_2_WFR____20240115T024000_20240115T024300_20240116T071200_0179_108_123______MAR_O_NT_002.SEN3 # 坐标是每个像素都有的必须读 geo xr.open_dataset(f{base}/geo_coordinates.nc) lat geo[latitude].values lon geo[longitude].values # 叶绿素注意变量名和文件名要以实际产品为准先 ls 一遍 chl xr.open_dataset(f{base}/chl_oc4me.nc) print(chl.data_vars) # 先看有什么别猜三个细节容易被忽略第一坐标是逐像素数组不是一维轴。latitude和longitude的维度都是(rows, columns)形如 4000×4865。这意味着它不是一个规则网格而是一个“不规则采样点集合”。后面所有的重投影都要围绕这个事实展开。第二有效值掩膜。产品里通常有一个_FillValue或missing_value属性噪声像元会被填成一个极大值比如 9.99e8 之类。如果不先处理就参与统计你的均值会被拉到天上。稳妥做法是读进来立刻替换为np.nan。第三行列顺序和图像直觉相反。数组的第一个维度是沿轨方向的扫描行第二个是跨轨方向的列。有些产品在存储时做了转置用shape和总像素数交叉验证一下别直接假设。3.2 非规则网格为什么不能直接当 GeoTIFF 用这是新手最容易翻车的地方。你把chl数组和lat/lon数组一起丢给常规栅格工具它会报错或者生成一张完全错位的图因为它没有仿射变换参数。成像类卫星的常规产品带GeoTransform告诉你左上角坐标、像素大小、旋转项所以任意像素的坐标可以用一个线性公式算出来。这个产品没有这个东西它给的是每个像素的真实经纬度而且因为推扫式成像、地球自转、轨道弯曲等因素行与行之间的几何关系并不完全均匀。为什么这么设计因为它保真。推扫式传感器在扫描带两端会有明显的几何拉伸用仿射模型近似会引入系统性位置误差。逐像素给坐标虽然占存储两个 float32 数组4000×4865 就是 300 MB 出头但保证了定位精度。对应的处理逻辑就变了必须做重采样重投影把不规则采样点映射到你指定的规则网格上。这一步不做后面没法做多时相叠加也没法和别的数据源对齐。3.3 三种重投影方案与选型对比实际可用的路子有三条我按踩坑次数从多到少排方案工具优点缺点适用场景图形化处理器重投影桌面处理软件 命令行批处理内置支持好参数少启动慢批处理脚本调试麻烦一次性处理、可视化探索Python 库重采样开源重采样库灵活、可编程、能嵌进管道大数据量要自己管分块生产管线、云端计算先转 VRT 再转栅格栅格工具链生态成熟后续处理方便中间文件大IO 压力大需要与已有栅格工作流对接我目前的默认选择是第二条。原因很实在分块和并行可以自己控制而且能直接接上对象存储的流式读取避免中间文件。图形化处理器那套其实也很好用尤其是命令行模式一条命令就能出结果gpt reproject.xml -Ssourceinput.SEN3 -Poutputout_wgs84.tif -c 4G -q 8其中reproject.xml是操作图定义核心参数就三个目标坐标系、重采样方法、输出像元大小。像元大小怎么定给你一个可用的换算方法目标投影下 1 度纬度约 111 公里要把 300 米分辨率大致保住像元取 0.003 度左右约 333 米。如果做全球长时间序列很多人会降到 0.01 度约 1.1 公里来换存储和速度这时候用降分辨率产品反而更合适不要用全分辨率产品硬降。重采样方法上分类量用最近邻连续量用双线性。叶绿素浓度这种跨数量级的量双线性插值在浓度梯度大的锋面附近会产生虚假的中间值我个人更倾向于最近邻海表温度这种空间连续的场双线性更平滑标志位数组绝对只能用最近邻插值出来的标志位是没有物理意义的。3.4 质量控制标志位与云掩膜所有地球物理量产品都带一个质量标志数组通常是uint32的位掩码bitmask每一位代表一种状况。用法是按位与flags xr.open_dataset(f{base}/wqsf.nc)[WQSF].values # 示意具体位定义必须查对应产品手册的版本说明 CLOUD_BIT 1 LAND_BIT 0 INVALID_BIT 25 cloud_mask (flags (1 CLOUD_BIT)) ! 0 land_mask (flags (1 LAND_BIT)) ! 0 invalid (flags (1 INVALID_BIT)) ! 0 valid ~(cloud_mask | land_mask | invalid)这里必须强调位定义会随处理基线变化。我在不同年份的产品里见过同一含义的位被挪到不同位置的情况。所以你的代码里不应该写死位号而应该写一个“版本到位号”的映射表或者至少在管道里加一步断言检查标志位数组的取值范围和分布是否异常。经验值上一景开阔海域数据里云像元比例通常在 20%–60% 之间如果算出来云占比 99%那八成是位读错了不是你运气差。还有一个常被忽略的点云的“膨胀”。位掩膜给出的是被判为云的像元但云的边缘往往有半透明的薄云和阴影这些像元可能没被标为云数值却已经被污染。稳妥做法是对云掩膜做几次形态学膨胀腐蚀掩膜扩大云的判定范围代价是损失一些有效像元收益是减少异常值。我在做内陆水体时一般膨胀 2–3 个像元开阔海域膨胀 1 个就够。4. 典型应用实操海色、海温与海面高度前面铺垫完了这章讲三个最主流的用法每个都给出可直接复现的思路。4.1 水体叶绿素与浑浊度从标准算法到自定义指数标准产品里给的叶绿素用的是最大波段比值法简而言之就是取几个波段反射率的比值组合再套一个多项式。它的设计前提是一类水体光学性质主要由浮游植物及其降解产物决定其他成分可以近似为常数。在内陆和近岸这个前提碎了。悬浮泥沙会让整体反射率抬升有色溶解有机物会强烈吸收短波蓝光两者都会把比值算歪。表现出来就是明明肉眼可见浓绿色藻华叶绿素反演值却只有个位数。我的处理路径是回到 1 级反射率自己做两个东西第一先做归一化水体指数把水陆分开。用近红外波段的低反射特性做阈值分割这个步骤简单但必须做不然陆地上的植被会被算成高叶绿素。第二用双波段比值指数代替标准算法。比如用约 665 nm 和约 708 nm 构造归一化差值指数def ndci(r665, r708): num r708 - r665 den r708 r665 with np.errstate(divideignore, invalidignore): out num / den return np.where(np.abs(den) 1e-6, np.nan, out)这个指数在富营养化水体里和叶绿素的相关性通常明显好于标准算法代价是它没有物理单位只能做相对分级不能直接报浓度。要用浓度就得有现场采样做回归标定通常准备 20–40 个同步采样点覆盖低中高浓度梯度拟合一条经验曲线。别指望别人的系数能直接搬过来——不同湖泊的悬浮物和溶解有机物组成差别很大系数搬过来误差能到几倍。第三做基线插值指数识别表层藻华。取约 665 nm、708 nm、753 nm 三个波段用前后两个波段在 708 nm 处做线性插值得到一条“无藻华时的基线反射率”实际反射率高出基线的部分就是藻华信号。这个方法的优势是对薄云和观测几何不敏感适合大范围筛查。提示所有基于反射率的指数都必须先完成大气校正。1 级数据是大气顶的辐亮度直接算比值意义不大。如果暂时没有好的大气校正方案退而求其次的做法是用两个邻近波段的比值因为气溶胶影响在两个邻近波段上高度相关比值能部分抵消——但这是妥协不是正解。4.2 海表温度双视角的价值在哪海表温度产品的核心难点从来不是测温度是剥离大气。水汽在热红外波段的吸收很强而海面水汽含量变化又大单视角观测无法同时确定“地表温度”和“大气水汽”两个未知数。双视角就是为这个设计的。同一点被看两次两次穿过的大气路径长度不同吸收量不同两个方程解两个未知数。理论上的收益是水汽含量高的热带海域单视角产品的偏差能到 1 K 以上双视角可以压到 0.4 K 以内。实际使用时的几个注意事项云检测宁严勿宽。热红外对亚像元云极其敏感一个 10% 覆盖的薄云就能让温度偏低好几度。产品的质量标志里通常有置信度分级我一般只保留最高等级。近岸要打问号。1 公里像元在近岸会混入陆地信号尤其是潮间带和浅水区。这不是算法能解决的只能靠掩膜剔除或者用更高分辨率数据补。白天和夜间产品要分开看。白天的太阳反射会引入偏差很多应用场景比如日较差分析必须按地方时筛选观测。验证不要用陆地气象站。这是一个很基础的错误但经常出现拿岸边气象站的气温去验证海表温度两者物理意义完全不同一个在陆地 2 米高度一个是水面皮层温度。合理的做法是用浮标数据且要匹配观测时间窗口一般取 ±2 小时内。4.3 测高波形重跟踪与轨道级时间序列高度计的处理思路和成像完全不同我单独拿出来讲因为这块的坑最隐蔽。基本几何关系是卫星到海面的距离由雷达往返时间测得卫星轨道高度已知两者相减得到海面相对参考椭球的高度。这个高度里混合了大量信号大地水准面起伏、潮汐、海况偏差、对流层和电离层延迟。2 级产品里这些校正项基本都提供好了你要做的是选择性地应用并检查每一项的合理性。主要校正项我自己列了一个核对顺序干对流层由大气压算、湿对流层由微波辐射计算这是那个双频辐射计的作用、电离层双频雷达可以自校正、逆气压效应、固体潮、海洋潮、极潮、海况偏差。每一项都是厘米到分米量级海面高度异常信号本身也就十几厘米所以少一项误差就等于把信号淹了。重跟踪是绕不过去的一步。在开阔大洋标准的波形拟合算法足够好但在近岸、海冰边缘、内陆水域回波波形会被陆地或冰面污染标准算法给出的距离偏差可达几十厘米。解决办法是用专门的重跟踪算法重新估计波形前沿。目前的实践是用现成的重跟踪器参数上用窄的跟踪窗口和专门针对近岸的拟合模型。时间序列有一个关键技巧按轨道号筛选不要按经纬度框筛选。高度计的采样是沿轨稀疏的如果你用一个地理范围去裁不同重访周期裁出来的采样点位置会漂移导致时间序列里混进了空间差异。正确做法是固定相对轨道号这样每次过境走的是几乎完全相同的星下点轨迹空间位置一致差异才真的是时间差异。4.4 一个能跑的最小流程把上面的东西串起来给一段可以直接改参数用的骨架代码。我自己写管线时习惯用它起步验证数据读取和掩膜是否正常再往里塞业务逻辑。import numpy as np import xarray as xr import rioxarray from pyproj import CRS def load_olci_l2(base_dir, var_filechl_oc4me.nc, var_nameCHL_OC4ME): 读取 OLCI 2 级产品返回经纬度和目标变量缺失值统一成 nan geo xr.open_dataset(f{base_dir}/geo_coordinates.nc) lat geo[latitude].values lon geo[longitude].values ds xr.open_dataset(f{base_dir}/{var_file}) raw ds[var_name].values # 关键一步把填充值换成 nan fill ds[var_name].attrs.get(_FillValue, 9.99e8) data np.where(raw fill * 0.9, np.nan, raw) return lat, lon, data def make_valid_mask(base_dir, exclude_bits(0, 1, 25)): 按位排除陆地、云、无效像元bit 定义需按产品版本核对 flags xr.open_dataset(f{base_dir}/wqsf.nc)[WQSF].values bad np.zeros_like(flags, dtypebool) for b in exclude_bits: bad | (flags (1 b)) ! 0 return ~bad def build_point_table(lat, lon, data, valid, step4): 降采样抽稀成点位表交给后续 GIS 或数据库处理 sub (slice(None, None, step), slice(None, None, step)) lat, lon, data, valid lat[sub], lon[sub], data[sub], valid[sub] ok valid np.isfinite(data) (lat -89.9) (lat 89.9) return lon[ok], lat[ok], data[ok]step4这个抽稀参数不是随手定的。300 米分辨率的全分辨率产品一景有效像元可能有几百万个直接做点密度图会把内存吃穿。降采样 4 倍后像元间距约 1.2 公里对做区域态势图完全够用内存压力降一个数量级。如果你要做逐像元的精细统计那就别抽稀改成按块流式写入文件。5. 批处理工程化与性能优化单景处理跑通之后真正的难题才开始你有几百上千景数据怎么让它在可接受的时间里跑完。5.1 内存、分块与并行一景全分辨率色仪的反射率21 个波段全部读进来按 4000×4865 的尺寸算是 800 MB 量级加上坐标、标志位、中间结果内存峰值轻松破 3 GB。看着不大但如果你用多进程跑 8 个并发就是 25 GB普通工作站直接被打爆。三个可用的手段分块读取。用支持懒加载的数组结构把数据切成 1024×1024 的块处理完一块写一块内存占用就固定在一个块的大小上。这个改动通常能把内存峰值降到几百兆。限制并发数而不是死磕单进程速度。我的经验值是按“单景峰值内存 × 并发数 ≤ 物理内存的 70%”倒推并发数剩下的 30% 留给系统缓存和垃圾回收。宁可少开两个进程也不要让系统开始换页——一旦开始换页整体速度会慢十倍。IO 与计算分离。从对象存储拉数据是网络瓶颈重投影是 CPU 瓶颈。最有效的办法是先把需要的变量、按空间窗口裁一遍落到本地临时盘再启动计算任务。这一步多花的时间在计算阶段能十倍地省回来。5.2 存储格式与中间结果管理中间结果用什么格式存长期看影响很大。GeoTIFF通用性最好任何工具都能开适合最终交付。缺点是带压缩之后随机读取性能一般而且一个文件只能存一套网格。NetCDF / Zarr适合多维数据支持分块和压缩能直接对接并行读取。做多时相数据立方体时我几乎只用这个因为它能把“时间×纬度×经度”当成一个整体来切片读某一天的数据不需要打开所有文件。中间结果的命名必须带上下文。我见过最惨的一次事故是管线里有两个版本的重投影参数像元大小不同输出文件名一样跑到一半被覆盖最后一批结果和其他批次的网格差了 3 倍图形上看起来“分辨率突然变粗”排查了两天才定位到问题。从那之后我的输出文件名强制包含数据版本、处理基线、投影代码、像元大小、处理时间戳。长是长了点但省下的调试时间远超命名成本。6. 常见问题与排查技巧实录这一节讲我实际遇到过的问题都是文档里不太会写、但一定会碰上的。6.1 六个高频故障现场故障一读出来的图整体错位海岸线对不上。八成是把逐像素经纬度当成了规则网格的一维坐标轴。表现是图形整体平移几十公里或者出现斜向条纹。排查方法很简单把lat数组的第一行和最后一行打印出来如果同一行内的数值不是单调的、或者行内变化量明显不均匀就说明它是不规则网格必须重投影。故障二云掩膜之后几乎没有有效像元。先检查位号是否正确。快速验证的办法是统计标志位数组的取值分布看是不是只有少数几个值出现。如果所有像元都命中了你排除的某一位那肯定是位定义搞错了。另一个可能是你把“水体”位和“云”位搞反了这种情况在跨基线版本混用数据时特别常见。故障三时间序列出现台阶状跳变。三个可能处理基线混用、定标系数更新、传感器状态变化。用元数据按基线分组画箱线图如果跳变点和分组边界完全重合就确认是版本问题。处理方式要么追溯同基线的完整序列数据量可能不够要么在跳变点前后做相对量分析而不是绝对量分析。故障四批处理跑到一半内存爆掉。除了并发数还要检查是不是有某几景数据尺寸异常。少量产品因为扫描带边缘的填充会比其他景大很多。我的做法是先跑一遍“只读元数据”的预扫描记录每个文件的行列数和文件大小把异常大的挑出来单独限流处理。故障五极区附近图形被严重拉伸。等经纬投影在高纬度会极度变形这是投影本身的特性不是数据问题。解法是换投影比如极射赤面或者兰勃特或者干脆按纬度设置分析范围。做全球产品时很多团队会做两套中低纬用等经纬高纬单独用极区投影最后在分析阶段再拼。故障六多源数据叠加后位置差一个像元。像素坐标的定义方式不同有的产品把坐标定义为像元左上角有的定义为像元中心。差半个像元在 300 米分辨率下就是 150 米的系统性偏移做岸线变化分析时会非常致命。叠加之前先在一个已知地物比如一个形状规整的水库边界上目视核对一次花五分钟省两天。6.2 常见问题速查表现象最可能原因验证方法处理方式图形错位、斜条纹不规则网格未重投影打印经纬度行内变化做重采样后使用有效像元几乎为零标志位定义读错统计标志位取值分布按产品版本映射位号数值全是填充值缺失值未替换检查最大值是否接近 9.99e8读入后立刻置 nan叶绿素系统性偏低二类水体用了标准算法与现场采样比对改自定义指数 本地标定温度整体偏高 1–3 K云污染未剔除干净目视叠加原图掩膜膨胀 只保留高置信度时间序列有台阶处理基线混用按元数据分组画图分组分析或统一基线批处理内存溢出并发过高或单景异常记录峰值内存限流 分块多源叠加偏移像素坐标定义差异已知地物目视核对统一到像元中心定义6.3 几条踩坑经验第一条也是最重要的一条先做一景做到位再谈批处理。我见过太多人拿到数据就直接写循环处理几千景跑到第 800 景发现前面的掩膜逻辑错了全部重来。一景做对的标准不是“程序没报错”是“在一景上把结果和已知事实核对过”——比如找一片公认的清水海域看叶绿素是否落在合理区间找一个已知的高温区看温度场分布是否符合预期。第二条保留中间产物至少保留一个批次。管线调优时你会反复怀疑“到底是输入错了还是算法错了”有中间产物就能一层层往回查。我现在习惯在关键节点各落一份原始读取、掩膜之后、重投影之后、最终结果。占空间但省命。第三条不要迷信产品自带的算法但也不要轻易否定它。我一开始对标准叶绿素算法很不信任因为在研究区表现很差。后来发现是我的研究区本身不适合它换到开阔海域它的表现相当不错而且跨年份一致性好。判断顺序应该是先确认我的应用场景是否符合算法假设再决定要不要自己写。第四条关于验证数据。现场采样的成本很高很多人会拿“别人的论文里的数值”当参考。这个做法在我这里是不通过的——不同时间、不同方法、不同深度的采样值根本不可比。要么自己做同步采样要么用公开的浮标或浮台长时间序列数据做统计比对至少保证时间匹配和量纲一致。最后分享一个我自己常用的小技巧在正式分析前先做一个“裸数据体检”。写一段几十行的脚本把一景数据的行列数、经纬度范围、各变量的有效像元比例、数值分位数、标志位分布全部打印出来。这个动作只花几分钟但能提前发现八成以上的数据质量问题。我现在拿到任何新产品第一件事都是跑这个体检脚本已经成了肌肉记忆。这个内容后续还可以往两个方向扩展一是多源融合和更高分辨率的光学数据做时空融合二是长时序趋势分析配合同一区域的多年代理数据做变化检测——后者对处理基线一致性的要求比我上面讲的还要苛刻动工前务必先把元数据梳理清楚。
