简介：CE318太阳光度计是大气科学中测量气溶胶与水汽的重要仪器，这份资源面向大气遥感、环境监测方向的研究者和学生，提供了一套基于C++的CE318观测数据处理代码，核心目标是反演气溶胶光学厚度（AOD）与大气水汽含量（WV）。压缩包为zip格式，共5个文件，全部为cpp源码，体积仅9KB，虽然精简却覆盖了主要计算环节，包括原始数据读取、太阳位置与几何订正、AOD反演以及水汽吸收带计算等模块，适合直接阅读、编译和二次开发。目前已有1083人学习，代码结构清晰，配合观测数据即可快速跑通从原始记录到AOD、WV输出的完整流程，免去从零搭建的繁琐。相比零散笔记，这份源码能帮助使用者理解每一步的数据处理逻辑，并针对自身站点数据修改算法参数，实用性强。
1. 太阳光度计数据处理：CE318 从原始 DN 值到 AOD 与水汽的两条反演链路
坐标某区域站点，任务是把 CE318 太阳光度计攒了一个月的原始辐亮度数据，处理成气溶胶光学厚度 AOD 和大气水汽含量 WV。手里的东西是一堆以日期、时间和太阳几何角度开头的文本文件，信号值还是十进制的 DN 数，距离科研需要的物理量中间隔着定标、大气质量数、瑞利散射订正、云检测和逐通道回归五六道工序。这套东西最反直觉的地方在于：AOD 和水汽看似是两个完全独立的物理量，在 CE318 上却共用同一批太阳直射信号，区别只在选哪个波长做“分母”。本资源就是一套从原始文件读入到最终反演结果导出的完整处理流程，附带可改参数的脚本和逐段排错记录，适合遥感、大气环境方向的从业者和研一新生，也适合手头攒着几年老数据却不知道从哪下手的站点管理员。
2. 前置条件与数据链路：Langley 定标、通道指认和云检测的三道门槛
2.1 认识 CE318 的通道结构与数据文件
CE318 是太阳光度计里最常见的一类型号，多数配置下有 6 个太阳直射通道和多个散射通道。太阳直射通道里，440nm、675nm、870nm、1020nm 用于反演 AOD，936nm 专用于水汽，500nm 常作为辅助可见光通道。散射通道（如 440nm、870nm 等）主要做大气气溶胶相函数和粒径分布，但不下放给本流程。
数据文件的典型行格式如下（不同批次仪器导出列数可能有差异，但顺序大体一致）：
20250315_084532 38.21 146.52 34122 41008 35644 28810 22103 24750 15122 14800 20250315_084702 38.24 146.60 34088 40961 35592 28744 22055 24701 15089 14762 20250315_084832 38.27 146.68 34011 40845 35498 28691 21984 24610 15033 14699逻辑说明：第一列是时间戳，第二列是太阳天顶角（单位度），第三列是太阳方位角，后续从第四列开始的数值依次对应 CH1 到 CH8 的原始 DN 值（或电压值）。部分仪器导出时会把 CH7、CH8 这类偏振通道塞在末尾，处理前先进行通道指认，否则后续反演时会把偏振值当成直射信号参与运算，结果整段作废。
参数说明：这里 CH6（936nm）是水汽专用通道，CH7、CH8 如果仪器没有偏振观测模块，默认不参与反演。读取时建议统一保留原始列，不要只保留需要的通道。很多站点数据是从老机器里拷出来的，列数会少一两列，保留原始列更容易对齐。
我一般会在读文件后立即打印各列均值，检查 CH1-CH6 的量级是否在几千到几十万之间。如果某一列出现整段恒为零或者量级比其他通道小两个数量级，先查仪器探头，不要直接进入反演。
2.2 定标：Langley 斜率法的原理与适用窗口
AOD 反演的物理基础是比尔-朗伯定律：太阳辐射穿过大气后的衰减与大气光学厚度成正比。仪器测得的 DN 值与大气外太阳辐照度的定标常数 V0 之间存在指数衰减关系。定标就是要确定这个 V0。实际操作中，AOD 反演公式可以写成：
AOD(λ) = (lnV0(λ) - lnDN(λ) - τ_rayleigh(λ) × m - τ_其它气体(λ) × m) / m
其中 m 是大气质量数。如果我们在一天中选取气溶胶稳定的时段，对 lnDN 和 m 做线性回归，截距就是 lnV0，斜率反映了大气总光学厚度。这就是 Langley 定标。
注意一个关键前提：Langley 定标要求大气状态在观测窗口内保持“清洁且平稳”，通常选在山区站点、连续几天无云且气溶胶浓度极低的日子里做。城市站点数据直接拿某一天的整段数据回归，会因为气溶胶和云的快速变化产生巨大误差。常见做法是把多天的清晨数据点绘在同一坐标系里，选取 lnDN 随 m 线性度最好的一段进行拟合。
import numpy as np def langley_calibrate(m, ln_dn, channel): # m 为大气质量数序列，ln_dn 为对应通道的对数辐亮度序列 fit = np.polyfit(m, ln_dn, 1) slope, intercept = fit v0 = np.exp(intercept) return slope, v0 m_selected = m[(m >= 2.0) & (m <= 6.0)] dn_selected = ln_dn[(m >= 2.0) & (m <= 6.0)] slope, v0 = langley_calibrate(m_selected, dn_selected, channel='440nm') print(slope, v0)逻辑说明：这里强制限制大气质量数在 2.0 到 6.0 之间，是因为 m 太小时太阳高度角大，对指向精度极敏感；m 太大时光程长，容易被薄云和低层气溶胶污染。截取中间段做回归，是朗格莱法的经典做法。
参数说明：m 下限和上限可视站点海拔调整。高海拔站点可以把下限放到 1.5，污染严重的城市站点建议下限拉到 2.5 以上。计算得到的 v0 会和仪器出厂定标文件里的值做对照，偏差在 3% 以内视为合格。
2.3 云检测与数据质量控制
云是太阳光度计反演里最大的威胁。单层薄云会导致信号缓慢衰减，反演出的 AOD 会明显偏高；碎云会造成信号陡峭跳动，结果则随机漂移。不知道什么时候该删点，是新手处理 CE318 数据时最常见的翻车点。
常用做法是分两级检测。第一级：对每分钟的 3 个采样点求标准差，如果同一通道的变异系数超过 1%，则该分钟数据标记为可疑。第二级：对滑动窗口内的 10 分钟序列做线性趋势检验，如果相邻时刻的信号跳变超过整体均值的 5%，则定位为碎云穿过。
def cloud_filter(df, channel='CH4', coeff_threshold=0.01, jump_threshold=0.05): df['cv'] = df[f'{channel}_std'] / df[f'{channel}_mean'] df.loc[df['cv'] > coeff_threshold, 'flag'] = 'cloud' df['diff'] = df[f'{channel}_mean'].diff().abs() df.loc[df['diff'] / df[f'{channel}_mean'].shift(1) > jump_threshold, 'flag'] = 'broken_cloud' return df[df['flag'].isna()]逻辑说明：先按变异系数过滤，再按时间连续性过滤，两级顺序不能反。因为碎云出现时变异系数可能不大，但时间序列上会出现陡升陡降；反过来，如果只做时间连续性检测，又会把仪器瞬时指向偏差当成云剔除，导致可用数据量锐减。
参数说明：变异系数阈值 1% 适用于大多数晴天。沙尘天或低太阳高度角场景下，辐亮度本身波动大，可以放宽到 2%。跳变阈值 5% 是经验值，我一般会用半天数据的 95 分位来动态设定，自适应效果比固定值好得多。
3. 反演实现：AOD 逐通道计算与水汽含量的改进 Langley 法
3.1 大气质量数与瑞利散射订正
AOD 计算里的 m（大气质量数）并不是简单的 1/cosZ。低太阳高度角下，大气曲率效应会让实际光程比几何光程短，必须用 Kasten-Young 公式修正：
m = 1 / (cosZ + 0.50572 × (96.07995 - Z)^(-1.6364))
其中 Z 是太阳天顶角。用这个公式替代 1/cosZ，在 Z=80° 时两者的偏差可以超过 10%，直接体现在 AOD 本底上。
瑞利散射订正项也不可跳过。大气分子散射在 440nm 处贡献的光学厚度大约在 0.2 以上，大于很多干净地区的气溶胶光学厚度。简化形式里，常用如下经验公式：
τ_rayleigh(λ) = 0.008569 × λ^(-4) × (P / 1013.25)
其中 P 为站点气压（单位 hPa），λ 单位微米。这个公式对 440nm 到 1020nm 通道的精度已经足够，不需要动用逐线积分。
import math def airmass(zenith_deg): z = math.radians(zenith_deg) cosz = math.cos(z) return 1.0 / (cosz + 0.50572 * (96.07995 - zenith_deg) ** (-1.6364)) def rayleigh_od(wavelength_nm, pressure_hpa): lam_um = wavelength_nm / 1000.0 return 0.008569 * lam_um ** (-4) * (pressure_hpa / 1013.25)逻辑说明：airmass 函数把度转换为弧度后先算余弦，再用 Kasten-Young 公式修正，避免低角度失真。rayleigh_od 则是把波长从纳米换算成微米，再按气压比例伸缩，返回无量纲的光学厚度。
参数说明：如果手头只有海拔数据没有实测气压，用标准大气近似：P = 1013.25 × (1 - 2.25577e-5 × H)^5.25588，H 为海拔（单位米）。把这两个函数封装好后，所有通道共用一套，不需要逐波长改参数。
3.2 AOD 反演：对数比值法与波长指数
有了定标常数 V0、大气质量数 m 和瑞利订正项之后，AOD 的反演就退化成逐通道的代数运算。代码里最需要注意的地方是通道对应关系：440nm 对应 CH1，675nm 对应 CH3，870nm 对应 CH4，1020nm 对应 CH5。写脚本时用一个字典维护通道波长映射，比硬编码列号安全得多，处理多台仪器数据时尤其省事。
def calculate_aod(dn, v0, m, wavelength_nm, pressure_hpa): tau_rayleigh = rayleigh_od(wavelength_nm, pressure_hpa) tau_total = (math.log(v0) - math.log(dn)) / m aod = tau_total - tau_rayleigh return max(aod, 0.0) channels = {'CH1': 440, 'CH3': 675, 'CH4': 870, 'CH5': 1020} for ch, wl in channels.items(): df[f'aod_{wl}'] = df.apply( lambda row: calculate_aod(row[f'{ch}_mean'], v0_dict[wl], row['airmass'], wl, pressure), axis=1 )逻辑说明：calculate_aod 里先用定标常数和实测信号的对数差除以大气质量数，得到大气总光学厚度，再减去瑞利散射贡献，剩余量就是气溶胶光学厚度。max 函数把负值强制归零，实际工作中偶尔会出现 AOD 恰好等于零的极端干净情况，但大多数负值源于定标不准或云检测漏网，建议同时记录原始负值，不要直接抹掉。
参数说明：v0_dict 是各波长的定标常数集合，由第 2 章的 Langley 定标结果生成。pressure 参数只对瑞利项有影响，在高原站和海平面站之间，该项差异可达 20%。
反演完成后，顺手计算 440nm 与 870nm 之间的 Angstrom 波长指数 α：
α = ln(AOD₁₄₄₀ / AOD₈₇₀) / ln(870 / 440)
α 值落在 0.8 到 1.8 之间表示以细粒子为主，接近 0 则表示以沙尘等粗粒子为主。这个指数可以作为下一个质量控制点：如果 α 出现负值，说明两个通道之一的反演结果异常，优先回到云检测步骤排查，而不是继续强行解释。
3.3 936nm 水汽反演：改进 Langley 法与透过率转换
水汽反演和 AOD 的不同之处在于：936nm 通道同时受气溶胶衰减和水汽吸收影响，没法只靠一条比尔-朗伯定律一次解出。标准做法是先用其它通道插值出 936nm 处的气溶胶光学厚度，把它从总衰减里全部扣除，剩下的衰减就归因于水汽吸收；再用改进的 Langley 法处理吸收通道，得到可降水量 WV。
这里有一个极其重要的坑：936nm 落在水汽吸收带内，仪器响应对水汽总量高度敏感，所以不能用整段数据做 Langley 回归。改进 Langley 法只选取水汽吸收较弱的观测点——通常是通过大气质量数和水汽总量的双条件筛选，只保留大气质量数适中的数据点参与回归，避免水汽饱和段污染斜率。
def water_vapor_retrieval(dn_936, aod_936, m, v0_936): # 扣除气溶胶衰减后的“净水汽透过率”的对数值 ln_T_wv = math.log(dn_936) + m * aod_936 - math.log(v0_936) # 近似线性关系：ln(T_wv) 与 (m * sqrt(WV)) 成正比 sqrt_wv = abs(ln_T_wv) / (m * 0.93) return sqrt_wv ** 2 df['wv'] = df.apply( lambda row: water_vapor_retrieval( row['CH6_mean'], row['aod_936'], row['airmass'], v0_936 ), axis=1 )逻辑说明：代码里把总信号取对数后，加上 m × aod_936（等价于把被气溶胶吃掉的那部分“补回去”），再减去定标常数，剩下 ln(T_wv) 就纯粹是水汽吸收的贡献。0.93 是从 HITRAN 谱线数据库按 CE318 936nm 通道带宽做平均得到的吸收系数近似值，在 0.5 到 5cm 的可降水量范围内误差可控制在 10% 以内。
参数说明：v0_936 是从改进 Langley 定标中获得的，它的物理意义是“大气外水汽通道等效辐亮度”，并不对应真实大气顶辐亮度，所以同一个仪器上，v0_936 和其它通道的 V0 量级可能差很多，属于正常现象。
aod_936 的获得方式是从相邻通道插值：AOD₉₃₆ = AOD₈₇₀ × (936 / 870)^(-α)，其中 α 来自 440nm 与 870nm 的波长指数。如果当天数据 α 无意义，直接用 AOD₈₇₀ 代替 AOD₉₃₆，误差通常在可接受范围内。
4. 避坑：CE318 数据处理五个高发错误与排查
4.1 现象一：8 个通道里无缘无故多出负值
某天反演后 AOD 出现一批负值，集中在早晨和傍晚。检查发现 CH1（440nm）的信号在太阳天顶角大于 75° 时极低，甚至低于仪器的暗电流噪声值，导致对数相减后变成负的“光学厚度增量”。
原因：低太阳高度角时，大气质量数大，440nm 的瑞利散射和气溶胶衰减联合起来把信号压到了仪器的灵敏度下限以下。这是物理极限，不是算法缺陷。
解决：处理代码里对 m > 6 的数据直接标记为无效，不参与 AOD 反演。同时保留原始负值日志，便于做数据质量统计时判断截断比例是否过高。如果某天超过 30% 的数据被标记无效，建议只保留 10:00 到 15:00 的时段。
4.2 现象二：Langley 定标的截距离谱
在站点做 Langley 定标，得到 440nm 通道的 V0 比仪器出厂证书值高了 15%，重做两次依然如此。
原因：定标时段选在午后，站点处在下风向，正在经历一次轻度沙尘输送。气溶胶光学厚度并非恒定，而是在缓慢上升，Langley 曲线整体弯曲，线性拟合的截距被抬高。
解决：换成连续三天的清晨数据进行联合定标。数据点不满足线性段标准时，肉眼判断曲线形态：如果 lnDN 对 m 散点图出现系统性穹顶形状，说明气溶胶在加速变化，丢弃这段数据。从那以后，我每次做 Langley 定标都会先画散点图，再谈拟合。
4.3 现象三：水汽反演跳变，结果忽大忽小
某日 936nm 水汽反演结果从 2.5cm 突然跳到 4.8cm，又跌回 2.3cm，逐分钟曲线呈锯齿状。
原因：碎云干扰未被云检测识别出来。936nm 通道对水汽吸收敏感，云中液态水同样会吸收辐射，信号变化极大，碎云导致透过率被低估，反演出偏高的水汽量。
解决：将云检测的变异系数阈值从 1% 调整到 0.5%，同时对水汽通道单独做了一个额外检测：当相邻 5 分钟的 WV 变化超过 20% 时，强制标记该点。简单粗暴，但能有效去掉碎云点。
4.4 现象四：波长指数几乎恒等于零
一批 440nm 和 870nm 的 AOD 反演结果都比正常值偏大，算出的波长指数 α 长时段接近 0，看似是沙尘事件，但同期地面观测并无沙尘。
原因：500nm 通道与 440nm 通道的信号被写反了（数据文件列错位），440nm 通道使用的是 500nm 的信号与定标常数，导致总光学厚度偏大。列错位问题在数据导出环节最容易发生。
解决：把反演前的各通道定标常数 V0 打印出来对照检查。V0 随波长变化应该单调，440nm 的 V0 必然大于 500nm 的 V0，一旦出现交叉，立刻检查数据列映射。
4.5 现象五：数据文件时间戳错位，反演结果错峰
所有通道反演结果看起来正常，但和相邻站点的对比序列整体偏移了 15 分钟左右。
原因：部分 CE318 数据文件里的时间戳是仪器内部时钟的“采集开始时间”，而信号值是采集结束时刻的读数，两者之间存在一个扫描周期的时间差。处理时直接用时间戳和太阳几何角度匹配，导致整段序列错位。
解决：读取文件时提取时间戳和太阳天顶角的对应关系，对天顶角做插值，以时间戳为准重新对齐采集结束时间。对于固定扫描周期为 3 分钟的站点，直接把时间戳加上周期值即可恢复真实时刻。
5. 验证与进阶：用 Langley 斜率检验反演品质并输出最终 AOD
5.1 用定标残差和光谱平滑性做三位一体验证
反演完成后，第一件事不是直接出图，而是用三种廉价方式做交叉验证：
- 检查 500nm 通道和 440nm 通道的 AOD 差值。整套流程不反演 500nm，但可以用插值法从 440nm 和 675nm 推出 500nm 的 AOD，再和实测 500nm 仪器通道（如果数据里有）对比，偏差超过 0.02 就要怀疑瑞利订正或定标常数。
- 对全天的 AOD 逐小时均值做调和分析，正常晴天气溶胶日变化呈平滑 U 形或倒 V 形，出现高频抖动时优先怀疑云检测。
- 把当天的 440nm 反演 AOD 和 Langley 定标回归残差放在同一张图上。残差大的点反演结果一般不靠谱。
5.2 批量输出与水汽结果的质量标识
处理完整月数据时，我会把反演结果统一输出为 CSV，同时附带一列 quality_flag 标志位，取值为 0（正常）、1（云污染）、2（低太阳高度角）、3（定标可疑）。这一步对后续做气候统计分析特别有用——没有质量标识的 AOD 序列，在按月平均时会混入大量污染点。
output = pd.DataFrame({ 'timestamp': df['time'], 'AOD_440': df['aod_440'], 'AOD_675': df['aod_675'], 'AOD_870': df['aod_870'], 'AOD_1020': df['aod_1020'], 'WV_cm': df['wv'], 'alpha': df['alpha'], 'quality_flag': df['quality_flag'] }) output.to_csv('CE318_processed_month.csv', index=False)逻辑说明：输出时把波长指数也一起保存，因为后续计算细粒子占比时需要用到，重新反演一次成本高且容易重复踩坑。quality_flag 来源于云检测、大气质量数和定标偏差三个环节的组合逻辑：任何一个环节标记异常，最终标志位就取递归最小值。
参数说明：CSV 输出时建议保留原始时间戳但转成 ISO 格式；如果站点后续要和 MODIS 或 CALIPSO 数据对比，时间格式统一可以减少大量不必要的麻烦。
验收时还有一个习惯：把反演结果和同一时段的气溶胶自动监测站点的 PM10 数据做散点对比。AOD 与 PM10 在相对湿度稳定的日子里应当存在正相关，相关性很低时，优先检查云检测是否漏掉了一批薄云。从那以后，我每次处理完 CE318 数据，都会强制走一遍“散点图-Langley残差-时间序列平稳度”的三步检查，确认没有问题才敢把结果交出去。希望这套流程和踩坑记录能帮你在 CE318 数据处理上少走几个来回。
本文还有配套的精品资源，点击获取
