简介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 用于反演 AOD936nm 专用于水汽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 这类偏振通道塞在末尾处理前先进行通道指认否则后续反演时会把偏振值当成直射信号参与运算结果整段作废。参数说明这里 CH6936nm是水汽专用通道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, channel440nm) 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, channelCH4, coeff_threshold0.01, jump_threshold0.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在 Z80° 时两者的偏差可以超过 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.25588H 为海拔单位米。把这两个函数封装好后所有通道共用一套不需要逐波长改参数。3.2 AOD 反演对数比值法与波长指数有了定标常数 V0、大气质量数 m 和瑞利订正项之后AOD 的反演就退化成逐通道的代数运算。代码里最需要注意的地方是通道对应关系440nm 对应 CH1675nm 对应 CH3870nm 对应 CH41020nm 对应 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[faod_{wl}] df.apply( lambda row: calculate_aod(row[f{ch}_mean], v0_dict[wl], row[airmass], wl, pressure), axis1 )逻辑说明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 ), axis1 )逻辑说明代码里把总信号取对数后加上 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 出现一批负值集中在早晨和傍晚。检查发现 CH1440nm的信号在太阳天顶角大于 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 斜率检验反演品质并输出最终 AOD5.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, indexFalse)逻辑说明输出时把波长指数也一起保存因为后续计算细粒子占比时需要用到重新反演一次成本高且容易重复踩坑。quality_flag 来源于云检测、大气质量数和定标偏差三个环节的组合逻辑任何一个环节标记异常最终标志位就取递归最小值。参数说明CSV 输出时建议保留原始时间戳但转成 ISO 格式如果站点后续要和 MODIS 或 CALIPSO 数据对比时间格式统一可以减少大量不必要的麻烦。验收时还有一个习惯把反演结果和同一时段的气溶胶自动监测站点的 PM10 数据做散点对比。AOD 与 PM10 在相对湿度稳定的日子里应当存在正相关相关性很低时优先检查云检测是否漏掉了一批薄云。从那以后我每次处理完 CE318 数据都会强制走一遍“散点图-Langley残差-时间序列平稳度”的三步检查确认没有问题才敢把结果交出去。希望这套流程和踩坑记录能帮你在 CE318 数据处理上少走几个来回。本文还有配套的精品资源点击获取