简介面向大气遥感与气溶胶观测的学习者这份CE318太阳光度计数据处理资料围绕AOD与大气水汽含量wv的反演展开适合具备一定大气光学基础、希望用C实现从原始观测到参数计算的读者。压缩包共5个文件均为cpp源码整体约9KB按功能覆盖CE318数据读取、气溶胶光学厚度与水汽含量的对比计算、仪器定标以及角度拟合等环节可作为课程实验或科研预处理的参考实现。已有1084人学习下载说明其在同类小众数据处理资料中具有一定参考价值。读者可从中了解观测数据清洗、辐射量比较与反演结果校核的组织方式并借助源码梳理AOD与wv计算的模块划分、关键参数传递和常见异常处理思路便于移植到Python或其他分析流程中。包体轻量、主题集中更适合作为算法原型和代码框架来查阅而非完整业务系统。1. 从一台 CE318 的观测文件说起AOD 和 wv 到底怎么落地凌晨三点自动观测站的 CE318 又传回一包 .dat 文件里面是八个波段、每几分钟一次的太阳直射电压。对做气溶胶和大气水汽的人来说这些数字真正的价值只有两个气溶胶光学厚度 AOD以及整层水汽含量 wv。CE318 的太阳光度计数据看起来门槛不高无非是电压除以定标常数、取对数、再除以大气质量数但真要把 AOD 做到能和公开基准产品对上、把 940 nm 通道的水汽反演稳定在可发布水平中间隔着定标、云筛选、瑞利订正、臭氧订正和查找表系数这几道坎。这篇笔记按我实际处理观测文件的顺序拆先讲原始数据和 triplet 质量控制再写 AOD 反演和 Angstrom 指数然后落到 940 nm 水汽通道最后把踩过的坑和交叉验证方法摊开。适合手头有 CE318 观测文件、准备自己写处理脚本的同行也适合刚接手站点数据、需要判断现有流程靠不靠谱的人。2. 数据预处理从原始电压到可用的观测量2.1 先认清 CE318 的测量序列和文件结构CE318 的一次观测不是单帧而是一组按时间排列的测量序列。典型文件里包含日期、时间、太阳天顶角、各波段直射电压以及仪器内部温度、气压等辅助量。不同固件导出的列名和分隔符不一样常见的是空格或逗号分隔的文本字段顺序也随版本变化。我一般先做一件事拿一天的文件把头十行和尾十行打印出来确认时间格式、波段顺序、有没有表头、有没有以#开头的注释行。这一步不做后面列名对不上AOD 会整列错位而且错得很隐蔽。处理之前要把三个辅助量准备好观测点的气压、臭氧柱总量、日地距离修正因子。气压影响瑞利光学厚度臭氧影响短波波段日地距离影响定标常数的换算。CE318 文件里通常有气压字段但如果没有就要从站点气象数据补。臭氧可以用当天的卫星或再分析资料精度不需要极高但 440 nm 附近不能忽略。日地距离因子按儒略日计算常见做法是用简化的椭圆轨道公式。这些量不需要手工算写成函数挂在读取流程后面即可。另一个容易被忽略的点是时间。CE318 记录的是 UTC 还是本地时间不同站点配置不同。AOD 反演本身对时间不敏感但一旦要和公开基准产品逐点对比时区错一小时就会让匹配失败。我的习惯是读入后先统一转成 UTC并在数据框里保留原始时间列方便回溯。时间戳解析用pandas.to_datetime时显式指定utcTrue不要依赖默认本地时区。2.2 云筛选和 triplet 质量控制CE318 的 AOD 反演里云是最大的污染源。薄云会让电压偏高AOD 偏低厚云直接让信号饱和或掉零。单帧判断云很难所以厂家和公开处理流程都用 triplet 方法连续三次测量作为一个 triplet三个值之间的变异系数超过阈值就整组丢弃。常见阈值是 1% 到 5%具体看波段和站点气候。我一般先按 2% 跑一遍再看剩余数据量如果剩余太少再放宽到 5%但要在日志里注明。triplet 的实现不复杂按时间排序每三条一组计算每个波段电压的标准差除以均值。注意不是对 AOD 做而是对原始电压做因为电压更接近观测值本身。还要处理一种情况三次测量中间缺了一次不能硬凑直接跳过这组。另一个细节是太阳天顶角太低的太阳高度角对应太大的人气质量数AOD 误差会被放大通常把天顶角大于 75° 到 80° 的数据剔掉。这个阈值不要一刀切低纬度站点可以放到 80°高纬度冬季可能只能用到 70°。除了 triplet还要做电压范围检查。每个波段的定标常数 V0 是已知的测量电压如果大于 V0 的 1.2 倍基本可以判定异常小于 0 或等于 0 的直接丢。这个检查能挡住大部分仪器饱和和遮光问题。做完这两步剩下的数据才进入 AOD 计算。2.3 读取与清洗的 Python 实现下面这段代码假设你已经把原始文件整理成 CSV列名统一为datetime, airmass, v440, v675, v870, v940, v1020。实际文件列名不同改rename映射即可。代码里保留了原始时间列并做了 triplet 标记。import numpy as np import pandas as pd # 读取已导出的 CSV原始文件如果是空格分隔把 sep 改成 r\s df pd.read_csv(ce318_day.csv, parse_dates[datetime]) df[datetime] pd.to_datetime(df[datetime], utcTrue) df df.sort_values(datetime).reset_index(dropTrue) # 基本电压检查负值、零值、超过定标常数 1.2 倍的标记为无效 V0 {v440: 12500.0, v675: 9800.0, v870: 7600.0, v940: 6200.0, v1020: 5400.0} for col, v0 in V0.items(): df[f{col}_valid] (df[col] 0) (df[col] 1.2 * v0) # triplet 质量控制每三条一组计算每个波段的变异系数 df[triplet_id] np.arange(len(df)) // 3 triplet_stats df.groupby(triplet_id).agg( n(datetime, count), cv440(v440, lambda x: x.std() / x.mean() if x.mean() 0 else np.nan), cv675(v675, lambda x: x.std() / x.mean() if x.mean() 0 else np.nan), cv870(v870, lambda x: x.std() / x.mean() if x.mean() 0 else np.nan), cv940(v940, lambda x: x.std() / x.mean() if x.mean() 0 else np.nan), ) # 只有三条齐全且四个波段变异系数都小于 0.02 的 triplet 保留 good_triplet triplet_stats[ (triplet_stats[n] 3) (triplet_stats[[cv440, cv675, cv870, cv940]] 0.02).all(axis1) ].index df[triplet_ok] df[triplet_id].isin(good_triplet) # 太阳天顶角筛选这里用 airmass 反推airmass 4 约对应天顶角 75° 以上 df[geometry_ok] df[airmass] 4.0 df_clean df[df[triplet_ok] df[geometry_ok]].copy() print(f原始 {len(df)} 条清洗后 {len(df_clean)} 条)这段代码的逻辑分三层。第一层是电压合法性V0是示意值实际必须换成你仪器定标文件里的数不能直接抄。第二层是 triplet按行号每三条一组要求三条齐全且各波段变异系数小于 0.02cv计算的是标准差除以均值均值接近零时返回NaN避免除零。第三层是几何筛选airmass 4.0是经验阈值对应太阳天顶角大约 75°低纬度站点可以放到 4.5高纬度冬季建议收到 3.0。清洗后不要急着算 AOD先把每个波段剩余数据量打印出来如果 440 nm 只剩个位数说明当天云太多或者阈值太严这份数据不适合做日平均。3. AOD 反演Beer-Lambert、定标和 Angstrom 指数3.1 从电压到光学厚度的公式与定标常数CE318 测的是太阳直射电压满足 Beer-Lambert 定律的离散形式V(λ) V0(λ) * (d0/d)^2 * exp(-m * τ(λ))。其中V0是定标常数d0/d是日地距离修正m是大气质量数τ是总光学厚度。反解出来就是τ ln(V0 * (d0/d)^2 / V) / m。总光学厚度里包含气溶胶、瑞利散射、臭氧吸收、二氧化氮吸收水汽吸收在 940 nm 单独处理不放进 AOD 波段。要求 AOD就把非气溶胶项减掉τ_aer τ_total - τ_Ray - τ_O3 - τ_NO2。定标常数V0是整个流程的命门。它通常来自厂家定标文件也可能来自站点自己做的 Langley 定标。不同仪器、不同波段、不同定标日期V0都不一样。我见过最常见的翻车是脚本里硬编码了一组V0换了一台仪器还在用AOD 整体偏移 0.05 到 0.1肉眼看曲线还挺正常和基准产品一比就露馅。正确做法是把V0放在单独的配置文件里文件名带仪器序列号和定标日期脚本启动时校验序列号是否匹配。Langley 定标不是本文重点但要知道如果V0是自己拟合的必须用晴空稳定日、大气质量数跨度足够大的数据且拟合残差要检查。日地距离修正因子按观测时刻的儒略日算常见近似公式是(d0/d)^2 1.00011 0.034221*cos(g) 0.00128*sin(g) 0.000719*cos(2g) 0.000077*sin(2g)其中g 2π*(day_of_year-1)/365.25。这个因子在 1 月附近约 1.0347 月附近约 0.967量级 3%对 AOD 的影响在小光学厚度时不能忽略。很多脚本图省事直接取 1短波波段会带来 0.01 左右的系统偏差。3.2 气溶胶光学厚度计算与波长插值AOD 计算按波段逐个做。先算总光学厚度再减瑞利、臭氧、二氧化氮。瑞利光学厚度用近似公式τ_R 0.008735 * λ^(-4.08) * P/1013.25λ单位是微米P是站点气压。这个公式在 440 到 1020 nm 范围内够用短波 340 nm 附近误差会大一些如果做 340 nm 要换更精确的 Bodhaine 参数化。臭氧光学厚度用τ_O3 k_O3(λ) * ozone_DU / 1000吸收系数k_O3随波长查表440 nm 附近约 0.0003 每 DU 量级675 nm 以上基本可以忽略。二氧化氮在污染站点才需要认真订正清洁站点可以设零但要在文档里写明。下面是一个按波段计算 AOD 的函数定标常数和系数都从字典传入不硬编码在公式里。import numpy as np # 臭氧吸收系数示意单位每 DU 的光学厚度实际用查表值替换 O3_COEF {440: 0.00032, 675: 0.00004, 870: 0.00001, 1020: 0.00000} def rayleigh_od(wavelength_nm, pressure_hpa): 瑞利光学厚度近似波长单位 nm气压单位 hPa lam_um wavelength_nm / 1000.0 return 0.008735 * lam_um ** (-4.08) * pressure_hpa / 1013.25 def no2_od(wavelength_nm, no2_du): 二氧化氮光学厚度清洁站点可传 0 # 简化系数实际按波段查表 coef {440: 0.00012, 675: 0.00002, 870: 0.00000, 1020: 0.00000} return coef.get(wavelength_nm, 0.0) * no2_du def compute_aod(row, v0, pressure_hpa, ozone_du, no2_du0.0): 输入一行观测返回各波段 AOD 字典 out {} for wl in [440, 675, 870, 1020]: v row[fv{wl}] if not np.isfinite(v) or v 0: out[wl] np.nan continue # 日地距离修正已在 row 中预先算好列名 earth_sun_factor tau_total np.log(v0[fv{wl}] * row[earth_sun_factor] / v) / row[airmass] tau_aer (tau_total - rayleigh_od(wl, pressure_hpa) - O3_COEF.get(wl, 0.0) * ozone_du - no2_od(wl, no2_du)) out[wl] tau_aer return out这段代码里row[earth_sun_factor]是提前算好的(d0/d)^2airmass是大气质量数。v0字典必须来自当前仪器的定标文件pressure_hpa用观测时刻的站点气压ozone_du用当天臭氧柱总量。O3_COEF和no2_od里的系数只是示意实际要换成你所用波段宽度和滤光片对应的查表值。函数返回NaN而不是抛异常是为了让批量处理能继续跑但事后要统计NaN比例。如果 440 nm 的 AOD 出现负值且绝对值大于 0.05先别急着当噪声丢检查定标常数和气压是不是用错了。3.3 Angstrom 指数与波长依赖的校验AOD 本身是波长相关的常用 Angstrom 指数描述这种依赖τ_aer(λ) β * λ^(-α)α就是 Angstrom 指数。用 440 nm 和 870 nm 两个波段算α -ln(τ440/τ870) / ln(440/870)。α大于 1 通常表示细模态气溶胶占优小于 0.5 表示粗模态占优沙尘或海盐场景常见。这个指数不参与水汽反演但它是检查 AOD 质量的好工具如果某天α突然从 1.2 掉到 0.1或者出现负值大概率是某个波段的定标或云筛选出了问题。我一般会把α和τ440画成散点图正常站点会形成一个有物理意义的带状分布异常点会飘在外面。还有一个检查是波长插值用 440 和 870 的 AOD 去预测 675 的 AOD看残差是否在 0.02 以内。如果 675 nm 的残差系统性偏大可能是 675 波段的定标常数有问题或者臭氧订正用错了。这些检查不需要额外数据用已有的三个波段就能做。Angstrom 指数对水汽反演也有间接影响940 nm 通道的 AOD 不能直接用来算水汽但如果你要检查 940 nm 的信号是否合理可以用α外推一个 940 nm 的 AOD 参考值和 940 nm 的总光学厚度比看看水汽吸收项是不是大得离谱。这个技巧在后面水汽章节还会用到。4. 水汽含量 wv 反演940 nm 通道的查找表与经验系数4.1 940 nm 水汽通道的物理基础940 nm 落在水汽吸收带内CE318 这个波段的电压同时受气溶胶散射、瑞利散射和水汽吸收影响。水汽反演的核心思路是先假设 940 nm 的气溶胶光学厚度可以用邻近波段外推再扣掉瑞利和臭氧剩下的吸收光学厚度就归给水汽。实际处理里更常用的是经验公式或查找表因为 940 nm 的滤光片带宽、仪器响应函数和水汽吸收线型都会影响结果纯解析计算很难做准。常见做法有两种。第一种是厂家配套软件里的经验公式形式大致是wv (1/m) * (ln(V0_940/V_940) / a)有的还带二次项或气压修正。第二种是查找表用辐射传输模式预先算好不同水汽含量、不同天顶角、不同气溶胶条件下的大气透过率反演时查表插值。两种方法没有绝对优劣经验公式快、适合批量查找表准、但依赖模式输入。我一般先用经验公式跑一版再用查找表或公开基准产品抽检偏差超过 10% 就回头查系数。不管用哪种方法V0_940都必须来自同一次定标且不能和 AOD 波段的V0混用。940 nm 的定标常数对水汽反演的影响是系统性的V0偏 2%wv 可能偏 5% 以上因为公式里是取对数后再除以吸收系数误差会被放大。另一个关键是大气质量数水汽反演用的大气质量数和水汽吸收的质量数不完全一样经验公式里通常仍用几何大气质量数但低太阳高度角时误差明显所以水汽反演也要做天顶角筛选。4.2 查找表法和经验公式的实现下面给一个经验公式的批量实现包含二次项修正和异常值过滤。公式结构是wv (1/m) * (ln(V0/V) / a) b * (ln(V0/V))^2其中a和b是仪器和波段相关的系数必须从定标文件或厂家文档里取代码里的数值只是占位不能直接用于生产。import numpy as np def retrieve_wv(row, v0_940, a_coef, b_coef, max_wv8.0): 经验公式反演水汽含量。 row 需要包含v940, airmass, earth_sun_factor v0_940: 940 nm 定标常数 a_coef, b_coef: 经验系数来自定标文件 max_wv: 物理上限超过视为异常 v row[v940] m row[airmass] if not np.isfinite(v) or v 0 or m 0: return np.nan # 日地距离修正后的等效电压比 ratio v0_940 * row[earth_sun_factor] / v if ratio 1.0: # 电压高于定标常数可能是云或仪器异常 return np.nan ln_ratio np.log(ratio) # 经验公式一次项加二次项修正 wv (ln_ratio / a_coef b_coef * ln_ratio ** 2) / m if not np.isfinite(wv) or wv 0 or wv max_wv: return np.nan return wv代码里ratio 1.0的判断很关键正常情况下V0 * earth_sun_factor应该大于测量电压因为大气有消光如果反过来说明测量电压异常偏高常见原因是薄云增强了前向散射或者仪器刚做过定标但常数没更新。max_wv8.0是经验上限实际站点水汽很少超过 6 cm超过 8 基本可以判为异常。a_coef和b_coef的取值因仪器而异不要抄网上的通用值去找你仪器最近一次定标报告。如果找不到二次项系数可以令b_coef0先用线性项跑但要在结果里标注方法版本。批量跑的时候不要只输出 wv 一列。把ln_ratio、airmass、earth_sun_factor和原始电压一起保留后面排查异常时能直接看出是哪一项出的问题。我习惯再加一列wv_flag标记ratio1、wvmax_wv、airmass4等不同原因统计时按 flag 分组比只看NaN有用。4.3 结果对比与异常值处理水汽反演结果要和至少一种独立来源对比。常见的是公开观测网的 Level 2.0 产品或者同站点的微波辐射计、探空资料。对比时不能只比日均值要按时间窗口匹配比如 ±30 分钟内的均值否则日变化会被抹平。匹配后算偏差和均方根误差偏差在 10% 以内、均方根误差在 0.3 cm 以内一般可以接受如果偏差是系统性的优先查V0_940和吸收系数如果误差是随机的查云筛选和天顶角阈值。异常值处理要克制。水汽在时间序列上连续变化突然出现一个 0.2 cm 或 7 cm 的点大概率是云或仪器问题可以直接丢。但如果连续几个小时都偏低不要逐点丢先看是不是气压或臭氧订正用了常数、实际当天变化较大。还有一种情况是 940 nm 电压接近饱和ln_ratio很小反演出的 wv 接近零这种点会被wv0过滤掉但过滤比例高的时候要检查是不是仪器增益设置不对。5. 避坑与排查AOD 和 wv 反演里最容易翻车的五件事5.1 定标常数用错整条曲线偏移现象AOD 时间序列形状正常但和基准产品比整体偏高或偏低 0.05 以上Angstrom 指数也偏得不多。原因通常是V0用了别的仪器、别的定标日期或者日地距离因子没乘。解决把V0从脚本里挪到配置文件文件名带序列号和定标日期启动时校验日地距离因子单独写函数单元测试用一个已知日期对一遍。我现在的习惯是每次换仪器先跑一天数据和基准产品比日均值偏差超过 0.02 就不往下做。5.2 triplet 阈值一刀切晴空数据被误杀现象清洗后数据量骤减尤其是冬季或高纬度站点剩下不到原始量的 10%。原因是变异系数阈值 0.02 对低信号波段太严或者三条测量里有一条缺测被整组丢掉。解决按波段设不同阈值440 nm 可以严一点1020 nm 放宽到 0.05缺测的 triplet 不要硬凑但可以标记后单独评估。另一个办法是先按天统计剩余率连续多天低于 20% 就回头检查阈值不要直接认为天气差。5.3 水汽公式漏了二次项或气压修正现象wv 在低水汽时和基准产品对得上高水汽时系统性偏低偏差随 wv 增大而增大。原因是经验公式只用了线性项高水汽时吸收饱和线性近似不够。解决找定标报告里的二次项系数没有就改用查找表同时检查气压修正有些公式里的m需要乘P/1013.25。如果二次项系数来源不明宁可不加也不要用网上抄来的通用值否则低水汽段会被带偏。5.4 时间戳时区错位匹配全失败现象AOD 和 wv 自身看着正常但和基准产品按时间匹配后一条都对不上或者匹配上的都是错误时段。原因是原始文件用本地时间脚本按 UTC 解析差了几个小时。解决读入时保留原始时间字符串显式转换时区输出文件里同时写 UTC 和本地时间两列匹配前先打印前五条时间戳和基准产品对一眼。这个坑不复杂但一旦发生排查成本很高因为数据本身没错错的是索引。5.5 把云污染当晴空AOD 偏低 wv 偏高现象AOD 出现不合理的低值同时 wv 偏高Angstrom 指数跳变。原因是薄云增强了 940 nm 的前向散射电压偏高ln_ratio偏小wv 反而被算大AOD 波段电压偏高则让 AOD 偏低。解决除了 triplet再加一道“相邻波段一致性”检查用 440、675、870 的 AOD 拟合 Angstrom 指数残差超过 0.05 的整组标记水汽那边检查ratio是否接近 1接近 1 的直接丢。云筛选没有万能阈值只能多层过滤加人工抽检。6. 从单站到批量用公开基准产品做交叉验证和自动化单站跑通只是第一步真正要拿这份数据处理长期观测得把批量脚本和交叉验证固定下来。我的做法是每天定时跑前一天的文件输出三张表清洗后的观测表、AOD 表、wv 表每张表都带质量标记和版本号。版本号不是形式主义定标常数、系数、阈值任何一项改了都要升版本否则半年后回头看根本分不清哪批结果用的是什么参数。批量脚本里加一个config.yaml把仪器序列号、定标日期、V0字典、triplet 阈值、天顶角上限、水汽系数全写进去脚本启动时打印配置摘要日志里留一份。交叉验证我一般按周做不按天。按天做容易被个别云污染点带偏按周看偏差和趋势更稳。匹配窗口用 ±30 分钟先把基准产品的时间转成 UTC再按最近邻匹配匹配不上就跳过不要插值。对比指标至少看三个平均偏差、均方根误差、相关系数。AOD 的相关系数低于 0.95 就要查水汽的均方根误差大于 0.4 cm 就要查。查的顺序是先看时间匹配数量再看异常点分布最后才动公式。很多时候问题不在公式在匹配和筛选。下面是一个批量交叉验证的骨架重点在匹配和统计不在绘图。import pandas as pd import numpy as np def match_nearest(df_a, df_b, time_coldatetime, tolerance30min): 按最近时间匹配两个数据框返回匹配后的合并表 df_a df_a.sort_values(time_col).copy() df_b df_b.sort_values(time_col).copy() df_a[time_col] pd.to_datetime(df_a[time_col], utcTrue) df_b[time_col] pd.to_datetime(df_b[time_col], utcTrue) matched pd.merge_asof( df_a, df_b, ontime_col, directionnearest, tolerancepd.Timedelta(tolerance), suffixes(_mine, _ref) ) return matched.dropna(subset[AOD440_ref]) def summary_stats(matched, col_mine, col_ref): 计算偏差、均方根误差、相关系数 x matched[col_mine].astype(float) y matched[col_ref].astype(float) mask np.isfinite(x) np.isfinite(y) x, y x[mask], y[mask] bias (x - y).mean() rmse np.sqrt(((x - y) ** 2).mean()) corr np.corrcoef(x, y)[0, 1] if len(x) 2 else np.nan return {n: len(x), bias: bias, rmse: rmse, corr: corr}pd.merge_asof要求两个表都按时间排序directionnearest是最近邻tolerance30min控制最大允许间隔。匹配后一定要看n如果一周只剩几十个点先别分析偏差回去查时间格式和筛选阈值。summary_stats里先做有限值过滤避免NaN污染统计。这个骨架不绘图但输出可以直接写进日志或周报。我的习惯是每周一早上跑一遍上周的匹配把n、bias、rmse、corr四列贴在一张表里连续几周bias漂移就说明定标或系数该更新了。还有一个小技巧把 AOD 和 wv 的交叉验证结果按天气类型分组晴天、薄云、沙尘各看一遍。整体指标好看不代表每种场景都对沙尘天的 Angstrom 指数低水汽反演的云筛选阈值可能偏松分组之后问题更容易暴露。从那以后我每次换定标常数都强制走一遍“单日对比 → 周统计 → 分组检查”这三步不再只看一条时间序列顺眼就发布。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/27 6:13:51 · 阅读 1.2k CE318太阳光度计AOD与水汽反演：质控、定标和交叉验证
简介面向大气遥感与气溶胶观测的学习者这份CE318太阳光度计数据处理资料围绕AOD与大气水汽含量wv的反演展开适合具备一定大气光学基础、希望用C实现从原始观测到参数计算的读者。压缩包共5个文件均为cpp源码整体约9KB按功能覆盖CE318数据读取、气溶胶光学厚度与水汽含量的对比计算、仪器定标以及角度拟合等环节可作为课程实验或科研预处理的参考实现。已有1084人学习下载说明其在同类小众数据处理资料中具有一定参考价值。读者可从中了解观测数据清洗、辐射量比较与反演结果校核的组织方式并借助源码梳理AOD与wv计算的模块划分、关键参数传递和常见异常处理思路便于移植到Python或其他分析流程中。包体轻量、主题集中更适合作为算法原型和代码框架来查阅而非完整业务系统。1. 从一台 CE318 的观测文件说起AOD 和 wv 到底怎么落地凌晨三点自动观测站的 CE318 又传回一包 .dat 文件里面是八个波段、每几分钟一次的太阳直射电压。对做气溶胶和大气水汽的人来说这些数字真正的价值只有两个气溶胶光学厚度 AOD以及整层水汽含量 wv。CE318 的太阳光度计数据看起来门槛不高无非是电压除以定标常数、取对数、再除以大气质量数但真要把 AOD 做到能和公开基准产品对上、把 940 nm 通道的水汽反演稳定在可发布水平中间隔着定标、云筛选、瑞利订正、臭氧订正和查找表系数这几道坎。这篇笔记按我实际处理观测文件的顺序拆先讲原始数据和 triplet 质量控制再写 AOD 反演和 Angstrom 指数然后落到 940 nm 水汽通道最后把踩过的坑和交叉验证方法摊开。适合手头有 CE318 观测文件、准备自己写处理脚本的同行也适合刚接手站点数据、需要判断现有流程靠不靠谱的人。2. 数据预处理从原始电压到可用的观测量2.1 先认清 CE318 的测量序列和文件结构CE318 的一次观测不是单帧而是一组按时间排列的测量序列。典型文件里包含日期、时间、太阳天顶角、各波段直射电压以及仪器内部温度、气压等辅助量。不同固件导出的列名和分隔符不一样常见的是空格或逗号分隔的文本字段顺序也随版本变化。我一般先做一件事拿一天的文件把头十行和尾十行打印出来确认时间格式、波段顺序、有没有表头、有没有以#开头的注释行。这一步不做后面列名对不上AOD 会整列错位而且错得很隐蔽。处理之前要把三个辅助量准备好观测点的气压、臭氧柱总量、日地距离修正因子。气压影响瑞利光学厚度臭氧影响短波波段日地距离影响定标常数的换算。CE318 文件里通常有气压字段但如果没有就要从站点气象数据补。臭氧可以用当天的卫星或再分析资料精度不需要极高但 440 nm 附近不能忽略。日地距离因子按儒略日计算常见做法是用简化的椭圆轨道公式。这些量不需要手工算写成函数挂在读取流程后面即可。另一个容易被忽略的点是时间。CE318 记录的是 UTC 还是本地时间不同站点配置不同。AOD 反演本身对时间不敏感但一旦要和公开基准产品逐点对比时区错一小时就会让匹配失败。我的习惯是读入后先统一转成 UTC并在数据框里保留原始时间列方便回溯。时间戳解析用pandas.to_datetime时显式指定utcTrue不要依赖默认本地时区。2.2 云筛选和 triplet 质量控制CE318 的 AOD 反演里云是最大的污染源。薄云会让电压偏高AOD 偏低厚云直接让信号饱和或掉零。单帧判断云很难所以厂家和公开处理流程都用 triplet 方法连续三次测量作为一个 triplet三个值之间的变异系数超过阈值就整组丢弃。常见阈值是 1% 到 5%具体看波段和站点气候。我一般先按 2% 跑一遍再看剩余数据量如果剩余太少再放宽到 5%但要在日志里注明。triplet 的实现不复杂按时间排序每三条一组计算每个波段电压的标准差除以均值。注意不是对 AOD 做而是对原始电压做因为电压更接近观测值本身。还要处理一种情况三次测量中间缺了一次不能硬凑直接跳过这组。另一个细节是太阳天顶角太低的太阳高度角对应太大的人气质量数AOD 误差会被放大通常把天顶角大于 75° 到 80° 的数据剔掉。这个阈值不要一刀切低纬度站点可以放到 80°高纬度冬季可能只能用到 70°。除了 triplet还要做电压范围检查。每个波段的定标常数 V0 是已知的测量电压如果大于 V0 的 1.2 倍基本可以判定异常小于 0 或等于 0 的直接丢。这个检查能挡住大部分仪器饱和和遮光问题。做完这两步剩下的数据才进入 AOD 计算。2.3 读取与清洗的 Python 实现下面这段代码假设你已经把原始文件整理成 CSV列名统一为datetime, airmass, v440, v675, v870, v940, v1020。实际文件列名不同改rename映射即可。代码里保留了原始时间列并做了 triplet 标记。import numpy as np import pandas as pd # 读取已导出的 CSV原始文件如果是空格分隔把 sep 改成 r\s df pd.read_csv(ce318_day.csv, parse_dates[datetime]) df[datetime] pd.to_datetime(df[datetime], utcTrue) df df.sort_values(datetime).reset_index(dropTrue) # 基本电压检查负值、零值、超过定标常数 1.2 倍的标记为无效 V0 {v440: 12500.0, v675: 9800.0, v870: 7600.0, v940: 6200.0, v1020: 5400.0} for col, v0 in V0.items(): df[f{col}_valid] (df[col] 0) (df[col] 1.2 * v0) # triplet 质量控制每三条一组计算每个波段的变异系数 df[triplet_id] np.arange(len(df)) // 3 triplet_stats df.groupby(triplet_id).agg( n(datetime, count), cv440(v440, lambda x: x.std() / x.mean() if x.mean() 0 else np.nan), cv675(v675, lambda x: x.std() / x.mean() if x.mean() 0 else np.nan), cv870(v870, lambda x: x.std() / x.mean() if x.mean() 0 else np.nan), cv940(v940, lambda x: x.std() / x.mean() if x.mean() 0 else np.nan), ) # 只有三条齐全且四个波段变异系数都小于 0.02 的 triplet 保留 good_triplet triplet_stats[ (triplet_stats[n] 3) (triplet_stats[[cv440, cv675, cv870, cv940]] 0.02).all(axis1) ].index df[triplet_ok] df[triplet_id].isin(good_triplet) # 太阳天顶角筛选这里用 airmass 反推airmass 4 约对应天顶角 75° 以上 df[geometry_ok] df[airmass] 4.0 df_clean df[df[triplet_ok] df[geometry_ok]].copy() print(f原始 {len(df)} 条清洗后 {len(df_clean)} 条)这段代码的逻辑分三层。第一层是电压合法性V0是示意值实际必须换成你仪器定标文件里的数不能直接抄。第二层是 triplet按行号每三条一组要求三条齐全且各波段变异系数小于 0.02cv计算的是标准差除以均值均值接近零时返回NaN避免除零。第三层是几何筛选airmass 4.0是经验阈值对应太阳天顶角大约 75°低纬度站点可以放到 4.5高纬度冬季建议收到 3.0。清洗后不要急着算 AOD先把每个波段剩余数据量打印出来如果 440 nm 只剩个位数说明当天云太多或者阈值太严这份数据不适合做日平均。3. AOD 反演Beer-Lambert、定标和 Angstrom 指数3.1 从电压到光学厚度的公式与定标常数CE318 测的是太阳直射电压满足 Beer-Lambert 定律的离散形式V(λ) V0(λ) * (d0/d)^2 * exp(-m * τ(λ))。其中V0是定标常数d0/d是日地距离修正m是大气质量数τ是总光学厚度。反解出来就是τ ln(V0 * (d0/d)^2 / V) / m。总光学厚度里包含气溶胶、瑞利散射、臭氧吸收、二氧化氮吸收水汽吸收在 940 nm 单独处理不放进 AOD 波段。要求 AOD就把非气溶胶项减掉τ_aer τ_total - τ_Ray - τ_O3 - τ_NO2。定标常数V0是整个流程的命门。它通常来自厂家定标文件也可能来自站点自己做的 Langley 定标。不同仪器、不同波段、不同定标日期V0都不一样。我见过最常见的翻车是脚本里硬编码了一组V0换了一台仪器还在用AOD 整体偏移 0.05 到 0.1肉眼看曲线还挺正常和基准产品一比就露馅。正确做法是把V0放在单独的配置文件里文件名带仪器序列号和定标日期脚本启动时校验序列号是否匹配。Langley 定标不是本文重点但要知道如果V0是自己拟合的必须用晴空稳定日、大气质量数跨度足够大的数据且拟合残差要检查。日地距离修正因子按观测时刻的儒略日算常见近似公式是(d0/d)^2 1.00011 0.034221*cos(g) 0.00128*sin(g) 0.000719*cos(2g) 0.000077*sin(2g)其中g 2π*(day_of_year-1)/365.25。这个因子在 1 月附近约 1.0347 月附近约 0.967量级 3%对 AOD 的影响在小光学厚度时不能忽略。很多脚本图省事直接取 1短波波段会带来 0.01 左右的系统偏差。3.2 气溶胶光学厚度计算与波长插值AOD 计算按波段逐个做。先算总光学厚度再减瑞利、臭氧、二氧化氮。瑞利光学厚度用近似公式τ_R 0.008735 * λ^(-4.08) * P/1013.25λ单位是微米P是站点气压。这个公式在 440 到 1020 nm 范围内够用短波 340 nm 附近误差会大一些如果做 340 nm 要换更精确的 Bodhaine 参数化。臭氧光学厚度用τ_O3 k_O3(λ) * ozone_DU / 1000吸收系数k_O3随波长查表440 nm 附近约 0.0003 每 DU 量级675 nm 以上基本可以忽略。二氧化氮在污染站点才需要认真订正清洁站点可以设零但要在文档里写明。下面是一个按波段计算 AOD 的函数定标常数和系数都从字典传入不硬编码在公式里。import numpy as np # 臭氧吸收系数示意单位每 DU 的光学厚度实际用查表值替换 O3_COEF {440: 0.00032, 675: 0.00004, 870: 0.00001, 1020: 0.00000} def rayleigh_od(wavelength_nm, pressure_hpa): 瑞利光学厚度近似波长单位 nm气压单位 hPa lam_um wavelength_nm / 1000.0 return 0.008735 * lam_um ** (-4.08) * pressure_hpa / 1013.25 def no2_od(wavelength_nm, no2_du): 二氧化氮光学厚度清洁站点可传 0 # 简化系数实际按波段查表 coef {440: 0.00012, 675: 0.00002, 870: 0.00000, 1020: 0.00000} return coef.get(wavelength_nm, 0.0) * no2_du def compute_aod(row, v0, pressure_hpa, ozone_du, no2_du0.0): 输入一行观测返回各波段 AOD 字典 out {} for wl in [440, 675, 870, 1020]: v row[fv{wl}] if not np.isfinite(v) or v 0: out[wl] np.nan continue # 日地距离修正已在 row 中预先算好列名 earth_sun_factor tau_total np.log(v0[fv{wl}] * row[earth_sun_factor] / v) / row[airmass] tau_aer (tau_total - rayleigh_od(wl, pressure_hpa) - O3_COEF.get(wl, 0.0) * ozone_du - no2_od(wl, no2_du)) out[wl] tau_aer return out这段代码里row[earth_sun_factor]是提前算好的(d0/d)^2airmass是大气质量数。v0字典必须来自当前仪器的定标文件pressure_hpa用观测时刻的站点气压ozone_du用当天臭氧柱总量。O3_COEF和no2_od里的系数只是示意实际要换成你所用波段宽度和滤光片对应的查表值。函数返回NaN而不是抛异常是为了让批量处理能继续跑但事后要统计NaN比例。如果 440 nm 的 AOD 出现负值且绝对值大于 0.05先别急着当噪声丢检查定标常数和气压是不是用错了。3.3 Angstrom 指数与波长依赖的校验AOD 本身是波长相关的常用 Angstrom 指数描述这种依赖τ_aer(λ) β * λ^(-α)α就是 Angstrom 指数。用 440 nm 和 870 nm 两个波段算α -ln(τ440/τ870) / ln(440/870)。α大于 1 通常表示细模态气溶胶占优小于 0.5 表示粗模态占优沙尘或海盐场景常见。这个指数不参与水汽反演但它是检查 AOD 质量的好工具如果某天α突然从 1.2 掉到 0.1或者出现负值大概率是某个波段的定标或云筛选出了问题。我一般会把α和τ440画成散点图正常站点会形成一个有物理意义的带状分布异常点会飘在外面。还有一个检查是波长插值用 440 和 870 的 AOD 去预测 675 的 AOD看残差是否在 0.02 以内。如果 675 nm 的残差系统性偏大可能是 675 波段的定标常数有问题或者臭氧订正用错了。这些检查不需要额外数据用已有的三个波段就能做。Angstrom 指数对水汽反演也有间接影响940 nm 通道的 AOD 不能直接用来算水汽但如果你要检查 940 nm 的信号是否合理可以用α外推一个 940 nm 的 AOD 参考值和 940 nm 的总光学厚度比看看水汽吸收项是不是大得离谱。这个技巧在后面水汽章节还会用到。4. 水汽含量 wv 反演940 nm 通道的查找表与经验系数4.1 940 nm 水汽通道的物理基础940 nm 落在水汽吸收带内CE318 这个波段的电压同时受气溶胶散射、瑞利散射和水汽吸收影响。水汽反演的核心思路是先假设 940 nm 的气溶胶光学厚度可以用邻近波段外推再扣掉瑞利和臭氧剩下的吸收光学厚度就归给水汽。实际处理里更常用的是经验公式或查找表因为 940 nm 的滤光片带宽、仪器响应函数和水汽吸收线型都会影响结果纯解析计算很难做准。常见做法有两种。第一种是厂家配套软件里的经验公式形式大致是wv (1/m) * (ln(V0_940/V_940) / a)有的还带二次项或气压修正。第二种是查找表用辐射传输模式预先算好不同水汽含量、不同天顶角、不同气溶胶条件下的大气透过率反演时查表插值。两种方法没有绝对优劣经验公式快、适合批量查找表准、但依赖模式输入。我一般先用经验公式跑一版再用查找表或公开基准产品抽检偏差超过 10% 就回头查系数。不管用哪种方法V0_940都必须来自同一次定标且不能和 AOD 波段的V0混用。940 nm 的定标常数对水汽反演的影响是系统性的V0偏 2%wv 可能偏 5% 以上因为公式里是取对数后再除以吸收系数误差会被放大。另一个关键是大气质量数水汽反演用的大气质量数和水汽吸收的质量数不完全一样经验公式里通常仍用几何大气质量数但低太阳高度角时误差明显所以水汽反演也要做天顶角筛选。4.2 查找表法和经验公式的实现下面给一个经验公式的批量实现包含二次项修正和异常值过滤。公式结构是wv (1/m) * (ln(V0/V) / a) b * (ln(V0/V))^2其中a和b是仪器和波段相关的系数必须从定标文件或厂家文档里取代码里的数值只是占位不能直接用于生产。import numpy as np def retrieve_wv(row, v0_940, a_coef, b_coef, max_wv8.0): 经验公式反演水汽含量。 row 需要包含v940, airmass, earth_sun_factor v0_940: 940 nm 定标常数 a_coef, b_coef: 经验系数来自定标文件 max_wv: 物理上限超过视为异常 v row[v940] m row[airmass] if not np.isfinite(v) or v 0 or m 0: return np.nan # 日地距离修正后的等效电压比 ratio v0_940 * row[earth_sun_factor] / v if ratio 1.0: # 电压高于定标常数可能是云或仪器异常 return np.nan ln_ratio np.log(ratio) # 经验公式一次项加二次项修正 wv (ln_ratio / a_coef b_coef * ln_ratio ** 2) / m if not np.isfinite(wv) or wv 0 or wv max_wv: return np.nan return wv代码里ratio 1.0的判断很关键正常情况下V0 * earth_sun_factor应该大于测量电压因为大气有消光如果反过来说明测量电压异常偏高常见原因是薄云增强了前向散射或者仪器刚做过定标但常数没更新。max_wv8.0是经验上限实际站点水汽很少超过 6 cm超过 8 基本可以判为异常。a_coef和b_coef的取值因仪器而异不要抄网上的通用值去找你仪器最近一次定标报告。如果找不到二次项系数可以令b_coef0先用线性项跑但要在结果里标注方法版本。批量跑的时候不要只输出 wv 一列。把ln_ratio、airmass、earth_sun_factor和原始电压一起保留后面排查异常时能直接看出是哪一项出的问题。我习惯再加一列wv_flag标记ratio1、wvmax_wv、airmass4等不同原因统计时按 flag 分组比只看NaN有用。4.3 结果对比与异常值处理水汽反演结果要和至少一种独立来源对比。常见的是公开观测网的 Level 2.0 产品或者同站点的微波辐射计、探空资料。对比时不能只比日均值要按时间窗口匹配比如 ±30 分钟内的均值否则日变化会被抹平。匹配后算偏差和均方根误差偏差在 10% 以内、均方根误差在 0.3 cm 以内一般可以接受如果偏差是系统性的优先查V0_940和吸收系数如果误差是随机的查云筛选和天顶角阈值。异常值处理要克制。水汽在时间序列上连续变化突然出现一个 0.2 cm 或 7 cm 的点大概率是云或仪器问题可以直接丢。但如果连续几个小时都偏低不要逐点丢先看是不是气压或臭氧订正用了常数、实际当天变化较大。还有一种情况是 940 nm 电压接近饱和ln_ratio很小反演出的 wv 接近零这种点会被wv0过滤掉但过滤比例高的时候要检查是不是仪器增益设置不对。5. 避坑与排查AOD 和 wv 反演里最容易翻车的五件事5.1 定标常数用错整条曲线偏移现象AOD 时间序列形状正常但和基准产品比整体偏高或偏低 0.05 以上Angstrom 指数也偏得不多。原因通常是V0用了别的仪器、别的定标日期或者日地距离因子没乘。解决把V0从脚本里挪到配置文件文件名带序列号和定标日期启动时校验日地距离因子单独写函数单元测试用一个已知日期对一遍。我现在的习惯是每次换仪器先跑一天数据和基准产品比日均值偏差超过 0.02 就不往下做。5.2 triplet 阈值一刀切晴空数据被误杀现象清洗后数据量骤减尤其是冬季或高纬度站点剩下不到原始量的 10%。原因是变异系数阈值 0.02 对低信号波段太严或者三条测量里有一条缺测被整组丢掉。解决按波段设不同阈值440 nm 可以严一点1020 nm 放宽到 0.05缺测的 triplet 不要硬凑但可以标记后单独评估。另一个办法是先按天统计剩余率连续多天低于 20% 就回头检查阈值不要直接认为天气差。5.3 水汽公式漏了二次项或气压修正现象wv 在低水汽时和基准产品对得上高水汽时系统性偏低偏差随 wv 增大而增大。原因是经验公式只用了线性项高水汽时吸收饱和线性近似不够。解决找定标报告里的二次项系数没有就改用查找表同时检查气压修正有些公式里的m需要乘P/1013.25。如果二次项系数来源不明宁可不加也不要用网上抄来的通用值否则低水汽段会被带偏。5.4 时间戳时区错位匹配全失败现象AOD 和 wv 自身看着正常但和基准产品按时间匹配后一条都对不上或者匹配上的都是错误时段。原因是原始文件用本地时间脚本按 UTC 解析差了几个小时。解决读入时保留原始时间字符串显式转换时区输出文件里同时写 UTC 和本地时间两列匹配前先打印前五条时间戳和基准产品对一眼。这个坑不复杂但一旦发生排查成本很高因为数据本身没错错的是索引。5.5 把云污染当晴空AOD 偏低 wv 偏高现象AOD 出现不合理的低值同时 wv 偏高Angstrom 指数跳变。原因是薄云增强了 940 nm 的前向散射电压偏高ln_ratio偏小wv 反而被算大AOD 波段电压偏高则让 AOD 偏低。解决除了 triplet再加一道“相邻波段一致性”检查用 440、675、870 的 AOD 拟合 Angstrom 指数残差超过 0.05 的整组标记水汽那边检查ratio是否接近 1接近 1 的直接丢。云筛选没有万能阈值只能多层过滤加人工抽检。6. 从单站到批量用公开基准产品做交叉验证和自动化单站跑通只是第一步真正要拿这份数据处理长期观测得把批量脚本和交叉验证固定下来。我的做法是每天定时跑前一天的文件输出三张表清洗后的观测表、AOD 表、wv 表每张表都带质量标记和版本号。版本号不是形式主义定标常数、系数、阈值任何一项改了都要升版本否则半年后回头看根本分不清哪批结果用的是什么参数。批量脚本里加一个config.yaml把仪器序列号、定标日期、V0字典、triplet 阈值、天顶角上限、水汽系数全写进去脚本启动时打印配置摘要日志里留一份。交叉验证我一般按周做不按天。按天做容易被个别云污染点带偏按周看偏差和趋势更稳。匹配窗口用 ±30 分钟先把基准产品的时间转成 UTC再按最近邻匹配匹配不上就跳过不要插值。对比指标至少看三个平均偏差、均方根误差、相关系数。AOD 的相关系数低于 0.95 就要查水汽的均方根误差大于 0.4 cm 就要查。查的顺序是先看时间匹配数量再看异常点分布最后才动公式。很多时候问题不在公式在匹配和筛选。下面是一个批量交叉验证的骨架重点在匹配和统计不在绘图。import pandas as pd import numpy as np def match_nearest(df_a, df_b, time_coldatetime, tolerance30min): 按最近时间匹配两个数据框返回匹配后的合并表 df_a df_a.sort_values(time_col).copy() df_b df_b.sort_values(time_col).copy() df_a[time_col] pd.to_datetime(df_a[time_col], utcTrue) df_b[time_col] pd.to_datetime(df_b[time_col], utcTrue) matched pd.merge_asof( df_a, df_b, ontime_col, directionnearest, tolerancepd.Timedelta(tolerance), suffixes(_mine, _ref) ) return matched.dropna(subset[AOD440_ref]) def summary_stats(matched, col_mine, col_ref): 计算偏差、均方根误差、相关系数 x matched[col_mine].astype(float) y matched[col_ref].astype(float) mask np.isfinite(x) np.isfinite(y) x, y x[mask], y[mask] bias (x - y).mean() rmse np.sqrt(((x - y) ** 2).mean()) corr np.corrcoef(x, y)[0, 1] if len(x) 2 else np.nan return {n: len(x), bias: bias, rmse: rmse, corr: corr}pd.merge_asof要求两个表都按时间排序directionnearest是最近邻tolerance30min控制最大允许间隔。匹配后一定要看n如果一周只剩几十个点先别分析偏差回去查时间格式和筛选阈值。summary_stats里先做有限值过滤避免NaN污染统计。这个骨架不绘图但输出可以直接写进日志或周报。我的习惯是每周一早上跑一遍上周的匹配把n、bias、rmse、corr四列贴在一张表里连续几周bias漂移就说明定标或系数该更新了。还有一个小技巧把 AOD 和 wv 的交叉验证结果按天气类型分组晴天、薄云、沙尘各看一遍。整体指标好看不代表每种场景都对沙尘天的 Angstrom 指数低水汽反演的云筛选阈值可能偏松分组之后问题更容易暴露。从那以后我每次换定标常数都强制走一遍“单日对比 → 周统计 → 分组检查”这三步不再只看一条时间序列顺眼就发布。希望帮到你。本文还有配套的精品资源点击获取