📚 课程目录第4章数据清洗与质控的Python基础工具

4. 数据清洗与质控的Python基础工具

海洋结构监测系统（如导管架平台应变监测、浮式平台运动监测、海底管道振动监测）通常以固定采样频率连续记录数据，单日数据量可达数十万至数百万条。面对这样的数据规模，手工处理既不现实也不可靠。在正式进入缺失值填补、异常识别、平滑滤波等清洗任务之前，需要先搭建一套稳定的编程环境，并掌握三个核心库的基本操作：NumPy负责数值数组与向量化计算，Pandas负责带时间标签的表格数据组织与聚合，Matplotlib负责把数据形态和清洗效果画出来。本节以一段模拟的海洋平台应变监测数据为例，走通“读取—切片—聚合—可视化”的完整链路。

一、环境准备与数据对象约定

推荐使用 Anaconda 创建独立环境，避免不同项目的库版本互相干扰：

conda create -n marine_qc python=3.11
conda activate marine_qc
conda install numpy pandas matplotlib

本节使用的示例数据为合理假设的模拟数据：某导管架平台在 2024-06-01 00:00 至 2024-06-02 00:00 期间，以 1 Hz 采样频率记录 4 个应变测点（SG01～SG04）的应变值（单位 με），并附带一个环境温度通道。数据以 CSV 文件保存，字段为 timestamp, SG01, SG02, SG03, SG04, temp。需要强调：以下所有数值结果均为示例演示，不代表任何实际工程测量值。

二、NumPy：数值数组与向量化运算

NumPy 的核心是 ndarray，它把同类型数值连续存放在内存中，使整列运算不必写循环。读取 CSV 的数值部分后，可以先做几件清洗前必须确认的事：

import numpy as np

strain = np.loadtxt("platform_strain_20240601.csv",
                    delimiter=",", skiprows=1, usecols=(1, 2, 3, 4))
print(strain.shape)        # (86400, 4)
print(strain.dtype)        # float64

# 基本统计：逐通道均值、标准差、极值
print(np.mean(strain, axis=0))
print(np.std(strain, axis=0))
print(np.min(strain, axis=0), np.max(strain, axis=0))

# 布尔索引定位疑似异常（示例阈值：|应变| > 300 με）
mask = np.abs(strain) > 300
print(mask.sum(axis=0))    # 各通道超阈值点数

这里 axis=0 表示沿时间方向压缩，得到每个测点的统计量；mask 是一个与数据同形状的布尔数组，可直接用于计数或后续替换。向量化写法 np.abs(strain) > 300 一次完成 86400×4 个元素的比较，比 Python 循环快一到两个数量级。此外，np.isnan(strain).sum() 可快速统计缺失数量，np.diff(timestamp_seconds) 可检查采样间隔是否恒为 1 s——这两项是质控前的例行检查。

三、Pandas：时间索引、切片与聚合

NumPy 数组没有时间标签，而监测数据的清洗几乎都围绕时间展开，因此实际工作中以 Pandas 为主入口：

import pandas as pd

df = pd.read_csv("platform_strain_20240601.csv",
                 parse_dates=["timestamp"])
df = df.set_index("timestamp").sort_index()
print(df.index.freq)       # 若为 None，说明存在缺测或时间跳变

切片方面，Pandas 支持按时间字符串直接截取区间，例如 df.loc["2024-06-01 08:00":"2024-06-01 10:00"] 取出两小时数据，df["SG01"] 取出单列，df.iloc[0:3600] 按位置取前 3600 行。若需按固定窗口统计，用 resample：

# 10 分钟均值与标准差，用于观察趋势和波动
stat10 = df.resample("10min").agg(["mean", "std", "count"])
print(stat10["SG01"].head())

# 按小时统计各通道最大值与缺失数
hourly = df.resample("1h").agg({"SG01": "max", "SG02": "max",
                                "SG03": "max", "SG04": "max"})
print(hourly)

count 列尤其重要：若某 10 分钟窗口的 count 明显小于 600，说明该时段存在缺测，这正是后续清洗要处理的对象。agg 支持对不同列指定不同函数，也可传入自定义函数，例如用 lambda x: x.quantile(0.99) 统计高分位数，用于判断冲击型异常。

四、Matplotlib：把数据形态画出来

清洗前先画图，是判断数据“长什么样”的最快方式。下面同时画出原始时序与 10 分钟均值：

import matplotlib.pyplot as plt

fig, ax = plt.subplots(2, 1, figsize=(11, 6), sharex=True)
ax[0].plot(df.index, df["SG01"], lw=0.4, color="#2b7bba")
ax[0].set_ylabel("应变 / με")
ax[0].set_title("SG01 原始时序（示例数据）")

ax[1].plot(stat10.index, stat10[("SG01", "mean")],
           color="#e07b39", marker="o", ms=2)
ax[1].set_ylabel("10 min 均值 / με")
ax[1].set_xlabel("时间")
plt.tight_layout()
plt.savefig("sg01_overview.png", dpi=150)

观察要点：原始曲线若出现孤立尖峰（单点跳变后立即恢复），多为传感器噪声或瞬时干扰；若出现持续平台段或长时间恒定值，可能是信号饱和或采集卡死；若 10 分钟均值曲线在某时段突然抬升，则需回到原始数据核对该时段是否存在阶跃。这些判断都依赖图形，而非单纯统计量。

五、把四步串成一条可复用的检查脚本

将上述操作组织成一个函数，每次拿到新数据先跑一遍，能快速形成对数据质量的初步判断：

def quick_check(path, channels=("SG01", "SG02", "SG03", "SG04")):
    df = pd.read_csv(path, parse_dates=["timestamp"]).set_index("timestamp")
    df = df.sort_index()
    report = {}
    report["行数"] = len(df)
    report["时间跨度"] = df.index[-1] - df.index[0]
    report["缺失数"] = df[list(channels)].isna().sum().to_dict()
    report["超阈值数"] = (df[list(channels)].abs() > 300).sum().to_dict()
    win = df[list(channels)].resample("10min").count()
    report["不完整窗口数"] = int((win < 600).any(axis=1).sum())
    return report

print(quick_check("platform_strain_20240601.csv"))