3. 温盐剖面数据预处理实战
拿到一条CTD剖面,第一件事不是画T-S图,而是判断这条剖面“能不能用”。原始温盐剖面几乎不可能直接进入分析环节:传感器在下降过程中会出现短暂失锁、压力通道存在零点偏移、盐度在浅层因泵未启动而出现异常、个别采样点因电磁干扰形成尖峰。本节以一条典型的温盐深剖面(示例数据,深度0~1500 m,采样间隔1 dbar,含少量人为注入的缺失与尖峰)为对象,完整走一遍清洗链路,并说明每一步的判断依据。
先建立参照约定与数据体检
本节统一约定:深度向下为正,压力 p 以 dbar 为单位,1 dbar 约等于 1 m;温度 T 用 ITS-90 温标,单位 ℃;盐度 S 采用实用盐度 PSU。所有校正量都按“观测值 + 校正量 = 校正后值”的方向书写,避免正负号混淆。
清洗前先做体检,重点看三件事:一是缺失分布,是整段缺失还是零星缺测;二是压力序列是否单调,下降剖面中压力应基本单调递增,若出现回退说明存在船体起伏或采样异常;三是温盐的垂向梯度是否连续,尖峰往往表现为单点跳变后立即回落。可以用下面的思路快速统计:
import numpy as np
p = prof["pressure"]; t = prof["temp"]; s = prof["sal"]
print("缺失数:", np.isnan(p).sum(), np.isnan(t).sum(), np.isnan(s).sum())
dp = np.diff(p)
print("压力回退点数:", (dp < 0).sum(), "最大回退:", dp.min())
# 一阶差分粗筛尖峰:相邻点跳变超过阈值即标记
for name, x, thr in [("T", t, 0.5), ("S", s, 0.1)]:
d = np.abs(np.diff(x))
print(name, "疑似尖峰点:", np.where(d > thr)[0])
示例剖面体检结果:温度缺失 6 点、盐度缺失 9 点,集中在 320~340 dbar 与 1180 dbar 附近;压力序列在 60 dbar 附近出现 3 个回退点;温度一阶差分在 745 dbar 处出现 0.83 ℃ 的跳变,盐度在 745 dbar 处跳变 0.21,属于典型单点尖峰。这些就是后续处理的靶点。
缺失值处理:先判断“该不该补”
缺失分两类。零星单点缺失(连续缺测不超过 3 个采样点)可以用垂向线性插值补齐,因为温盐在垂向上短距离内近似线性;整段缺失(如泵故障导致连续几十个点无盐度)不应插值,应标记为无效段并在后续分析中剔除,否则会人为制造出并不存在的梯度。
示例中 320~340 dbar 的盐度缺失属于连续 9 点,超过阈值,处理方式是保留 NaN 并打上质量标记;1180 dbar 附近的温度单点缺失则用相邻有效点线性插值。插值实现时以压力为自变量而非采样序号,因为采样间隔并不严格等距:
def fill_gaps(p, x, max_gap=3):
x = x.copy(); idx = np.where(np.isnan(x))[0]
if len(idx) == 0: return x
# 按连续缺失段分组
groups = np.split(idx, np.where(np.diff(idx) != 1)[0] + 1)
for g in groups:
if len(g) <= max_gap and g[0] > 0 and g[-1] < len(x) - 1:
x[g] = np.interp(p[g], [p[g[0]-1], p[g[-1]+1]],
[x[g[0]-1], x[g[-1]+1]])
return x
去尖峰:中值滤波加阈值回判
尖峰的判定不能只看差分阈值,否则在温跃层强梯度处会误杀真实信号。稳妥做法是“中值滤波残差 + 物理阈值”双重判据:先用窗口宽度 5 的中值滤波得到局部基线,计算残差 r = x − median(x),当 |r| 超过阈值且该点与前后点符号相反(即跳变后回落)时判为尖峰,再用相邻有效点的线性插值替换。
示例中 745 dbar 处温度残差 0.79 ℃、盐度残差 0.19,均超过阈值(T 取 0.2 ℃,S 取 0.05),且前后点残差符号相反,判定为尖峰并替换。需要提醒的是,阈值必须结合仪器噪声水平设定:本例假设噪声标准差约 0.02 ℃,阈值取 10 倍噪声较为合理;若换用低精度传感器,阈值应相应放宽,否则会把噪声当尖峰反复“削平”,反而损失真实微结构。
压力校正:零点偏移与动态响应
压力通道常见两类误差。第一类是零点偏移:CTD 出水后在空气中读数不为零,记作 p0,整条剖面的压力都应减去 p0。第二类是动态响应滞后:传感器响应时间常数造成压力读数滞后于真实深度,表现为同一等温面上行与下行压力不一致。零点偏移的校正直接做减法:
p0 = np.nanmedian(prof["pressure"][prof["pressure"] < 2]) # 空气中读数
p_corr = prof["pressure"] - p0
示例中空气中读数中位数为 0.6 dbar,校正后表层压力回到 0 附近。动态响应滞后若不做专门实验难以精确标定,工程上常用“上下行剖面比对”估计:同一水团在下降与上升段应具有相同温盐,若压力系统性偏移,则说明存在滞后,可据偏移量做一阶修正。本例数据仅含下降段,故只做零点校正,并在质量标记中注明“未做动态响应校正”。
传感器漂移校正:用稳定水团做参照
盐度漂移是长期观测中最棘手的问题。电导池结垢、电极老化都会让电导率读数缓慢偏移,表现为同一深度的盐度逐年或逐航次系统性变化。校正思路是找一个“温盐稳定”的参照层,例如深层水或温跃层以下的均匀水团,用历史气候值或同航次相邻剖面的参照值做比对,得到偏移量后整段平移。
示例中取 1200~1400 dbar 的深层水作为参照层,该层温度梯度小于 0.005 ℃/dbar,盐度稳定。假设历史气候值为 34.62,本剖面该层盐度均值为 34.58,则漂移量 ΔS = 34.62 − 34.58 = +0.04,整段盐度加 0.04。注意漂移校正应放在去尖峰和缺失处理之后,否则异常点会污染参照层均值;同时校正量应记录在元数据中,便于追溯。
顺序为什么不能颠倒
上述五步的顺序有讲究。缺失处理必须在去尖峰之前,否则插值会把尖峰“抹”进相邻点,形成更难识别的宽峰;去尖峰必须在漂移校正之前,否则参照层均值会被异常点拉偏;压力零点校正应尽早做,因为后续所有以压力为自变量的插值和分层统计都依赖正确的深度坐标。若先做漂移校正再做去尖峰,一个 0.2 的盐度尖峰足以让 200 dbar 参照层的均值偏移 0.001 量级,看似很小,但对精度要求 0.01 的深层水分析已不可忽略。
结果核验:三条硬指标
清洗完成后不能只看图“顺眼”,要核验三条硬指标。第一,压力单调性:校正后下降段压力回退点数应为 0,本例校正前 3 个回退点,扣除零点偏移后仍存在,需结合船体起伏记录判断是否剔除该段。第二,温盐垂向梯度连续性:去尖峰后一阶差分的最大值应回落到噪声量级,本例温度最大跳变由 0.83 ℃ 降至 0.06 ℃。第三,参照层一致性:漂移校正后参照层盐度均值应等于气候值,本例校正后为 34.62,与目标一致。
完成上述链路后,这条剖面才具备进入下一步分析的条件:压力坐标可信、缺失有标记、尖峰已剔除、漂移已补偿。清洗不是一次性动作,而是一套带参数、带日志、可复核的处理流程,后续的T-S图、水团划分与层结计算都建立在这一步的质量之上。