简介在大气气溶胶与温室气体研究中AOD气溶胶光学厚度和WV水汽含量是反映大气环境状态的重要参数CE318太阳光度计通过多波段太阳辐射观测为反演这些参数提供原始数据。一套聚焦CE318数据处理的C源码资源可用于实现从观测数据到AOD、WV反演的完整流程主要面向大气科学、遥感应用及环境监测方向的研究生、科研人员与相关课程学习者。压缩包共5个文件均为C源码大小仅9KB分别覆盖原始观测数据读取、太阳几何位置解算、Angström指数拟合、气溶胶光学厚度计算和水汽反演等功能模块代码紧凑、结构清晰按处理流程拆分既可对单环节单独调用也能串联成完整反演链。目前已有1082人学习下载。借助该资源读者不仅能从源码层面复现CE318观测数据的处理流程还能理解辐射传输、气溶胶散射吸收以及水汽吸收带选择等核心算法并结合自身站点数据调整参数后投入使用为科研实验、毕业设计或教学演示提供一份轻量且可扩展的参考实现。1. 从CE318原始电压到AOD和水汽WV数据处理链路到底卡在哪一个观测站放着一台CE318太阳光度计采集了一整年的原始电压数据报文里每行是十几个通道的电压读数看起来规整但真要把它变成逐分钟的气溶胶光学厚度AOD和水汽柱含量WV中间隔着一长串步骤定标、日地距离修正、云筛选、瑞利散射扣除、臭氧吸收订正最后还要对936nm水汽通道做非线性拟合。不少做气溶胶观测的从业者拿到一套CE318处理脚本和说明文档打开之后不清楚每一步的物理背景更不清楚参数为什么那么设。这份资源的核心价值是把链路完整串了起来——从报文解析到输出带质量标记的AOD/WV时间序列适合刚接手CE318数据、或者正在搭建处理流程的人参考。2. 数据基础与定标电压读数换算成物理量前的三个关键步骤2.1 CE318的原始数据长什么样通道、列格式与读取CE318的原始输出是二进制格式厂商配套软件可以导出为文本文件每行记录包含时间戳和9个通道的电压值。通道波长通常是1020nm、870nm、670nm、500nm、440nm、380nm、340nm外加一个水汽通道936nm部分型号是940nm。电压值单位一般是毫伏或DN数不同型号导出的列顺序略有差别。读取的第一件事是确认列名和波长一一对应否则后面算AOD时通道错位结果完全没有可解释性。import pandas as pd import numpy as np # CE318文本导出格式示例date time V1020 V870 V670 V500 V440 V380 V340 V936 df pd.read_csv(ce318_20240915.txt, sepr\s, parse_dates{dt: [0, 1]}) df df.sort_values(dt).reset_index(dropTrue) # 把电压列重命名方便后续按波长调用 df.rename(columns{ V1020: v_1020, V870: v_870, V670: v_670, V500: v_500, V440: v_440, V380: v_380, V340: v_340, V936: v_936 }, inplaceTrue)读进来之后要做三件事时间戳统一成UTC、电压值检查负值或零值、按连续观测段做分组。CE318是自动跟踪太阳的雨天或太阳被遮挡时会输出异常低的电压甚至负数这些在读取阶段就要标记出来不能直接送进定标或反演流程。我一般会在读取后加一行筛掉电压小于阈值的记录阈值取晴天正常值的五分之一。这个阈值不能定得太高否则日出日落时段的大气质量数较大、电压自然偏低会被误删。2.2 定标的两条路线Langley法与标准灯野外观测首选前者定标是整条链路上最容易被忽略的一步。CE318测量的物理量本质是太阳直射辐照度对应的电压把电压换算成大气光学厚度需要知道大气顶的电压值V0这个值就是定标常数。两条路线一是标准灯定标在实验室用溯源到计量基准的标准灯照射仪器得到V0二是Langley定标在晴空稳定天气下利用不同太阳天顶角下电压随大气质量数变化的规律做线性拟合外推得到V0。野外观测站的首选几乎都是Langley法。原因很直接标准灯定标需要把仪器拆下来送回实验室运输途中光学窗口清洁状态的变化会引入额外误差而且实验室环境和野外环境温度不同V0会有可观的系统差。Langley定标只要碰上半个上午或半个下午的稳定晴空就能做在观测现场算完直接更新V0对站点运维来说省掉一大圈物流和时间成本。# Langley定标核心步骤对每个通道做 ln(V) ~ m 线性拟合 # m 为大气质量数用 kasten_young 函数计算 def langley_calibration(m, voltage): mask (m 4) (m 8) np.isfinite(voltage) coeff np.polyfit(m[mask], np.log(voltage[mask]), 1) v0 np.exp(coeff[1]) tau_total -coeff[0] return v0, tau_total, mask拟合时质量数下限取4是为了让光程足够长气溶胶的微小变化对拟合斜率影响相对小上限取8是因为天顶角太大时大气折射和仪器指向误差开始变得不可忽略。拟合出来的V0要和上一次定标结果做对比变化超过3%要警惕光学窗口脏污或探测器温漂超过5%建议重新做一次标准灯检定。这个判断阈值不是拍脑袋定死的而是来自V0不确定性对AOD误差的传递V0偏差1%AOD大约偏差0.01对清洁大气来说已经很明显。2.3 日地距离修正与大气质量数一个系数一个角度错了全盘皆输Beer-Lambert定律里的电压要修正到平均日地距离处否则季节性的日地距离变化会给AOD带来约±3.4%的虚假波动。日地距离修正因子d用年积日计算公式虽然简单但年积日必须取自UTC时间不能拿地方时直接算。大气质量数用Kasten-Young公式比简单的余弦倒数更准尤其在太阳天顶角超过70度时两者差异明显。def day_fraction(doy): 日地距离修正因子doy为UTC年积日 return 1 - 0.01672 * np.cos(np.deg2rad(0.9856 * (doy - 4))) def kasten_young(sza_deg): Kasten-Young 大气质量数公式输入太阳天顶角度 sza np.deg2rad(sza_deg) return 1.0 / (np.cos(sza) 0.50572 * (96.07995 - sza_deg) ** (-1.6364))调用时注意两个参数sza_deg传入的是0到90度的浮点数数组不能是弧度值doy取UTC年积日避免地方时和UTC跨日导致修正因子跳变。这两步做完电压数据才真正变成可以参与反演的物理量序列。很多初学脚本只做其中一个修正或者把d乘了一次全年AOD都带着明显的季节周期排查的时候很难想到是这个原因。3. AOD反演逐通道扣除瑞利散射与气体吸收剩下的才是气溶胶3.1 Beer-Lambert定律分解总光学厚度怎么拆反演AOD的起点是V0、日地距离修正d和大气质量数m都确定之后的Beer-Lambert定律。对每个通道电压V和总光学厚度τ_total的关系是ln(V / V0) -m·τ_total于是τ_total -ln(V / (V0·d²)) / m。τ_total包含了气溶胶消光、瑞利散射、臭氧吸收、NO2吸收等多个成分AOD的定义是气溶胶的消光光学厚度所以必须把其他成分逐个扣除。def total_optical_depth(v, v0, d, m): 由电压反算总光学厚度 return -np.log(v / (v0 * d * d)) / m这里有一个容易踩的细节d必须乘两次因为太阳辐照度随日地距离平方反比变化。很多初始版本只乘一次d结果AOD全年带周期性偏差在冬季和夏季的反演值系统性偏高或偏低。常见做法是先对全部通道统一算total_optical_depth再按通道波长分别扣除瑞利和气体吸收。通道之间的电压列长度必须一致否则数组对位错位影响是隐性的画时间序列才看得出来。3.2 瑞利散射光学厚度参数化公式与波长单位瑞利散射是分子散射可以用标准大气模型参数化计算。常用形式是τ_R (P / P0) · a · λ^(-b)其中P是站点气压P0是标准大气压1013.25hPaa和b的取值随公式版本略有差别。一套实际处理中常用的参数是a0.008735b4.08。这个公式对λ的单位极敏感λ必须用微米如果直接填纳米进去λ^(-4.08)这个量级会完全失控。def rayleigh_tau(wavelength_um, pressure_hpa1013.25): 瑞利散射光学厚度参数化波长单位微米 lam np.asarray(wavelength_um, dtypefloat) return 0.008735 * lam ** (-4.08) * pressure_hpa / 1013.25站点气压可以用观测站的气压计读数。如果没有实时气压用标准大气压代替会带来约5%的误差对440nm来说相当于0.005左右的AOD偏差在清洁大气场景下不能忽略。我一般会在配置文件里放一个气压序列文件路径脚本按日读取并插值到观测时刻。如果连气压序列也没有至少把站点海拔填进去用标准大气剖面算一个平均气压比直接套标准大气压要好。3.3 臭氧与NO2吸收什么时候该扣、什么时候可忽略臭氧在Chappuis带约450-750nm有明显吸收对500nm和670nm通道的AOD有实际影响。CE318的通道里没有臭氧专用测量波段所以臭氧光学厚度要么用外部臭氧柱总量数据要么用气候平均值。NO2在城市污染场景下对440nm有一点点贡献但相对臭氧和瑞利散射小一个量级日常处理可以先忽略做高精度科研反演时再补。# 臭氧吸收光学厚度简化计算o3_du 为臭氧柱总量Dobson Unit def o3_tau(wavelength_um, o3_du): 臭氧Chappuis带吸收系数查表获得这里给典型值 coef {0.500: 0.030, 0.675: 0.045} # 单位DU^-1 if wavelength_um in coef: return o3_du * coef[wavelength_um] return 0.0这段代码是简化版实际使用时要按仪器通道波长插值吸收系数表不能只靠两个离散点。臭氧柱总量数据源可以选同站点的探空资料或者卫星再分析产品把数值按日期填进配置文件脚本按日读取即可。需要强调的是臭氧吸收对380nm和340nm通道的贡献很小主要影响落在500-670nm区间如果只反演短波通道臭氧项可以不参与计算节省一个数据源的维护成本。3.4 云筛选与AOD输出三倍标准差法的参数选择与质量标记云筛选是决定AOD数据质量的关键一环。常用的三倍标准差法以某时刻为中心取一个时间窗口计算窗口内AOD均值与标准差当前值偏离均值超过三倍标准差就标记为云污染。窗口长度选择直接影响判别效果窗口太长会把真实的气溶胶变化误判成云太短则样本量不足统计量本身不稳定。def cloud_filter(aod_series, window_min10, n_sigma3.0): 三倍标准差云筛选返回布尔掩码True表示通过 radius pd.Timedelta(minuteswindow_min // 2) mask np.zeros(len(aod_series), dtypebool) for i, ts in enumerate(aod_series.index): win aod_series[(aod_series.index ts - radius) (aod_series.index ts radius)] if len(win) 3: mean, std win.mean(), win.std() mask[i] np.abs(aod_series.iloc[i] - mean) n_sigma * std return mask窗口我一般取10分钟最少5分钟。低于5分钟在日出日落时段样本太少均值本身不稳定超过20分钟会把短时气溶胶事件平滑掉。筛选后建议再叠加相邻点突变检测相邻两个时间点的AOD差值超过0.02后一点直接标记为可疑。到这一步AOD时间序列和质量标记就可以写进输出CSV了。输出列建议包含UTC时间、太阳天顶角、各通道AOD、质量标记质量标记按0/1/2分三级0为晴空可靠1为可疑2为云污染。4. 水汽WV反演936nm通道的改良Langley拟合与误差控制4.1 为什么用水汽通道而不是直接做差分水汽反演依赖936nm通道对水汽吸收的敏感响应。水汽的转动-振动吸收带分布在900-1000nm之间936nm处于一个吸收强度适中、不容易饱和的区间透过率和水汽总量的关系近似满足指数律。如果选940nm仪器带宽和滤光片特性会影响吸收系数的标定936nm窄带通道的标定相对稳定这是很多处理流程把它作为默认水汽通道的原因。水汽反演的目标量是整层大气的水汽柱总量WV单位习惯用cm或g/cm²数值在温带地区通常介于0.1到5之间。4.2 Modified Langley法的数学形式非线性项的线性化处理水汽通道的消光关系里光学厚度不再和水汽总量W成线性关系而是近似正比于(m·W)^bb在0.5左右。直接套用普通Langley拟合会得到弯曲的ln(V)-m曲线所以要用Modified Langley法把非线性项当作拟合参数一起解算。实际操作中直接对电压模型做非线性最小二乘拟合拟合时段选取晴空稳定、且WV变化平缓的连续观测段。from scipy.optimize import curve_fit def water_vapor_model(m, wv_cm, v0_936, k): 936nm水汽通道透过率模型V V0 * exp(-k * (m * wv)^0.5) return v0_936 * np.exp(-k * np.sqrt(m * wv_cm)) # m_sel 为晴空段大气质量数数组v_sel 为对应936nm电压 popt, pcov curve_fit(water_vapor_model, m_sel, v_sel, p0[0.5, 1e4, 0.3], maxfev5000) wv_cm, v0_fit, k_fit poptp0的初值设置很关键wv_cm给0.5是温带地区的典型柱水汽量v0给1e4是前面Langley定标得到的V0数量级k给0.3是水汽弱吸收带的典型吸收系数。初值偏差太大会导致拟合不收敛或者收敛到局部极小值反演出的WV明显不合理。拟合完成后要看两个指标拟合残差是否在平稳区间内以及反演出的v0_fit和主通道定标得到的V0_936是否在同一量级。如果v0_fit偏离V0_936超过5%说明拟合时段选的不好或云筛不彻底。4.3 水汽反演的误差来源温度漂移、V0不确定性与时段选择936nm通道的探测器响应受温度影响明显仪器虽然有温控但野外温差大的环境下温控精度有限电压读数会叠加温度漂移。处理办法有两个方向一是用仪器记录的温度做线性校正把电压归一到一个参考温度二是直接挑选温度变化平缓的时段做拟合避开昼夜温差大的过渡段。第二个方向在站点数据里更容易落实因为拟合本身只需要一段连续晴空不是全天数据都要参与。注意V0_936的定标不确定度对WV影响很大经验上V0有1%的误差WV会偏约2%。所以水汽通道的定标常数建议每季度与主通道一起复查不要只做一次就长期使用。时段选择上我习惯优先用上午的观测段。上午大气通常更稳定气溶胶和WV的演变比下午平缓Modified Langley拟合的残差更小。下午对流发展起来后水汽垂直分布变化快拟合质量明显变差。如果数据里只有下午可用至少要确认拟合时段内WV没有明显趋势否则可以把时段拆成两段分别拟合并对比结果。5. 避坑CE318数据处理中我踩过的五个具体的坑5.1 定标系数过期未复查AOD整体系统性偏高现象同一时段反演出的AOD比同站另外两台仪器稳定高0.03-0.05波长越短越明显。原因仪器光学窗口在雨季被灰尘污染V0实际值比使用值低等于拿偏大的V0去反演导致AOD偏高。解决每个月用Langley法复查V0若连续两次变化超过3%就更新定标文件同时检查窗口清洁状态。从那以后我把定标检查脚本直接挂进了月度批处理流程不跑完定标质量检查就不准出数据。5.2 云筛选窗口设成30分钟晴空数据混入云污染点现象AOD时间序列午后频繁出现尖峰尖峰持续几分钟就消失而周边点看起来正常。原因窗口过长均值把云的瞬时贡献平均掉了当前点对比一个被抬高的均值时没能被识别出来。解决窗口改成10分钟并对每个点叠加相邻点差值检测差值超过0.02直接标记为可疑。宁可牺牲少量晴空点也要保住数据可靠性AOD的平滑度比点数更重要。5.3 936nm缺测用线性插值补WV时间序列出现中凹陷现象WV曲线正午出现一个明显下凹但周边时段正常而水汽在晴空下本应平稳变化。原因936nm通道在正午有一小段缺测线性插值引入了几分钟的水汽“假变化”后续Modified Langley拟合对这个假变化非常敏感。解决缺测点不超过5分钟就直接剔除超过5分钟用半小时中值代替插值后重新检查WV曲线是否平滑。这里最稳妥的办法是把936nm通道的缺测处理独立出来不要跟其他通道共用插值逻辑。5.4 瑞利散射公式里的波长单位填错短波AOD翻车现象440nm和380nm的AOD全部偏大0.2以上870nm和1020nm正常。第一个反应是定标有问题重做Langley定标后依旧如此。原因瑞利散射公式里λ用纳米代入λ^(-4.08)在短波部分被严重放大短波通道偏差远大于长波通道。解决在脚本入口加一个波长范围断言如果输入值大于1.5就判定为纳米并除以1000从源头堵住这个坑。这个断言我后来一直在用虽然简单但确实拦住了不少粗心操作。5.5 大气质量数用简单余弦倒数日出日落时段数据全废现象太阳天顶角大于75度时AOD系统性偏大越接近地平线越夸张。原因余弦倒数公式在低太阳高度角时没有考虑大气折射和地球曲率按平直大气算出的m偏大乘到总光学厚度上就放大了。解决改用Kasten-Young公式并且SZA超过80度的数据直接不参与反演和定标。日出日落那几分钟的数据看似多实际上误差大不值得保留。6. 结果验证与批量处理波长指数自检和脚本化收尾6.1 AOD波长指数α作为质量自检指标反演出的AOD能不能用除了跟参考数据对比之外还有一个不依赖外部数据的快速自检手段——Ångström波长指数α。用440nm和870nm两个通道的AOD计算α -ln(AOD440/AOD870) / ln(440/870)这个值能反映气溶胶粒径分布城市污染型气溶胶α在1.0-1.5之间沙尘粒子主导时α降到0.3以下海盐粒子大致在0.5左右。处理完一段数据后画出α的时间序列突然跳到负数或者急剧上升多半是某个通道的定标或云筛出了问题值得回查。def angstrom_alpha(aod_440, aod_870): 440-870nm 波长指数用于气溶胶类型和质量检验 return -np.log(aod_440 / aod_870) / np.log(440 / 870)使用这个自检指标时先过滤掉AOD小于0.005的数据点避免log运算在极低值处不稳定。α本身没有绝对阈值关键是看它是否在合理气溶胶类型范围内以及时间序列是否平滑。如果α在半小时内从1.3跳到0.4再跳回来几乎可以肯定是数据质量问题。6.2 批量处理脚本把一个月的数据自动跑完并归档实际观测经常是按月攒数据逐天手工跑不现实。批量脚本的骨架是用glob把目录下所有观测文件收集起来循环执行读取、定标检查、反演、云筛和质量标记最后把每天的AOD和WV分别写进按日期命名的CSV文件里。每条记录都要带质量标记字段后续画图时用颜色区分不同等级的数据点。for f in ce318_2024*.txt; do python process_ce318.py --input $f --site 某观测站 --config config.yaml done资源附带的处理脚本就是按这个流程组织的配置项单独放在config.yaml里站点纬度、气压、定标文件路径、臭氧柱总量数据路径、云筛窗口和阈值全部集中一处修改。这个安排的好处是换站点或换年份时只需要修改配置文件不用动处理代码。脚本跑完后自动生成一个summary.csv汇总每天的有效数据时长、云筛剔除比例和WV均值方便快速判断这个月的数据质量。前面说过那一次漏更新定标系数导致整个季度的AOD全部系统性偏低返工花了两天。从那以后我每次跑批量前都强制走一遍定标质量检查脚本对比V0序列之后再决定是否继续。希望这份资源里的处理流程和避坑记录能帮你少走这段弯路。本文还有配套的精品资源点击获取 CE318太阳光度计数据处理：从电压到AOD与水汽WV反演全链路解析
简介在大气气溶胶与温室气体研究中AOD气溶胶光学厚度和WV水汽含量是反映大气环境状态的重要参数CE318太阳光度计通过多波段太阳辐射观测为反演这些参数提供原始数据。一套聚焦CE318数据处理的C源码资源可用于实现从观测数据到AOD、WV反演的完整流程主要面向大气科学、遥感应用及环境监测方向的研究生、科研人员与相关课程学习者。压缩包共5个文件均为C源码大小仅9KB分别覆盖原始观测数据读取、太阳几何位置解算、Angström指数拟合、气溶胶光学厚度计算和水汽反演等功能模块代码紧凑、结构清晰按处理流程拆分既可对单环节单独调用也能串联成完整反演链。目前已有1082人学习下载。借助该资源读者不仅能从源码层面复现CE318观测数据的处理流程还能理解辐射传输、气溶胶散射吸收以及水汽吸收带选择等核心算法并结合自身站点数据调整参数后投入使用为科研实验、毕业设计或教学演示提供一份轻量且可扩展的参考实现。1. 从CE318原始电压到AOD和水汽WV数据处理链路到底卡在哪一个观测站放着一台CE318太阳光度计采集了一整年的原始电压数据报文里每行是十几个通道的电压读数看起来规整但真要把它变成逐分钟的气溶胶光学厚度AOD和水汽柱含量WV中间隔着一长串步骤定标、日地距离修正、云筛选、瑞利散射扣除、臭氧吸收订正最后还要对936nm水汽通道做非线性拟合。不少做气溶胶观测的从业者拿到一套CE318处理脚本和说明文档打开之后不清楚每一步的物理背景更不清楚参数为什么那么设。这份资源的核心价值是把链路完整串了起来——从报文解析到输出带质量标记的AOD/WV时间序列适合刚接手CE318数据、或者正在搭建处理流程的人参考。2. 数据基础与定标电压读数换算成物理量前的三个关键步骤2.1 CE318的原始数据长什么样通道、列格式与读取CE318的原始输出是二进制格式厂商配套软件可以导出为文本文件每行记录包含时间戳和9个通道的电压值。通道波长通常是1020nm、870nm、670nm、500nm、440nm、380nm、340nm外加一个水汽通道936nm部分型号是940nm。电压值单位一般是毫伏或DN数不同型号导出的列顺序略有差别。读取的第一件事是确认列名和波长一一对应否则后面算AOD时通道错位结果完全没有可解释性。import pandas as pd import numpy as np # CE318文本导出格式示例date time V1020 V870 V670 V500 V440 V380 V340 V936 df pd.read_csv(ce318_20240915.txt, sepr\s, parse_dates{dt: [0, 1]}) df df.sort_values(dt).reset_index(dropTrue) # 把电压列重命名方便后续按波长调用 df.rename(columns{ V1020: v_1020, V870: v_870, V670: v_670, V500: v_500, V440: v_440, V380: v_380, V340: v_340, V936: v_936 }, inplaceTrue)读进来之后要做三件事时间戳统一成UTC、电压值检查负值或零值、按连续观测段做分组。CE318是自动跟踪太阳的雨天或太阳被遮挡时会输出异常低的电压甚至负数这些在读取阶段就要标记出来不能直接送进定标或反演流程。我一般会在读取后加一行筛掉电压小于阈值的记录阈值取晴天正常值的五分之一。这个阈值不能定得太高否则日出日落时段的大气质量数较大、电压自然偏低会被误删。2.2 定标的两条路线Langley法与标准灯野外观测首选前者定标是整条链路上最容易被忽略的一步。CE318测量的物理量本质是太阳直射辐照度对应的电压把电压换算成大气光学厚度需要知道大气顶的电压值V0这个值就是定标常数。两条路线一是标准灯定标在实验室用溯源到计量基准的标准灯照射仪器得到V0二是Langley定标在晴空稳定天气下利用不同太阳天顶角下电压随大气质量数变化的规律做线性拟合外推得到V0。野外观测站的首选几乎都是Langley法。原因很直接标准灯定标需要把仪器拆下来送回实验室运输途中光学窗口清洁状态的变化会引入额外误差而且实验室环境和野外环境温度不同V0会有可观的系统差。Langley定标只要碰上半个上午或半个下午的稳定晴空就能做在观测现场算完直接更新V0对站点运维来说省掉一大圈物流和时间成本。# Langley定标核心步骤对每个通道做 ln(V) ~ m 线性拟合 # m 为大气质量数用 kasten_young 函数计算 def langley_calibration(m, voltage): mask (m 4) (m 8) np.isfinite(voltage) coeff np.polyfit(m[mask], np.log(voltage[mask]), 1) v0 np.exp(coeff[1]) tau_total -coeff[0] return v0, tau_total, mask拟合时质量数下限取4是为了让光程足够长气溶胶的微小变化对拟合斜率影响相对小上限取8是因为天顶角太大时大气折射和仪器指向误差开始变得不可忽略。拟合出来的V0要和上一次定标结果做对比变化超过3%要警惕光学窗口脏污或探测器温漂超过5%建议重新做一次标准灯检定。这个判断阈值不是拍脑袋定死的而是来自V0不确定性对AOD误差的传递V0偏差1%AOD大约偏差0.01对清洁大气来说已经很明显。2.3 日地距离修正与大气质量数一个系数一个角度错了全盘皆输Beer-Lambert定律里的电压要修正到平均日地距离处否则季节性的日地距离变化会给AOD带来约±3.4%的虚假波动。日地距离修正因子d用年积日计算公式虽然简单但年积日必须取自UTC时间不能拿地方时直接算。大气质量数用Kasten-Young公式比简单的余弦倒数更准尤其在太阳天顶角超过70度时两者差异明显。def day_fraction(doy): 日地距离修正因子doy为UTC年积日 return 1 - 0.01672 * np.cos(np.deg2rad(0.9856 * (doy - 4))) def kasten_young(sza_deg): Kasten-Young 大气质量数公式输入太阳天顶角度 sza np.deg2rad(sza_deg) return 1.0 / (np.cos(sza) 0.50572 * (96.07995 - sza_deg) ** (-1.6364))调用时注意两个参数sza_deg传入的是0到90度的浮点数数组不能是弧度值doy取UTC年积日避免地方时和UTC跨日导致修正因子跳变。这两步做完电压数据才真正变成可以参与反演的物理量序列。很多初学脚本只做其中一个修正或者把d乘了一次全年AOD都带着明显的季节周期排查的时候很难想到是这个原因。3. AOD反演逐通道扣除瑞利散射与气体吸收剩下的才是气溶胶3.1 Beer-Lambert定律分解总光学厚度怎么拆反演AOD的起点是V0、日地距离修正d和大气质量数m都确定之后的Beer-Lambert定律。对每个通道电压V和总光学厚度τ_total的关系是ln(V / V0) -m·τ_total于是τ_total -ln(V / (V0·d²)) / m。τ_total包含了气溶胶消光、瑞利散射、臭氧吸收、NO2吸收等多个成分AOD的定义是气溶胶的消光光学厚度所以必须把其他成分逐个扣除。def total_optical_depth(v, v0, d, m): 由电压反算总光学厚度 return -np.log(v / (v0 * d * d)) / m这里有一个容易踩的细节d必须乘两次因为太阳辐照度随日地距离平方反比变化。很多初始版本只乘一次d结果AOD全年带周期性偏差在冬季和夏季的反演值系统性偏高或偏低。常见做法是先对全部通道统一算total_optical_depth再按通道波长分别扣除瑞利和气体吸收。通道之间的电压列长度必须一致否则数组对位错位影响是隐性的画时间序列才看得出来。3.2 瑞利散射光学厚度参数化公式与波长单位瑞利散射是分子散射可以用标准大气模型参数化计算。常用形式是τ_R (P / P0) · a · λ^(-b)其中P是站点气压P0是标准大气压1013.25hPaa和b的取值随公式版本略有差别。一套实际处理中常用的参数是a0.008735b4.08。这个公式对λ的单位极敏感λ必须用微米如果直接填纳米进去λ^(-4.08)这个量级会完全失控。def rayleigh_tau(wavelength_um, pressure_hpa1013.25): 瑞利散射光学厚度参数化波长单位微米 lam np.asarray(wavelength_um, dtypefloat) return 0.008735 * lam ** (-4.08) * pressure_hpa / 1013.25站点气压可以用观测站的气压计读数。如果没有实时气压用标准大气压代替会带来约5%的误差对440nm来说相当于0.005左右的AOD偏差在清洁大气场景下不能忽略。我一般会在配置文件里放一个气压序列文件路径脚本按日读取并插值到观测时刻。如果连气压序列也没有至少把站点海拔填进去用标准大气剖面算一个平均气压比直接套标准大气压要好。3.3 臭氧与NO2吸收什么时候该扣、什么时候可忽略臭氧在Chappuis带约450-750nm有明显吸收对500nm和670nm通道的AOD有实际影响。CE318的通道里没有臭氧专用测量波段所以臭氧光学厚度要么用外部臭氧柱总量数据要么用气候平均值。NO2在城市污染场景下对440nm有一点点贡献但相对臭氧和瑞利散射小一个量级日常处理可以先忽略做高精度科研反演时再补。# 臭氧吸收光学厚度简化计算o3_du 为臭氧柱总量Dobson Unit def o3_tau(wavelength_um, o3_du): 臭氧Chappuis带吸收系数查表获得这里给典型值 coef {0.500: 0.030, 0.675: 0.045} # 单位DU^-1 if wavelength_um in coef: return o3_du * coef[wavelength_um] return 0.0这段代码是简化版实际使用时要按仪器通道波长插值吸收系数表不能只靠两个离散点。臭氧柱总量数据源可以选同站点的探空资料或者卫星再分析产品把数值按日期填进配置文件脚本按日读取即可。需要强调的是臭氧吸收对380nm和340nm通道的贡献很小主要影响落在500-670nm区间如果只反演短波通道臭氧项可以不参与计算节省一个数据源的维护成本。3.4 云筛选与AOD输出三倍标准差法的参数选择与质量标记云筛选是决定AOD数据质量的关键一环。常用的三倍标准差法以某时刻为中心取一个时间窗口计算窗口内AOD均值与标准差当前值偏离均值超过三倍标准差就标记为云污染。窗口长度选择直接影响判别效果窗口太长会把真实的气溶胶变化误判成云太短则样本量不足统计量本身不稳定。def cloud_filter(aod_series, window_min10, n_sigma3.0): 三倍标准差云筛选返回布尔掩码True表示通过 radius pd.Timedelta(minuteswindow_min // 2) mask np.zeros(len(aod_series), dtypebool) for i, ts in enumerate(aod_series.index): win aod_series[(aod_series.index ts - radius) (aod_series.index ts radius)] if len(win) 3: mean, std win.mean(), win.std() mask[i] np.abs(aod_series.iloc[i] - mean) n_sigma * std return mask窗口我一般取10分钟最少5分钟。低于5分钟在日出日落时段样本太少均值本身不稳定超过20分钟会把短时气溶胶事件平滑掉。筛选后建议再叠加相邻点突变检测相邻两个时间点的AOD差值超过0.02后一点直接标记为可疑。到这一步AOD时间序列和质量标记就可以写进输出CSV了。输出列建议包含UTC时间、太阳天顶角、各通道AOD、质量标记质量标记按0/1/2分三级0为晴空可靠1为可疑2为云污染。4. 水汽WV反演936nm通道的改良Langley拟合与误差控制4.1 为什么用水汽通道而不是直接做差分水汽反演依赖936nm通道对水汽吸收的敏感响应。水汽的转动-振动吸收带分布在900-1000nm之间936nm处于一个吸收强度适中、不容易饱和的区间透过率和水汽总量的关系近似满足指数律。如果选940nm仪器带宽和滤光片特性会影响吸收系数的标定936nm窄带通道的标定相对稳定这是很多处理流程把它作为默认水汽通道的原因。水汽反演的目标量是整层大气的水汽柱总量WV单位习惯用cm或g/cm²数值在温带地区通常介于0.1到5之间。4.2 Modified Langley法的数学形式非线性项的线性化处理水汽通道的消光关系里光学厚度不再和水汽总量W成线性关系而是近似正比于(m·W)^bb在0.5左右。直接套用普通Langley拟合会得到弯曲的ln(V)-m曲线所以要用Modified Langley法把非线性项当作拟合参数一起解算。实际操作中直接对电压模型做非线性最小二乘拟合拟合时段选取晴空稳定、且WV变化平缓的连续观测段。from scipy.optimize import curve_fit def water_vapor_model(m, wv_cm, v0_936, k): 936nm水汽通道透过率模型V V0 * exp(-k * (m * wv)^0.5) return v0_936 * np.exp(-k * np.sqrt(m * wv_cm)) # m_sel 为晴空段大气质量数数组v_sel 为对应936nm电压 popt, pcov curve_fit(water_vapor_model, m_sel, v_sel, p0[0.5, 1e4, 0.3], maxfev5000) wv_cm, v0_fit, k_fit poptp0的初值设置很关键wv_cm给0.5是温带地区的典型柱水汽量v0给1e4是前面Langley定标得到的V0数量级k给0.3是水汽弱吸收带的典型吸收系数。初值偏差太大会导致拟合不收敛或者收敛到局部极小值反演出的WV明显不合理。拟合完成后要看两个指标拟合残差是否在平稳区间内以及反演出的v0_fit和主通道定标得到的V0_936是否在同一量级。如果v0_fit偏离V0_936超过5%说明拟合时段选的不好或云筛不彻底。4.3 水汽反演的误差来源温度漂移、V0不确定性与时段选择936nm通道的探测器响应受温度影响明显仪器虽然有温控但野外温差大的环境下温控精度有限电压读数会叠加温度漂移。处理办法有两个方向一是用仪器记录的温度做线性校正把电压归一到一个参考温度二是直接挑选温度变化平缓的时段做拟合避开昼夜温差大的过渡段。第二个方向在站点数据里更容易落实因为拟合本身只需要一段连续晴空不是全天数据都要参与。注意V0_936的定标不确定度对WV影响很大经验上V0有1%的误差WV会偏约2%。所以水汽通道的定标常数建议每季度与主通道一起复查不要只做一次就长期使用。时段选择上我习惯优先用上午的观测段。上午大气通常更稳定气溶胶和WV的演变比下午平缓Modified Langley拟合的残差更小。下午对流发展起来后水汽垂直分布变化快拟合质量明显变差。如果数据里只有下午可用至少要确认拟合时段内WV没有明显趋势否则可以把时段拆成两段分别拟合并对比结果。5. 避坑CE318数据处理中我踩过的五个具体的坑5.1 定标系数过期未复查AOD整体系统性偏高现象同一时段反演出的AOD比同站另外两台仪器稳定高0.03-0.05波长越短越明显。原因仪器光学窗口在雨季被灰尘污染V0实际值比使用值低等于拿偏大的V0去反演导致AOD偏高。解决每个月用Langley法复查V0若连续两次变化超过3%就更新定标文件同时检查窗口清洁状态。从那以后我把定标检查脚本直接挂进了月度批处理流程不跑完定标质量检查就不准出数据。5.2 云筛选窗口设成30分钟晴空数据混入云污染点现象AOD时间序列午后频繁出现尖峰尖峰持续几分钟就消失而周边点看起来正常。原因窗口过长均值把云的瞬时贡献平均掉了当前点对比一个被抬高的均值时没能被识别出来。解决窗口改成10分钟并对每个点叠加相邻点差值检测差值超过0.02直接标记为可疑。宁可牺牲少量晴空点也要保住数据可靠性AOD的平滑度比点数更重要。5.3 936nm缺测用线性插值补WV时间序列出现中凹陷现象WV曲线正午出现一个明显下凹但周边时段正常而水汽在晴空下本应平稳变化。原因936nm通道在正午有一小段缺测线性插值引入了几分钟的水汽“假变化”后续Modified Langley拟合对这个假变化非常敏感。解决缺测点不超过5分钟就直接剔除超过5分钟用半小时中值代替插值后重新检查WV曲线是否平滑。这里最稳妥的办法是把936nm通道的缺测处理独立出来不要跟其他通道共用插值逻辑。5.4 瑞利散射公式里的波长单位填错短波AOD翻车现象440nm和380nm的AOD全部偏大0.2以上870nm和1020nm正常。第一个反应是定标有问题重做Langley定标后依旧如此。原因瑞利散射公式里λ用纳米代入λ^(-4.08)在短波部分被严重放大短波通道偏差远大于长波通道。解决在脚本入口加一个波长范围断言如果输入值大于1.5就判定为纳米并除以1000从源头堵住这个坑。这个断言我后来一直在用虽然简单但确实拦住了不少粗心操作。5.5 大气质量数用简单余弦倒数日出日落时段数据全废现象太阳天顶角大于75度时AOD系统性偏大越接近地平线越夸张。原因余弦倒数公式在低太阳高度角时没有考虑大气折射和地球曲率按平直大气算出的m偏大乘到总光学厚度上就放大了。解决改用Kasten-Young公式并且SZA超过80度的数据直接不参与反演和定标。日出日落那几分钟的数据看似多实际上误差大不值得保留。6. 结果验证与批量处理波长指数自检和脚本化收尾6.1 AOD波长指数α作为质量自检指标反演出的AOD能不能用除了跟参考数据对比之外还有一个不依赖外部数据的快速自检手段——Ångström波长指数α。用440nm和870nm两个通道的AOD计算α -ln(AOD440/AOD870) / ln(440/870)这个值能反映气溶胶粒径分布城市污染型气溶胶α在1.0-1.5之间沙尘粒子主导时α降到0.3以下海盐粒子大致在0.5左右。处理完一段数据后画出α的时间序列突然跳到负数或者急剧上升多半是某个通道的定标或云筛出了问题值得回查。def angstrom_alpha(aod_440, aod_870): 440-870nm 波长指数用于气溶胶类型和质量检验 return -np.log(aod_440 / aod_870) / np.log(440 / 870)使用这个自检指标时先过滤掉AOD小于0.005的数据点避免log运算在极低值处不稳定。α本身没有绝对阈值关键是看它是否在合理气溶胶类型范围内以及时间序列是否平滑。如果α在半小时内从1.3跳到0.4再跳回来几乎可以肯定是数据质量问题。6.2 批量处理脚本把一个月的数据自动跑完并归档实际观测经常是按月攒数据逐天手工跑不现实。批量脚本的骨架是用glob把目录下所有观测文件收集起来循环执行读取、定标检查、反演、云筛和质量标记最后把每天的AOD和WV分别写进按日期命名的CSV文件里。每条记录都要带质量标记字段后续画图时用颜色区分不同等级的数据点。for f in ce318_2024*.txt; do python process_ce318.py --input $f --site 某观测站 --config config.yaml done资源附带的处理脚本就是按这个流程组织的配置项单独放在config.yaml里站点纬度、气压、定标文件路径、臭氧柱总量数据路径、云筛窗口和阈值全部集中一处修改。这个安排的好处是换站点或换年份时只需要修改配置文件不用动处理代码。脚本跑完后自动生成一个summary.csv汇总每天的有效数据时长、云筛剔除比例和WV均值方便快速判断这个月的数据质量。前面说过那一次漏更新定标系数导致整个季度的AOD全部系统性偏低返工花了两天。从那以后我每次跑批量前都强制走一遍定标质量检查脚本对比V0序列之后再决定是否继续。希望这份资源里的处理流程和避坑记录能帮你少走这段弯路。本文还有配套的精品资源点击获取