第一次拿到这份海冰浓度数据的时候我盯着文件名里那串“siconc”看了半天——Nimbus-7 SMMR、DMSP SSM/I-SSMIS、V002这些词拆开都认识但合在一起就是一份横跨四十多年的被动微波海冰浓度序列。当时我正要做北极海冰面积的长时序趋势分析最要紧的就是搞清楚这份数据是怎么把原始亮温变成“百分之多少海冰覆盖”的以及不同卫星传感器之间到底怎么拼接才不出错。这篇文章就把我实际用这份V002产品时的理解和踩过的坑一次性写清楚适合刚开始接触极区遥感、或者已经拿到数据但不太放心直接开算的研究生和工程师参考。1. 为什么海冰浓度要靠被动微波一场关于“看不见”的遥感接力1.1 为什么不是可见光极夜和云层决定了技术路线先说一个很多人刚入坑时会困惑的问题海冰看着清清楚楚的为什么不用光学卫星拍照片然后数像素原因其实特别直白——极地一年里有半年是极夜太阳都不出来可见光遥感直接抓瞎。就算在极昼的夏季北极地区也是出了名的多低云MODIS那一类光学影像经常十天半月找不到一块干净的晴空。对做气候监测的人来说等不起这个天气窗口。需要有东西能穿透云层、不依赖太阳光照而且白天黑夜都能用。这就是被动微波遥感登场的根本原因它接收的是地表和大气自身辐射的微波不是反射太阳光。1.2 被动微波的物理依据发射率差和极化差微波辐射计测量的是亮温brightness temperature单位是开尔文。对一个黑体而言亮温就是它的物理温度对真实地表来说亮温等于物理温度乘以发射率。海水和纯冰在物理温度上差别并不大海面零下1.8°C左右冰面可能零下十几度真正拉开距离的是发射率。在10到90GHz这个频段开阔海面的微波发射率大约只有0.4到0.5所以海水亮温通常在180到200K之间。海冰发射率则能到0.9以上亮温直接跳到240到250K。这一下就拉开了四五十五十度的亮温差比物理温差的贡献大得多。也正因为发射率差海冰浓度反演成了被动微波应用里最稳的一个方向。这里还有个关键细节是极化。辐射计一般有两个极化通道垂直极化和水平极化同时观测。光滑海面在微波波段就像一个镜面它对垂直极化和水平极化的发射差异非常大而海冰表面粗糙、内部还有大量气泡和盐晶结构辐射经过体散射后极化差异被抹平。所以海冰“看得见”除了靠绝对的亮温高还靠极化差小。这个特性直接为后面算法里的极化比提供了物理基础。1.3 线性混合假设浓度是怎么被“算”出来的海冰浓度Sea Ice ConcentrationSIC的定义很简单一个像元里海冰覆盖面积占像元总面积的比例0到100%。被动微波遥感器无法识别像元内部细节它测到的亮温是一个“混合”信号。最简单最朴素的做法就是线性混合假设像元亮温 冰区亮温 × 浓度 水区亮温 ×1 - 浓度然后反解浓度。但实际地表哪有这么听话。大气水汽、云中液态水、海面风浪、海冰上的积雪、融池都会改变亮温。所以从1970年代到现在几十种反演算法本质上都是围绕一个问题打转如何在亮温里把海冰信号和环境噪声分开。NSIDC这类海冰浓度产品里最经典的两套解法一个是NASA Team算法一个是Bootstrap算法都要在消除天气影响这件事上做文章。到了下一节你就能理解为什么需要三台不同类型的传感器组合起来才能拼出一条从1978年延续到现在的序列。2. 横跨45年的传感器家谱SMMR、SSM/I、SSMIS 到底差在哪2.1 三个传感器的技术参数对比这份产品名里的三个传感器是三个时代的三台主力。把它们的关键参数放在一起看差异就很清楚了。传感器搭载卫星工作时段海冰产品常用主要频段GHz极化入射角SMMRNimbus-71978.10 - 1987.086.6, 10.7, 18, 21, 37垂直/水平V/H约50.3°SSM/IDMSP F8-F151987.07 - 200919.35, 22.235, 37, 85.519/37/85为V/H22仅V约53.1°SSMISDMSP F16-F192003 - 至今19.35, 22.235, 37, 91.655等19/37/91为V/H22仅V约53.1°注意一个细节SMMR的频段是18GHz而SSM/I和SSMIS是19.35GHz差了1.35GHz。频率偏移直接导致在相同地表条件下测到的亮温有区别——频率越高大气吸收和发射越强同一种地表在不同频率下亮温表现也会不同。所以跨传感器拼接时不能直接把这几个传感器的相同“号位”通道当成同一个通道来用。2.2 拼接一个连续时间序列的动机你可能会问直接用其中一台传感器不就行了为什么要拼因为气候研究需要的是几十年尺度的趋势没有任何一台卫星能连续工作四十五年。SMMR从1978年干到1987年SSM/I接棒到2000年代SSMIS从2003年逐步接班并一直运行到现在。只有把这三位的数据拼接起来才能得到1978年至今的完整记录。这份产品本质是三个传感器接力跑马拉松的结果。另外SMMR还有个不太方便的特点由于搭载卫星的功耗限制它并非每天都开机扫极区早期是隔一天甚至隔两天的覆盖。所以在1978到1987年这段数据时间分辨率明显低于后来的SSM/I。做趋势分析时最好用月平均产品或者对日数据先聚合避免因为早期数据稀疏造成假信号。2.3 传感器之间最麻烦的差异过境时间与扫描几何传感器拼接最麻烦的不是频段数而是过境时间和扫描几何不一致。DMSP卫星是极轨气象卫星每颗星都有固定的地方时过境时间。F8大概在早上六点到晚上六点过境F11在下午五点半左右F17则在下午五点到五点半之间。海冰本身移动变化不算快但海冰表面的温度日变化、雪面状态、甚至海面风浪都会跟着地方时变。当传感器换了你看到的海冰浓度除了真实变化还混入了观测时刻不同的系统偏差。扫描几何也同样关键。这三台传感器都是圆锥扫描但有约3度的入射角差异。入射角不同地表发射率和大气路径长度都不同。SMMR更严重的是不同频段的足印分辨率相差很大37GHz只有大约27km分辨率而10.7GHz的足印将近100km。在做多频段联合反演时高分频率和低分频率数据必须匹配到同一个尺度这本身就是产品预处理里非常容易出错的一环。2.4 交叉校准是怎么做的为了把这些差异压到最小产品制作方在两代传感器重叠期做了交叉校准。所谓重叠期就是老传感器和新传感器同时还在工作的那段时间最典型的是1987年7月到8月的SMMR与SSM/I重叠以及2003年之后SSM/I与SSMIS的重叠。校准思路并不复杂把同一时间段、同一地点的亮温或海冰浓度拿来做回归拟合出经验偏移量再用这个偏移量把前一代数据统一到后一代传感器的参考框架下。但我必须提醒交叉校准不完美。1987年SMMR到SSM/I的交接至今仍然是海冰时间序列里一个公认的“台阶”位置很多研究文章都报告过海冰面积在那一两个星期出现台阶式跳变。V002版本花了大力气重新做这一段的修正但这也意味着它和V1产品在交接期的表现明显不同——如果你把两个版本的1987年到1990年数据拿来接在一起比较趋势大概率会发现人为的断点。3. 从亮温到百分比NASA Team 与 Bootstrap 的算法原理拆解3.1 为什么算法离不开37GHz和19GHz这两个频段看产品变量名可能只会看到一个简单的siconc但它的背后是一整个频段选择逻辑。反演海冰浓度最常用的是19GHz和37GHz外加22GHz用来做大气校正。为什么偏偏是这两个频率原因有两点。第一频率不能太低。低于10GHz的频段虽然大气影响小但分辨率太差一个像元几百平方公里海冰边缘的细节全被抹平而且低频辐射对海面粗糙度的响应太强海况一变亮温就乱跳。第二频率不能太高。到了85GHz以上大气水汽和云的影响变得非常重虽然分辨率好但亮温里掺了太多天气噪声。19GHz和37GHz是两个折中比较好的频段37GHz能提供较好的空间分辨率19GHz对大气和液态水的敏感性相对低、又能和37GHz形成差异信号。所以几乎所有经典算法都围绕这两个频段打转。3.2 NASA Team极化比、梯度比和那个经典公式NASA Team算法后面就叫NT算法是1984年由Cavalieri等人提出的也是这份产品里最核心的算法。它的输入只有三个量37GHz垂直和水平极化亮温T37V、T37H、19GHz垂直极化亮温T19V。先定义两个参数极化比 PR (T37V - T37H) / (T37V T37H)梯度比 GR (T37V - T19V) / (T37V T19V)极化比反映的是地表极化特性的强弱。开阔水面的PR很高差不多0.2到0.3因为水的水平极化发射率显著低于垂直极化海冰则因为体散射和粗糙面PR掉到0.05以下。所以PR直接携带了“这里是水还是冰”的信息。梯度比反映的是两个频率之间的亮温差异。这里有个巧妙之处大气中的云滴和水汽会强烈影响37GHz而19GHz相对更“干净”。如果GR异常高说明37GHz被大气“抬高”了数据很可能受了天气污染。因此GR不仅参与浓度反演还被当作天气滤波器来用典型阈值在0.05上下超过这个值就认定亮温被污染浓度结果打入无效或标记低质量。经典NT公式把PR和GR放进一个统计回归表达式里总浓度 C (a0 a1·PR a2·GR a3·PR·GR) / (b0 b1·PR b2·GR b3·PR·GR)其中a和b是一组经验系数来自辐射传输模型模拟和实测数据的标定。Cavalieri等人在1984年文章里的经典系数大概是a01014.6a11869.0a2-1560.3a3-343.7b0887.5b11843.6b2-1366.6b3-296.6。你不需要手算产品里已经替你算好了但理解这个结构还是有用的它本质上是用PR和GR两个参数在二维空间里画了一个曲面把冰水混合比例映射到0到100。3.3 Bootstrap二维特征空间里的线性内插NT算法是主算法但产品还会提供另一个算法做交叉验证常用的是Bootstrap算法。Bootstrap的思路跟NT完全不同它用19V和37V两个通道的亮温构成一个二维特征空间先标定出这个空间里的“纯水点”和“纯冰点”然后把每个像元的亮温点放到这两个参考点之间按距离做线性内插得到浓度。具体操作里开阔水域因为海温和盐度变化会在特征空间里聚成一个团海冰由于类型和雪厚不同会形成从一年冰到多年冰的一条延伸带。部分冰覆盖的像元应该落在水和冰这两个端元之间的连线附近位置越靠近冰端浓度越高。Bootstrap的优点是它对大气水汽的敏感度相对低因为它在处理时会把天气影响较大的样本先剔除缺点是它更依赖端元标定的准确性如果某个区域冰结构特殊就容易被低估或高估。一个在实测里非常明显的区别在夏季融冰期融池水会把海冰表面发射率压下去NT算法有时会把这种像元当成“开放水”导致浓度偏低Bootstrap在特定区域表现相对稳定。但反过来在多年冰比例很高的区域Bootstrap容易把多年冰当成更高浓度NT则可能混入积雪的影响。所以两个算法各有所长交叉使用更稳妥。3.4 两个算法结果应该信谁碰到这个问题我一般这样处理冬季、无雪或少雪时期NT算法和Bootstrap算法差距不大取其一即可如果要和国际上其他研究直接对比优先用NT算法结果因为文献里引用的“SSM/I海冰浓度”默认就是NT结果。夏季和有融池的季节两个算法差异会拉大这时候别急着选一个就完事至少要看看它们差多少。如果差异超过10%大概率是天气污染或融池影响对应的像元应该被特殊对待。另外要注意这份产品里的变量实际是用亮温直接算出来的单位是百分比。有些衍生产品的单位是0到1的小数做面积统计之前一定要先确认单位不然面积会差一百倍。4. V002版本升级的核心网格、校准和岸线掩膜做了哪些事4.1 说了很久的EASE-Grid 2.0终于换上了V002最直观的变化是投影网格从老一代的EASE-Grid 1.0换到了EASE-Grid 2.0。听起来只是换个名字实际影响不小。EASE-Grid是一套等积投影网格最大特点是在网格面积上保持准确这对海冰面积和趋势统计非常重要。2.0版本重新定义了投影椭球体、标准纬线等参数网格点的实际地理位置和老版本相比有几公里的平移。不要小看这几公里。在北极海冰边缘像元位置只要偏上几公里就可能把格陵兰东岸的冰间水道归错类在做历史序列对比时如果一个研究用了V1、另一个用了V2直接比较像元值是危险的。幸好EASE-Grid 2.0的投影参数是公开的各家GIS软件也普遍支持建议把数据处理流程统一到V2的网格上不要一张图混两种网格。V002的南北半球网格同样都是25km分辨率。北半球网格维度一般是304×448南半球是316×332经纬度已经是预计算好的二维坐标数组方便直接绘图和统计。拿到手里你会发现数据里除了siconc还有对应的latitude、longitude、x、y坐标处理起来比老的极坐标文件友好得多。4.2 天气滤波和海岸线掩膜是两项主要实质改动版本说明里不太会大张旗鼓宣传的是天气滤波阈值和海岸线掩膜的重新处理而这两项恰恰是我实测中感受最明显的改动。先说天气滤波。老版本对天气污染的剔除比较保守在北大西洋和北太平洋一侧经常能看到云中液态水造成的“虚假海冰”——一片开阔水域显示成带状或斑状的高浓度冰。V002把天气滤波的判定逻辑重新调了一遍GR阈值和37V阈值都有微调假冰现象明显减少。代价是少数真实的海冰边缘像元会被误删这种取舍在算法里很常见没有完美的滤波器。再说岸线掩膜。老版本最大的问题在海岸附近陆地像元的高亮温会污染近岸海洋像元导致沿海海冰浓度被低估甚至出现一条沿海“无冰带”现实中却是冰情严重的快速冰区。V002用更精细的海岸线数据重新做了掩膜同时增加了对陆地邻近像元的处理。我在波弗特海和鄂霍次克海沿岸做验证时近岸像元的数据质量明显提升和SAR影像的吻合度也更好。4.3 早期SSMIS的传感器校正被重新处理前面说过DMSP F16卫星上SSMIS早期数据有个著名的天线温度漂移问题。F16在2003年发射后的前几年高纬度亮温出现了随时间缓慢变化的系统误差来源是冷空反射器受到非均匀温度影响。这个问题直接污染了原来直接用原始亮温反演的海冰浓度尤其在2005到2009年这段有些地区的海冰面积异常偏高或偏低和邻近传感器对不上。V002对这一段数据专门做了重新校准。实际使用中我建议在2003到2009年之间要多长个心眼即使经过了校正这段的F16数据仍然是全序列里不确定性最高的区间。如果研究窗口正好落在这里最好用F15或F17的SSM/I/SSMIS数据做一次对比验证。数据文件里一般会保留传感器来源标记做质量筛选时别把这些元数据丢掉。4.4 如何确认你拿到的是V002下载时容易在文件名上栽跟头。NSIDC这类数据中心的命名规则里通常包含版本号有的叫v002、v2.0、v02r00有的是在XML元数据里标注。我建议下载完文件之后不要只看文件名直接打开NetCDF或HDF的全局属性找版本号字段确认。另外还要检查文件的最后修改时间——很多数据中心会把新版数据重新处理后在同一个目录下静默更新文件名一样、版本号一样但内容变了。做科研要留好数据溯源记录别嫌麻烦。5. 打开数据前的5个关键参数网格、单位、坐标和读图示例5.1 数据集的时空覆盖和分辨率这份V002产品的总体时间范围是从1978年10月26日Nimbus-7 SMMR开始工作起一直延伸到最近后续会滚动更新。产品里有日数据和月数据两个层次。日数据直接反映每天的冰况但受天气噪声影响大月数据是逐像元平均信噪比更好做趋势分析时强烈推荐先看月数据。空间分辨率是25km的EASE-Grid像元。25km听起来比MODIS那种几百米的细节差了很远但对海冰面积估算来说25km网格完全够用——北极海冰总面积的量级是几百万平方公里25km像元的面积误差摊到总面积上是可以接受的。但如果你要研究狭窄的冰间水道、峡湾或小岛周围的冰况25km是远远不够的应该换冰芯SAR或者更高分辨率的微波产品。5.2 单位、缺测值和编码规则打开变量后要立刻检查三个属性单位、有效值范围、缺测值。siconc的单位是百分比%范围是0到100。低于0或高于100的都是无效值或经过处理后的极端值统计分析前必须做掩膜剔除。还有一个容易踩的坑有些版本里用特殊数值表示“陆地”“外部区域”“极洞”。比如陆地的掩膜值可能是255byte类型或NaN。如果你直接用255去算平均浓度结果会被拉得极其离谱。正确做法是先读取数据质量掩膜或陆海掩膜单独处理。5.3 Python读取与快速出图下面是一段我经常用的最小示例代码直接读取V002的NetCDF并快速画一张北半球海冰浓度图。前提是已经装好xarray、matplotlib和cartopy。import xarray as xr import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature import numpy as np # 文件名按实际文件替换 ds xr.open_dataset(nt_20200301_f17_north_v2.0.nc) # V002的siconc通常是三维[time, y, x]先squeeze掉单例维度 sic ds[siconc].squeeze() # 从数据文件里取经纬度坐标EASE-Grid 2.0通常提供二维lon/lat lons ds[longitude].values lats ds[latitude].values # 北极投影示例 proj ccrs.LambertAzimuthalEqualArea(central_latitude90, central_longitude0) ax plt.axes(projectionproj) # 用pcolormesh画像元色块注意transform参数指明原始坐标是经纬度 mesh ax.pcolormesh( lons, lats, sic, transformccrs.PlateCarree(), cmapice, vmin0, vmax100 ) ax.add_feature(cfeature.COASTLINE, linewidth0.5) ax.set_extent([-180, 180, 50, 90], crsccrs.PlateCarree()) plt.colorbar(mesh, labelSea Ice Concentration (%)) plt.title(SMMR-SSM/I-SSMIS SIC V002 (2020-03-01)) plt.show()如果你的数据文件里没有直接的longitude和latitude变量就需要根据EASE-Grid 2.0投影参数用pyproj现场计算。这个方法更稳因为不同机构对网格坐标的命名不完全一致但投影参数是标准化的。用pyproj的Proj类定义好投影后把x、y网格坐标转成经纬度即可。5.4 关于投影转换的两个提醒第一个提醒25km EASE-Grid像元面积不是处处固定。虽然EASE-Grid是等积投影但那是针对投影面上的像元而言换算回经纬度坐标后一个“像元”在不同纬度的实际面积略有差异。北极附近由于投影极点处聚集面积误差最明显。统计总面积时要么用投影坐标下像元面积乘以像元数要么用官方提供的像元面积变量不要简单用等经纬度格点面积公式近似。第二个提醒国际日期变更线附近要小心。EASE-Grid网格的经度范围通常在-180到180之间在绘图和做时间序列叠加时日界线附近的像元可能因为经度跳变而错位。解决办法是用投影坐标直接处理或者在插值前把经度转换到0到360的统一框架不要在PlateCarree坐标下盲目做空间插值。6. 实测最容易翻车的四个地方北极洞、天气污染、传感器交接和夏季融池6.1 北极洞那个永远没有数据的极点第一次画北极海冰图时你会发现北极点周围有一个圆形的数据空洞格陵兰以北的北冰洋中心区域在图上就像一个洞。这不是数据缺失而是传感器扫描幅宽根本覆盖不到极点。SMMR、SSM/I和SSMIS都是倾斜圆锥扫描轨道倾角大约98.8度卫星轨迹最北只能到北纬84到87度附近极点周围就永远处于扫描带之外。这个洞在SMMR时代更大SSM/I时代稍微小一点。处理北极洞有几种常见策略。有些机构会直接用最靠近极点的一圈观测值把洞填平方便算海冰总面积更严格的做法是像NSIDC的某些CDR产品那样先用传感器数据和气候学插值把洞补上并生成一个“极洞掩膜”说明哪些像元是插值出来的。用之前一定要先搞清楚产品里填洞没有如果没填你要么放弃计算极点附近面积要么用插值填补并如实说明方法。做海冰面积趋势时如果前后期产品的填洞策略不一样这个洞的大小变化本身就会产生虚假趋势。6.2 天气污染冬天看GR阈值夏天看融池被动微波海冰浓度最大的敌人是云中液态水。低温的云滴在37GHz频段有很强的发射会把“海面”亮温抬到接近冰的亮温水平于是反演结果把开阔水域误判成海冰。这在气旋频繁经过的冬季阿留申群岛、挪威海、拉布拉多海一带特别常见。判断方法很简单打开该区域的GR变量如果GR显著超过天气滤波阈值浓度数据基本不可信。产品里往往已经做了这一步掩膜但漏网的仍不少。夏季翻车的机制完全不同。海冰表面开始融化后表面会形成融池融池水的发射率接近海水于是本应是100%浓度的像元亮温大幅下降反演出50%、60%。这就导致夏季海冰浓度系统性偏低特别是在边缘海。做夏季趋势分析时别把南半球冬季北半球夏季的数据直接拿来和北半球冬季一起用那会混入严重的季节性算法偏差。一个实用办法是只用高质量标记的像元或者用多种算法互相核对浓度差距太大的像元宁可剔除也不要硬用。6.3 传感器交接期的时间序列“台阶”前面提到1987年的SMMR到SSM/I交接是海冰面积序列里最知名的台阶。实际作图你会看到在1987年7月前后海冰面积曲线往往出现一个几万到几十万平方公里的跳变。这个跳变不是天气导致的而是传感器设计不同留下的痕迹。2003年SSM/I到SSMIS的交接期也有类似现象程度较轻。应对办法是第一做趋势分析时明确写出自己用了哪个版本、跨了哪些传感器不要假装数据是完全均一的第二如果研究不依赖1978到1987年的早期数据可以直接从1988年开始算趋势避开最严重的交接期第三如果一定要用全长数据可以引入传感器作为协变量做分段回归或者参考专门做海冰CDR气候数据记录的产品那些产品在交接处理上花过更多小心思。6.4 时间序列分析前我建议先做的事最后分享一个我自己的处理习惯。拿到V002数据后我不会急着算全北极平均浓度而是先按下面几步走用月平均数据不用日数据先把每天的天气噪声平滑掉。把siconc小于15%的像元单独看待。微波算法在低浓度区间的误差急剧增大15%以下像元只能定性使用不宜参与定量面积统计。对每个像元生成“有效数据天数”统计有些像元在某些月份常年被天气污染应该在年际比较时剔除。保存一份处理后的掩膜文件和中间结果万一以后换了数据版本还能溯源排查。这样一套流程走下来再去做北极海冰面积趋势结果会稳定很多。我自己在实际操作中感受最深的一点是再好的反演产品也只能告诉你“亮温信号大致对应多少冰”它无法替代你对某个区域海冰物理过程的理解。拿到V002数据后花点时间把三个传感器、两个算法、北极洞和天气滤波这些底层细节弄清楚远比直接跑一个回归有意义得多。