简介：面向大气遥感与太阳光度计数据处理方向的学习者，这份资源围绕CE318型太阳光度计的观测数据展开，重点解决如何从原始测量中反演气溶胶光学厚度（AOD）与大气水汽含量（WV）这一实际问题，适合具备一定大气光学基础、希望打通从数据读取到参数反演完整链路的高年级本科生、研究生及科研人员参考。压缩包共5个文件，均为cpp源码，体积约9KB，分别对应CE318数据读取、AOD对比计算、水汽对比计算、定标处理与角度拟合等环节，可直接作为反演流程的代码骨架阅读与复用。资源已有1084人学习下载，具备一定参考基础。读者可从中了解多通道滤光片观测数据的组织方式、异常值与缺失值清洗思路，以及AOD与水汽反演中辐射传输、吸收带选取、气象因素修正和定标校准的实现要点，并借助对比计算模块开展结果质量控制与趋势分析，为后续科研与环保监测应用提供可落地的脚本起点。
1. CE318 反演 AOD 与 wv：从一天的原始观测到一个能用的数据点
CE318 太阳光度计到手后，很多人的第一反应是打开随机软件点导出，拿到一列 AOT 就以为完事了。真正做反演的人清楚，从原始 DN 到一条可信的 AOD 曲线，中间隔着暗电流、V0 定标、空气质量数和臭氧订正四道关；再到 940 nm 的水汽总量 wv，还得额外把气溶胶贡献从吸收里剥出来。这条链路适合两类人：手上有仪器原始数据、想自己控一遍流程的观测人员，以及需要 AOD 和 wv 作为辐射传输或气溶胶反演输入量的算法工程师。下面按读文件、定标、AOD、wv、排错的顺序，把每一步的代码和判据都写死。
2. 读懂 CE318 观测文件与 V0 定标：反演前的两笔账
2.1 一次观测记录里到底有什么
CE318 的原始数据按天存放，一天一个二进制文件，里面是一串按时间排列的 measurement sequence。每个直射太阳 sequence 大致包含：一次暗电流采样、八个波段的太阳电压（DN）、太阳天顶角与方位角、四象限瞄准信号、探头温度、箱体温度、气压、GPS 时间戳。同一个 sequence 通常连续测三次，间隔几十秒，业内叫 triplet。这三次的离散度是后面判断有没有薄云最便宜的依据，别丢掉。
天空辐射（almucantar、principal plane）的记录混在同一个文件里，用测量类型字段区分。反演 AOD 只需要直射太阳那一类，先过滤再谈后面的事。八个标准通道的分工大致如下：
| 中心波长 | 主要用途 | 备注 |
|---|---|---|
| 340 nm | 短波端气溶胶、细模态 | 信号弱，清洁天噪声大 |
| 380 nm | 与 340 配合判吸收 | 受臭氧影响可忽略 |
| 440 nm | 气溶胶主通道 | 有微量二氧化氮吸收 |
| 500 nm | 与 440 构成短波对 | 报告 AOD 的常用参考波长 |
| 675 nm | 与 440 算 Ångström 指数 | 中段最稳的一档 |
| 870 nm | 与 440 配对、给 940 插值 | 长波端主力 |
| 940 nm | 水汽吸收通道 | 也出 AOD，但一般不用于气溶胶分析 |
| 1020 nm | 长波端、粗模态 | 信号弱，最容易被暗电流坑 |
940 nm 这一档要单独说：软件也会给你一个 AOD 值，但它对水汽吸收系数的误差极其敏感，拿它去做气溶胶谱反演基本是自找麻烦。它的正确用途是算 wv。
2.2 导出 CSV 与字段对齐：别写死 skiprows
用厂家那套处理软件（常见的是 ASTPWin 那一套）把二进制导成 CSV，选 Level 1.0 就行。导的时候务必勾上原始电压列和暗电流列，只导 AOT 列的话，后面自己重新定标、自己重算 Langley 都没法做，等于把黑匣子焊死了。
读文件这一步有个小坑：不同版本导出的表头行数不一样，写死skiprows=5的脚本换个版本就全线报错。稳妥的写法是找表头：
from pathlib import Path import pandas as pd def load_level1_csv(path): """读 CE318 处理软件导出的 Level 1.0 CSV，自动定位表头行""" lines = Path(path).read_text(encoding="utf-8", errors="ignore").splitlines() # 表头不一定在第一行，找同时含 Date 和 Time 的那一行 header_idx = next(i for i, line in enumerate(lines) if "Date" in line and "Time" in line) df = pd.read_csv(path, skiprows=header_idx, encoding="utf-8", errors="ignore") # 列名里的空格会让 df.AOT_440 这种写法失效，先统一掉 df.columns = [c.strip().replace(" ", "_") for c in df.columns] # 时间统一存 UTC，中间环节一律不碰本地时 df["dt_utc"] = pd.to_datetime( df["Date(dd:mm:yyyy)"] + " " + df["Time(hh:mm:ss)"], format="%d:%m:%Y %H:%M:%S", errors="coerce") return df.dropna(subset=["dt_utc"]).reset_index(drop=True)这段代码的逻辑很直白：先按行扫描找列名行，再按列名解析。参数上唯一需要改的是时间格式串——如果你的软件导出的日期是yyyy-mm-dd，把format换成对应写法即可，但一定要显式指定，让 pandas 去猜格式在跨年份数据上会慢十倍且容易猜错。
单位也要在这一步统一：气压用 hPa，波长用 nm，时间用 UTC，电压用 DN（无量纲计数）。这三样只要有一处混用，后面所有订正项都会带上一个说不清来源的偏差。
2.3 暗电流、饱和与指向：三个一秒筛掉坏观测的量
暗电流是每个 sequence 开头单独采的那一帧，探测器在无光条件下的输出。它不是常量，随探头温度漂移，夏天午后能到几十 counts。1020 nm 通道的太阳信号本身只有几千 counts，暗电流扣不扣，差别能到百分之一以上，落到 AOD 上就是 0.01 量级的系统性偏差——在清洁站，这个数字已经和真实气溶胶信号同一量级了。
第二个是饱和。DN 是十六位计数，理论上限 65535，但探测器的线性区到不了那么高，接近上限时响应会压平，算出来的 AOD 假性偏低甚至为负。我一般在满量程的 5% 处就设闸。
第三个是瞄准。四象限信号不平衡说明 seek/aim 没跟上太阳，低太阳高度角时这个误差会被气团质量放大。天顶角超过 80° 的观测，气团质量已经到 5.8 以上，折射和站址误差都开始显形，除非专门做 Langley，否则建议直接丢掉。
import numpy as np def screen_observations(df, dn_col="DN_1020", dark_col="DN_dark", sza_col="SZA", triplet_std_col="AOD_triplet_std"): keep = ( (df[dark_col].abs() < 50) & # 暗电流绝对值超 50 counts 的丢掉 (df[dn_col] < 60000) & # 留 5% 线性余量，防饱和 (df[sza_col] < 80) & # 气团质量 < 5.8 (df[triplet_std_col] < 0.02) # 三连测 AOD 标准差，粗筛云 ) return df[keep].copy()这几个阈值不是拍脑袋：暗电流 50 counts 大致对应 1020 nm 上 0.005 的 AOD 偏差，属于可接受上限；triplet 标准差 0.02 是文献里最常用的云筛门槛，薄卷云往往能把三次观测拉开 0.05 以上。如果你的站点气溶胶本底很低，可以把 0.02 收紧到 0.01，代价是丢掉不少有效点，这个取舍得看下游用途。
2.4 Langley：把 V0 从自己站点的数据里拟合出来
定标系数 V0 是整条链路精度的天花板。厂家证书给的 V0 会随时间老化，尤其是 940 nm 和 1020 nm 这两档，滤光片和探测器一起漂，几年下来偏几个百分点很正常。
Langley 法的出发点很简单：直射太阳的电压满足
V = V0 · d⁻² · exp(-m · τ)
其中 d 是日地距离（天文单位），m 是气团质量，τ 是整层大气光学厚度。把 V 折算到 1 AU 之后取对数，ln(V·d²) 对 m 就是一条直线，截距是 ln V0，斜率是 -τ。前提是观测窗口内 τ 恒定——所以只能挑清晨那半小时到一小时，或者高海拔站天气特别稳的时段。
import numpy as np def earth_sun_distance_sq(doy): """返回 (R/R0)^2，即 d^-2，doy 为年积日""" g = 2 * np.pi * (doy - 1) / 365.25 return (1.00011 + 0.034221 * np.cos(g) + 0.00128 * np.sin(g) + 0.000719 * np.cos(2 * g) + 0.000077 * np.sin(2 * g)) def kasten_young_airmass(sza_deg): """Kasten-Young 大气质量数，sza_deg 为太阳天顶角（度）""" sza = np.radians(sza_deg) return 1.0 / (np.cos(sza) + 0.50572 * (96.07995 - sza_deg) ** -1.6364) def langley_fit(dn, sza_deg, doy, min_m=1.5, max_m=5.0): m = kasten_young_airmass(sza_deg) d2 = earth_sun_distance_sq(doy) # d^-2 ok = (m >= min_m) & (m <= max_m) y = np.log(dn[ok] * d2[ok]) # 折算到 1 AU 的 lnV x = m[ok] A = np.vstack([x, np.ones_like(x)]).T slope, intercept = np.linalg.lstsq(A, y, rcond=None)[0] resid = y - (slope * x + intercept) r2 = 1 - np.sum(resid ** 2) / np.sum((y - y.mean()) ** 2) return np.exp(intercept), -slope, r2, np.std(resid)参数说明：min_m和max_m圈定参与拟合的气团质量区间，1.5 到 5.0 是常见取法，下限防止太阳太高时信号饱和，上限防止低高度角的折射误差。返回值依次是 V0、整层光学厚度、拟合 R² 和残差标准差。
判据怎么读：R² 一般要求 0.9995 以上，残差标准差小于 0.01（对应约 1% 的电压精度），拟合出的 τ 要和站点气候背景值对得上——清洁站 500 nm 应该在 0.02 到 0.1 之间，拟合出 0.5 说明这段数据里混了云。
同一天不同日期重复做 Langley，会得到一组 V0。把它们按时间画出来，一年内漂移超过 1% 就要查滤光片老化和温度系数；和证书值比偏差超过 2%，优先怀疑站址气压填错、暗电流没扣干净。至于 940 nm，普通 Langley 在它身上基本做不出直线——水汽变化太快，只能用改进 Langley 或者传递定标，这个后面单独说。
3. AOD 反演链路：气团质量、瑞利与臭氧三项订正怎么落到代码
3.1 从电压到光学厚度：公式是怎么拆的
把 2.4 节那个式子反过来用。观测电压对应的是整层大气所有消光成分的总和：
V = V0 · d⁻² · exp(-m · (τ_aer + τ_R + τ_O3 + τ_NO2))
取对数整理，气溶胶光学厚度就是总光学厚度减掉其余几项：
τ_aer = ln(V0 · d⁻² / V) / m - τ_R - τ_O3 - τ_NO2
四项里，τ_R 是瑞利散射，干净天在 440 nm 能到 0.2 以上，是最大的一项订正；τ_O3 是臭氧吸收，在 500 到 700 nm 之间最显著，量级 0.01 到 0.03；τ_NO2 只在 440 nm 附近有点存在感，通常几个 0.001，要求不高可以忽略。这三项都不是气溶胶，但都会被算进总消光里，少减一项，AOD 就整体偏高一个常数。
3.2 瑞利与臭氧：两项订正的参数从哪来
瑞利光学厚度有个经典近似：
τ_R(λ) = 0.008735 · λ^(-4.08) · (P / P0)
λ 用微米，P 是站址气压，P0 取 1013.25 hPa。这里最容易翻车的是气压：高海拔站点如果不做气压订正，直接用海平面气压，瑞利项会系统性偏大，AOD 被减多，整条曲线偏低 0.01 到 0.02。这个量级在清洁地区已经足以改变结论。站址气压要么用仪器自带的测量
本文还有配套的精品资源，点击获取
