简介面向大气科学、遥感与环保监测领域的学习者和研究者这份资源针对CE318太阳光度计观测数据提供从原始数据读取到AOD气溶胶光学厚度与水汽含量WV反演的实现方案。压缩包共5个文件均为C源码整体仅9KB按数据读取、仪器定标、角度拟合、AOD计算和水汽计算等环节拆分为独立模块覆盖太阳几何订正与大气衰减修正等关键处理步骤便于按需调用与对照理解。已有1082人学习适合具备大气光学基础和一定C编程能力的读者。资源通过源码直接展示完整反演链路从原始数据解析与清洗开始经定标和角度拟合最终得到AOD与水汽含量结果各模块可独立修改或组合复用适合在此基础上搭建自己的批量处理流程也可作为学习CE318数据反演方法的参考资料。1. 从CE318的DN值到AOD一条绕不开的数据链路手上攒了一整年的CE318太阳光度计原始文件想反演AOD和气柱水汽含量却发现网上的教程要么只讲原理要么只给个AERONET下载入口。真正从.dat数据走到逐小时的AOD时间序列中间还隔着定标系数、云筛选、气压订正、臭氧订正这四道关卡。这篇笔记把我实际跑通CE318数据反演AOD和936nm水汽含量的链路拆开讲重点放在可以直接复用的Python流程和参数边界上顺便把容易翻车的几个点都指出来。适合大气遥感、环境监测方向的研究生和一线数据工程师手里有CE318原始数据但不想把时间全耗在摸索上的读者可以直接照着走。2. CE318数据预处理文件结构、定标系数与云筛选三关2.1 CE318原始文件里到底有什么CE318是法国CIMEL公司的多波段太阳光度计AERONET全球网络的标准仪器。它测的不是AOD本身而是太阳直射辐射经过大气衰减后的仪器响应值也就是DN值。常用配置有8个通道不同通道对应的科学用途很有必要先记清楚后面反演时要一一对应。通道号中心波长(nm)主要用途1340臭氧和气溶胶边界信息2380细粒子气溶胶参考3440AOD标准通道4500AOD标准通道5675AOD标准通道6870AOD标准通道水汽反演参考7936水汽通道81020粗粒子气溶胶原始文件的行格式因固件版本不同差异很大常见的有两种一种是纯ASCII文本每行包含时间戳、8组DN值、太阳天顶角、方位角、仪器温度、环境气压、湿度另一种是类似NASA Ames格式的带列头的表格。拿到文件的第一件事不是看AOD而是把每一列对应到物理量。我一般会在读取脚本里把列索引列成一个诊断表格逐个通道打印前五行DN值确认没有错位。这个环节最容易出问题的场景是老仪器升级后多了一个1640nm通道导致后续所有通道的索引整体往后挪了一位。如果按原顺序去读440nm通道算出来的AOD会变成负值而且很难排查。读取代码建议先把时间解析和DN值提取做成两个独立函数因为后续的云筛选、Langley定标、AOD反演都要反复调用它们。时间解析尤其要注意有的固件输出UTC有的输出本地时如果后续要和AERONET Level 2.0对比统一转成UTC是底线。2.2 Langley定标为什么反演之前必须先求V0AOD反演的本质是比尔-朗伯定律的变形太阳辐射穿过大气后仪器测到的DN值与大气总光学厚度呈指数关系。假设V0是仪器在大气层顶测到的DN值那么对某一波长有V V0 × R² × exp(-m × τ_total)其中R是日地距离修正因子m是大气质量数τ_total是大气总光学厚度。把这个式子两边取对数再对m做线性回归截距就是lnV0斜率就是总光学厚度。这个过程叫Langley定标。Langley定标听起来简单坑却很多。首先它要求大气在观测时段内保持稳定所以一般选清晨或傍晚太阳天顶角较大、气团数在2.5到6之间的连续时段。其次单日数据回归的R²至少要大于0.995才可信。如果某一天云层不稳定回归点的离散度会明显变大。实际操作时我会连续选三天晴朗清晨的数据分别做回归把三天得到的V0取平均。不要只依赖单日结果因为即使肉眼看到是晴天高层薄云也可能让截距偏移几个百分点。对CE318的非水汽通道这个流程是通用的但936nm水汽通道不能直接这样算原因在第四章细说。下面是V0标定的参考代码输入是一天的原始数据输出该通道的V0和拟合质量。这段代码也适合用来快速筛查哪些天的数据适合做定标。import numpy as np import pandas as pd def langley_calibration(df, channel_col, sza_col, press_colNone): 对单个通道做Langley回归 df: 包含DN值、太阳天顶角的DataFrame channel_col: 目标通道DN值列名 sza_col: 太阳天顶角列名 press_col: 气压列名可选用于订正瑞利光学厚度 # 只取清晨或傍晚大气质量数在2.5~6之间 # 太阳天顶角越大大气质量数越大 sza df[sza_col].values mask (sza 60) (sza 85) # 粗略对应 m2.5~6 if mask.sum() 10: return None # 计算大气质量数Kasten Young公式 cosz np.cos(np.radians(sza[mask])) m_air 1.0 / (cosz 0.50572 * (96.07995 - sza[mask]) ** -1.6364) # 日地距离修正DOY为年积日 # R 1 - 0.0167 * cos(2*pi*(DOY-3)/365)V/R^2后再取对数 doy df[doy].values[mask] R 1 - 0.0167 * np.cos(2 * np.pi * (doy - 3) / 365.0) y np.log(df[channel_col].values[mask] / R ** 2) # 线性拟合x是大气质量数y是ln(V/R^2) coeffs np.polyfit(m_air, y, 1) slope, intercept coeffs y_fit np.polyval(coeffs, m_air) resid y - y_fit ss_res np.sum(resid ** 2) ss_tot np.sum((y - np.mean(y)) ** 2) r2 1 - ss_res / ss_tot return { v0: np.exp(intercept), slope: -slope, # 斜率对应总光学厚度 r2: r2, n_points: int(mask.sum()) }这段代码的关键参数是天顶角范围。m2.5对应天顶角约66.4°m6对应约84.1°。我实际用了60到85度这个窗口能覆盖大部分有效定标时段。R²低于0.995的结果我一般不采信宁可多等一个晴朗清晨。如果现场有气压传感器还可以把大气质量数换成“气压订正后的大气质量数”让回归更稳但对V0截距的影响不大日常定标可以暂时忽略。需要强调一点V0会随时间漂移尤其是滤光片受潮或仪器运输后。CE318的现场维护手册通常建议每三个月做一次Langley定标如果站点条件不允许至少也要每半年一次。用一年前甚至三年前的V0反演AOD误差会直接体现在AOD系统偏高或偏低。2.3 云筛选把数据里的云影先踢出去云是CE318反演最大的干扰源。云不像气溶胶那样缓慢变化它会短时间内让DN值剧烈波动如果不筛选反演出的AOD会偏大且跳动剧烈。AERONET把云筛选放在Level 1.5的质控里常规做法是用440nm、675nm、870nm三个通道的信号稳定性来判断。我的做法分两步。第一步是单点筛选检查相邻两次观测之间的DN值跳变幅度。CE318默认的观测间隔在晴天模式约15分钟在特殊观测模式约1分钟。更常用的基线是约2到3分钟一次的“太阳直射天空漫射”配合观测模式。def cloud_screen(df, ch_list, jump_limit0.0015, window3): 基于通道DN值稳定性和时间连续性的云筛选 ch_list: 需要检查的通道列名列表 jump_limit: 相邻两次观测DN值的最大相对变化率 window: 滑动窗口内变异系数的阈值判断 df df.sort_values(time).reset_index(dropTrue) keep np.ones(len(df), dtypebool) # 检查相邻跳变 for ch in ch_list: v df[ch].values.astype(float) rel_change np.abs(np.diff(v)) / np.maximum(v[:-1], 1e-6) # 如果变化率超过0.15%后一个点视为云影 jump_flag np.zeros(len(v), dtypebool) jump_flag[1:] rel_change jump_limit keep ~jump_flag # 检查滑动窗口内变异系数窗口大小5 for ch in ch_list: v df[ch].values.astype(float) # 用滚动窗口算变异系数超过阈值则标记 for i in range(2, len(df) - 2): window_vals v[i-2:i3] cv np.std(window_vals) / np.mean(window_vals) if cv 0.01: keep[i] False return df[keep]这段代码的思路是先对每个通道做相邻跳变检测再做滑动窗口的变异系数检测。跳变阈值0.15%是我在多数站点调出来的经验值晴天稳定大气下几分钟内DN值变化一般低于0.1%一旦有云影扫过变化率会瞬间到1%以上。窗口变异系数阈值0.01则用来捕捉那种“变化不大但一直抖”的薄云情况。筛选完成之后建议再叠加一个太阳天顶角限制。天顶角大于84度时大气质量数超过6订正误差被放大反演出的AOD噪声很大直接截掉更省心。AERONET Level 2.0实际上还包含一个“平滑度检验”会进一步剔除那些虽然瞬时稳定但在更长时间尺度上波动的数据。我们自己做预处理时做到上面的两步已经能满足大部分科研和业务需求。3. 朗利定标到AOD反演一份可复用的Python实现3.1 AOD反演公式的完整展开把比尔-朗伯定律整理成AOD的显式表达式所有反演流程都从这里出发。总光学厚度由三部分组成大气分子瑞利散射、气溶胶消光、痕量气体吸收。因此气溶胶光学厚度τ_a等于总光学厚度减去瑞利散射τ_R和臭氧、二氧化氮等气体吸收项。实际工程计算中每一项都有对应的处理方式。瑞利散射光学厚度必须做气压订正因为分子散射强度与大气柱分子数成正比。臭氧吸收在440nm以下不可忽略在500nm以上影响较小如果站点没有臭氧柱总量数据我会用0.1 atm-cm的经验值先顶着并在结果说明里标注这一点。NO₂吸收在可见光波段影响小日常反演多数站点直接忽略。输出的AOD还有个自检方法多波段AOD之间应该满足Ångström指数关系。如果440nm和870nm的AOD比值明显违反这个谱分布规律要么是定标有问题要么是云筛选没做干净。3.2 逐通道AOD反演代码与参数设置下面是完整的AOD反演函数输入是经过云筛选的DataFrame和V0字典输出是逐通道AOD和Ångström指数。这段代码我把它当模板用换站点时只需要改V0字典和气压来源。import numpy as np import pandas as pd # 波长表单位微米与CE318通道顺序对应 WAVE [0.340, 0.380, 0.440, 0.500, 0.675, 0.870, 0.936, 1.020] def rayleigh_od(wave_um, pressure_hpa): Bodhaine等1999年的瑞利光学厚度公式 wave_um: 波长单位微米 pressure_hpa: 站点气压单位hPa p0 1013.25 x wave_um tau_r 0.00864 * x ** (-3.916 0.074 * x 0.05 / x) * (pressure_hpa / p0) return tau_r def airmass(sza_deg): Kasten Young大气质量数 cosz np.cos(np.radians(sza_deg)) return 1.0 / (cosz 0.50572 * (96.07995 - sza_deg) ** -1.6364) def aod_from_raw(df, v0_dict, pressure_hpa1013.25, ozone_cm0.1): 从CE318原始DN值反演逐通道AOD df: 云筛选后的数据至少包含时间、天顶角、各通道DN列 v0_dict: {波长微米: V0值} pressure_hpa: 站点气压若无气压传感器可传默认值 ozone_cm: 臭氧柱总量单位atm-cm doy df[doy].values sza df[sza].values R 1 - 0.0167 * np.cos(2 * np.pi * (doy - 3) / 365.0) m airmass(sza) ozone_abs_coeff { 0.340: 0.0, 0.380: 0.0, 0.440: 0.02, 0.500: 0.03, 0.675: 0.01, 0.870: 0.0, 0.936: 0.0, 1.020: 0.0 } # 瑞利光学厚度 tau_r_dict {w: rayleigh_od(w, pressure_hpa) for w in WAVE} results {} for w in WAVE: # 假设列名是 dn_ 波长纳米例如 dn_440 col fdn_{int(w * 1000)} v df[col].values.astype(float) v0 v0_dict[w] ln_term np.log(v0 / (v * R ** 2)) tau_total ln_term / m tau_r tau_r_dict[w] tau_o3 ozone_abs_coeff[w] * ozone_cm tau_a tau_total - tau_r - tau_o3 results[w] tau_a aod_df pd.DataFrame(results) aod_df[sza] sza aod_df[time] df[time].values # 计算Ångström指数用440/870两个标准通道 alpha -np.log(aod_df[0.440] / aod_df[0.870]) / np.log(0.440 / 0.870) aod_df[angstrom_440_870] alpha return aod_df关键参数要交代清楚臭氧吸收系数我用的单位是cm⁻¹乘上臭氧柱总量atm-cm后是无量纲吸收光学厚度。臭氧通道340nm和380nm在CE318数据里常被用来辅助判断对流层臭氧但反演AOD时如果臭氧数据不准确宁可把440nm以下的AOD结果打上“仅供参考”的标记。气压参数pressure_hpa如果站点有气压传感器优先用实时值如果没有用标准大气压会导致瑞利订正在高海拔站点明显偏高AOD偏低0.01到0.03不等。日地距离修正因子R在代码里每年初给一个平均值即可但如果是长时间序列建议逐日计算。上面代码里的cos项是最常见的高精度近似相比查表法误差小于0.01%。3.3 反演结果的合理性自检AOD反演出来之后不要急着画图先做三个自检。第一440nm AOD是否明显大于870nm正常情况下城市和大陆站点440nm AOD约是870nm的1.5到2倍如果出现倒挂检查V0是不是某个通道写反了。第二全部波长AOD是否大于零出现负值通常意味着V0偏小、气压订正过度或者云筛选漏掉了云影响。第三晴空条件下AOD的时间序列应该平滑变化任何相邻时次超过0.05的跳变都要回看原始DN值。我一般会把反演结果和同站点的能见度或PM2.5观测做定性对比。能见度差、PM2.5高发的时段AOD应该明显抬升。如果AOD与这些外部数据完全不相关问题大概率出在V0上。4. 936nm水汽含量反演改良朗利法与查表法的Python实现4.1 水汽通道反演的特殊性936nm位于水汽的近红外弱吸收带CE318用这个通道和870nm参考通道的透过率比值来反演大气柱水汽含量PWV。870nm基本不受水汽吸收影响936nm的额外衰减主要由水汽造成。因此把936nm总透过率扣除瑞利散射和气溶胶消光后剩下的就是水汽透过率Tw。水汽透过率和PWV之间的关系不满足简单的指数衰减常用半经验模型Tw exp(-a × (m × W)^b)其中W是大气柱水汽含量单位cm或g/cm²m是大气质量数a和b是经验参数。反解得到W ((-ln Tw) / a)^(1/b) / m问题是936nm通道的V0不能直接用标准Langley回归获取因为水汽在大气中的时空变化太剧烈总光学厚度并不稳定。常用的替代方案有两种一种叫“修正Langley法”选取水汽含量相对稳定的时段把936nm通道的信号做多元回归同时拟合出V0另一种是通过与探空或GNSS水汽数据对比标定一组a和b参数然后查表计算。我实际跑业务数据时更倾向第二种因为修正Langley法对数据时段的要求太苛刻。下面给出一个两段式的实现先改进Langley法求V0_936再用水汽透过率模型反演PWV。4.2 改良朗利法求936通道V0改良的思路是把936nm的透过率拆成非水汽部分和水汽部分。先用标准方法得到870nm的V0反演出870nm AOD然后通过Ångström指数外推936nm的气溶胶光学厚度。把936nm的总光学厚度减掉瑞利、气溶胶和臭氧剩下的就是水汽吸收的光学厚度。这时对“水汽透过率”取对数理论上可以表示为-ln Tw a × (m × W)^b在一个较短的时间窗口内如果水汽柱总量W近似不变那么-ln Tw对m的b次方做回归截距应该接近0但V0需要另行求解。更常见的实际操作是反过来把V0_936当作待定参数对一组选定的晴空数据调整V0_936使Tw满足0到1且平滑变化。这样做虽然带有一点“调参”的意味但只要选点合理结果比强行做线性回归稳定得多。下面这段代码实现了改良朗利法的主流程。def calibrate_v0_936(df, aod_870, alpha, v0_870, pressure_hpa): 反推936通道V0的改良Langley流程 df: 包含dn_936、dn_870、sza的DataFrame aod_870: 由870nm通道反演的AOD序列 alpha: 440-870nm的Ångström指数序列 v0_870: 870nm通道V0 # 936nm气溶胶光学厚度通过Ångström外推 aod_936 aod_870 * (0.936 / 0.870) ** (-alpha) # 936nm瑞利光学厚度 tau_r_936 rayleigh_od(0.936, pressure_hpa) # 大气质量数 m airmass(df[sza].values) # 日地距离修正 R 1 - 0.0167 * np.cos(2 * np.pi * (df[doy].values - 3) / 365.0) # 暂时假设V0_936后续用搜索方式确定 candidate_v0 np.linspace(0.8, 1.2, 201) * v0_870 best_std np.inf best_v0 None for v0_936 in candidate_v0: v936 df[dn_936].values # 水汽透过率 ln_tw np.log(v0_936 / (v936 * R ** 2)) / m - tau_r_936 - aod_936 tw np.exp(-ln_tw) # 物理上需要0tw1且标准差最小 if np.any(tw 0) or np.any(tw 1): continue tw_std np.std(tw) if tw_std best_std: best_std tw_std best_v0 v0_936 return best_v0, best_std这段代码用“标准差最小化”作为搜索V0_936的判据本质上是假设晴空水汽透过率在短时间内变化不大。候选搜索区间设在0.8到1.2倍的870nm V0之间这是因为同台仪器的936通道和870通道光电响应差异不会超过这个范围。如果最佳V0落在搜索边界就把区间扩大重跑一次。要注意的是候选V0要乘以R²吗这里我把日地距离订正放在分母位置的v×R²上v0_936本身也已经经过同样的R²换算因此两边保持一致。经验上成功标定后Tw的标准差应小于0.02如果标准差大于0.05基本可以判断这段数据里有云或气溶胶突变不要用来做定标。4.3 水汽含量反演与参数选择拿到V0_936之后水汽反演就变成直接计算了。def pwv_from_tw(df, v0_936, a, b): 从936nm水汽透过率反演大气柱水汽含量 a, b: 水汽透过率模型参数 sza df[sza].values m airmass(sza) doy df[doy].values R 1 - 0.0167 * np.cos(2 * np.pi * (doy - 3) / 365.0) v936 df[dn_936].values.astype(float) # 总透过率未扣除气溶胶和瑞利 t_total v936 * R ** 2 / v0_936 # 此处需要外部传入水汽透过率实际使用时和AOD外推一起算 ln_tw np.log(v936 * R ** 2 / v0_936) / m # 简化处理ln_tw已经包含了瑞利和气溶胶订正下面是完整版 tw np.exp(-ln_tw) pwv ((-np.log(tw)) / a) ** (1 / b) / m return pwv参数a和b的选取是这个环节最玄学的地方。公开文献中936nm通道常用的经验值在a0.08、b0.55附近但不同仪器的滤光片带宽不同直接套用会带来系统性偏差。最稳妥的做法是收集站点附近探空站的PWV数据和同期CE318反演的-lnTw做回归拟合出本地化的a和b。如果暂时没有探空数据GNSS气象观测的PWV产品也可以作为替代基准。反演结果必须看量级中纬度夏季PWV通常2到4cm冬季0.5到1.5cm热带站点全年5cm左右。如果算出来20cm甚至更大基本是V0_936偏小或参数b不合适。水汽通道数据在太阳天顶角大于70度时信噪比急剧下降最好提前截断。5. CE318反演避坑指南五个高频翻车点与排查流程5.1 现象AOD早上偏高、中午偏低系统性日变化明显原因通常是日地距离修正没有做或者修正因子用反了。CE318原始DN值是仪器直接输出的数字量不少固件版本并不会自动除以日地距离R²需要自己在反演时补上。如果没有做这一步夏季反演的AOD会偏高约0.02冬季偏低约0.02。解决方法是把日地距离修正写进反演代码的公共函数里确保Langley定标和逐日AOD反演用的是同一套R²。另一个常见来源是Langley定标数据选在了冬夏跨季V0的截距把季节性的R²差异吸收进了定标系数里导致另一个季节反演的系统偏差。5.2 现象某个通道AOD出现大量负值负AOD最常见的原因是V0不匹配。CE318每个通道的V0差异很大尤其是340nm和380nm短波通道对定标误差非常敏感。实际遇到过用户把V0文件里的列顺序和通道顺序搞反导致440nm和500nm两个通道对调反演出的AOD在440nm变成负值。解决方法是构建V0字典时明确键值为波长而不是通道序号。每个通道反演前先打印V0和使用波长人工核对一次。剩余偶发负值则检查云筛选是否到位薄云导致的DN值轻微偏低会把AOD压低。5.3 现象936nm水汽反演结果跳变时大时小现象是PWV时间序列里出现孤立尖峰和相邻观测差好几倍。原因多半是936nm通道信号在低太阳高度角时太弱信噪比不足水汽透过率被噪声抬高。解决方法是把天顶角大于70度的数据直接剔除同时增加一个PWV物理范围的硬过滤例如0.1到7cm之间。还有一点容易忽略936nm通道对仪器温度比较敏感如果原始文件里有仪器温度列建议对936nm的DN值做温度订正后再反演。订正系数各仪器不同可以通过对比同一时段不同温度下的DN值回归出来。5.4 现象云筛选后数据量锐减但仍有时段有云残留数据量一下砍掉60%而且剩下序列里还有云污染信号这是因为云筛选只用了一个通道的瞬时跳变没有利用多通道的一致性。薄云和高云对短波和长波通道的影响程度不同440nm和870nm的AOD比值在云影下会异常偏大。解决方法是加入Ångström指数阈值判断。晴空条件下440到870nm的Ångström指数一般在0.5到1.8之间云影响下这个值常常冲到2.5以上。把这个条件叠加到原来的DN值稳定判据上能很大程度减少云残留。5.5 现象和AERONET官方Level 2.0相差0.02以上且呈系统性偏差AERONET官方数据的定标和质控流程比手工流程严格很多如果差异是系统性而不是随机噪声首先要对比V0。AERONET站点会定期公布仪器校准系数下载他们的V0文件替换自己的值通常就能把差值缩小到0.01以内。另一个原因是没有做NO₂订正。污染站点NO₂柱总量高时500nm以下通道反演的AOD会偏高。解决方式是在反演代码里增加一个NO₂吸收订正项系数取经验值即可。6. 用AERONET交叉验证结果十分钟确认AOD和水汽是否可信交叉验证不需要多复杂把AERONET同站点同时段的Level 2.0数据下载下来和自己反演的结果做时间对齐和偏差统计就能快速定位问题。具体步骤先到AERONET官网选择站点下载Level 2.0的AOD和PWV数据文件时间跨度覆盖自己数据所在的日期。然后按时间戳最近邻匹配计算两个指标平均偏差MBE和均方根误差RMSE。指标计算公式可接受范围MBEmean(反演值 - 参考值)440nm: ±0.01以内RMSEsqrt(mean(差值²))440nm: 0.02以内R²线性回归决定系数0.95以上匹配时要注意时区一致建议全部转成UTC。如果MBE稳定在正0.02说明自己的V0偏低如果只在高AOD时偏差扩大多半是Langley定标数据里混入了薄云。我自己的习惯是先做一次十分钟快检再决定要不要把前一天的批量数据重新跑一遍。AERONET数据只取Level 2.0不要拿Level 1.5来验证因为Level 1.5的云筛选并不严格两者差异反而容易误导排查方向。另外建议把反演的时间序列和同一站点的太阳辐射计或GNSS水汽做并排图。即使数值有微小偏差只要趋势一致就可以放心用本地定标参数做业务输出。从那以后我每次批量反演之前都会强制走一遍“Langley回归画图、V0人工确认、云筛选参数复检”的流程花不了几分钟却能避开大半夜的返工。希望这篇笔记能帮你在CE318数据处理上少踩几个同样的坑。本文还有配套的精品资源点击获取