简介利用CE318太阳光度计观测数据反演气溶胶光学厚度AOD与水汽含量WV的C源码包面向大气科学、遥感反演方向的科研人员与学习者解决从原始辐射观测到大气参数反演的完整流程。压缩包共5个C源文件仅9KB代码简洁功能模块清晰涵盖原始数据读取、仪器标定、AOD拟合计算可参考Klett法、Fernald法等以及基于940nm吸收带的水汽反演并兼顾太阳位置修正等细节。目前已有1082人学习下载该资源适合作为算法参考帮助读者快速理解CE318数据处理的各环节或直接移植用于站点观测数据处理与教学实验。通过阅读代码可掌握从数据清洗、辐射传输计算到质量控制的关键实现是一份实用的大气光学编程范例。1. 为什么外场观测都在等这份AOD和wvCE318的定位与数据链路做气溶胶外场观测的人大概率都碰过CE318太阳光度计。它测的是太阳直射辐照度但真正想要的是AOD气溶胶光学厚度和大气水汽含量wv这两个反演量——一个决定消光与辐射强迫估算一个影响大气校正和天气系统分析。问题是CE318原始记录里只有各波段的DN值和太阳跟踪器的角度AOD和wv不会自动出现在屏幕上。这篇笔记要把从原始直射数据到AOD和wv的完整链路讲透定标怎么处理、瑞利散射和气体吸收怎么扣、940nm水汽通道为什么不能照搬AOD公式、以及批量处理时最容易翻车的几个细节。适合课题组刚接手光度计的学生也给环境监测站和遥感团队一条可复现的处理路径。2. 先把原始信号洗干净定标、暗电流与Langley截距2.1 CE318的直射观测里有什么字段、通道和可用数据类型CE318的观测模式分直射太阳Direct Sun、平纬环扫描Almucantar和主平面扫描Principal Plane。反演AOD和wv只用到直射数据其余两种是反演气溶胶粒径分布和相函数用的不要在AOD这段混进来处理。直射观测文件里常见的字段是时间戳、太阳天顶角SZA、太阳方位角AZM以及每个通道的DN值或电压值。老一代CE318的气溶胶通道集中在340、380、440、500、675、870、1020nm水汽通道在936nm附近新一代CE318-T把水汽通道放在940nm。拿到数据先做一件事把通道波长落到你手里的定标文件上不要想当然用统一波长型号不同、滤光片批次不同都会有零点几个纳米的差异。数据导出的常见做法是用光度计自带的处理软件把原始文件批量导出成CSV或者写个小脚本直接解析二进制记录。前者省事后者灵活。无论哪条路先确认两件事——DN值有没有扣暗电流、通道带没有带增益档标记。这两项没确认就往下算后面查负AOD会查到怀疑人生。CE318通道与用途大体如下表处理前先按这个把数据列对齐。通道中心波长主要用途反演中的角色340/380nm紫外-紫光细粒子敏感参与Angstrom指数计算440/500nm蓝-绿光标准AOD反演臭氧吸收修正重点675/870nm红-近红外粗粒子与插值基准瑞利扣除小信噪比好1020nm近红外粗粒子通道常与870做双波长对比936/940nm水汽吸收带反演wv不参与AOD需单独定标2.2 Langley定标从DN值到V0为什么晴稳天高山站最可靠AOD反演的理论根基是Beer-Lambert定律。仪器测到的直射信号V(λ)可以写成定标系数V0(λ)、日地距离项和大气透过率的乘积。V0是这台仪器在大气层顶、日地平均距离处对太阳辐射的响应值它不能从说明书里抄必须靠定标得到。常见的定标方式有两种。第一种是Langley定标选一个大气清洁、气溶胶光学厚度稳定且偏小的晴天在整个上午或下午连续对太阳直射观测此时气溶胶光学厚度随时间近似不变把ln(V)对大气质量m做线性回归截距就是ln(V0)。这条线拟合得好不好直接决定你后面所有AOD的准确性。第二种是交叉定标把待定标仪器和一台定标过的参考仪器并排架在同一个站点用同一时段直射数据反推V0。我一般优先用Langley除非场地条件实在达不到。Langley定标最怕两件事气溶胶光学厚度在观测时段内漂移以及天空中出现肉眼几乎看不见的薄云。所以定标通常选高山清洁站或者连续晴稳天。定标完的V0要注意时间有效期滤光片会老化、窗口会脏放半年的V0直接套用AOD负值的概率非常大。血泪经验是每次外场实验前把仪器发回定标或者并排观测交换定标系数。2.3 大气质量与日地距离修正两个逃不掉的几何参数从DN值到AOD中间要经过两次几何修正。第一是大气质量m也就是太阳直射光穿过大气的相对路径长度。常看到有人用最简单的1/cos(SZA)但天顶角超过60度时这个近似会带来明显偏差。我一般用Kasten-Young公式m 1 / (cos(SZA) 0.50572 * (96.07995 - SZA)^(-1.6364))注意这里的SZA是度数不是弧度。这个公式在SZA小于80度时都能给出可靠结果算完AOD再回头看晴天序列会平滑很多。第二是日地距离修正。地球绕太阳公转导致日地距离在一年里有大约3.5%的变化这部分如果不修正会直接叠加进AOD的系统偏差里。定义F (d_mean / d)^2其中d_mean是平均日地距离d是观测当天的实际距离。实际计算时d可以通过年积日N用近似公式估算d d_mean * (1 - 0.0167 * cos(2π * (N - 3) / 365))然后令D_factor (d_mean / d)^2。这个数值在1月初约1.0347月初约0.967处理全年数据时不是小量。2.4 用Python读取CE318直射数据并做基础清洗结合前几节的参数设置读数据进行基础清洗的代码我写成这样import pandas as pd import numpy as np def load_ce318_direct(csv_path): df pd.read_csv(csv_path, parse_dates[datetime]) # 只保留直射观测行tracking_status0 表示正常跟踪 df df[df[tracking_status] 0].copy() # 扣除暗电流CE318 输出的是带暗底的计数时这一步不能省 dark_cols [c for c in df.columns if dark in c.lower()] for col in dark_cols: df[col.replace(_dark, )] df[col.replace(_dark, )] - df[col] sza df[sza].values sza_rad np.deg2rad(sza) # Kasten-Young 大气质量SZA 超过 80 度直接丢弃 m_val 1.0 / (np.cos(sza_rad) 0.50572 * (96.07995 - sza) ** (-1.6364)) df[airmass] m_val df df[(df[airmass] 15) (df[airmass] 1)] day_of_year df[datetime].dt.dayofyear.values r 1.0 - 0.0167 * np.cos(2 * np.pi * (day_of_year - 3) / 365.0) df[distance_factor] (1.0 / r) ** 2 return df这段代码做的事情有三个。第一步按跟踪状态去掉跟踪异常的记录太阳光度计在云缝里来回找目标时DN值是不可信的这一步比任何滤波都重要。第二步做暗电流扣除CE318的探测器在低信号时暗底不可忽略特别是340和380nm两个短波通道。第三步算大气质量和日地距离因子并把SZA过大、路径太斜的数据挡在门外。参数上要注意airmass上限15对应的SZA大约在86度实际处理我常用8-10更保守distance_factor直接乘在后面的辐射量里即V_corrected V_raw * distance_factor这里要写清楚再往下算。3. 反演AOD逐通道扣除瑞利散射与气体吸收的完整计算3.1 总光学厚度的Beer-Lambert关系看懂反演公式再动手AOD反演的起点是把上一节清洗后的V(λ)和定标系数V0(λ)代入Beer-Lambert关系。总光学厚度τ_total定义为τ_total(λ) (1 / m) * [ ln(V0(λ)) - ln(V(λ)) ln(D_factor) ]这里的D_factor就是2.3节的日地距离修正项对应(d_mean/d)^2。注意这个公式里V0包含了大气层顶、日地平均距离下的定标响应所以ln(D_factor)要加回来符号反了会整体偏一个固定量。总光学厚度是瑞利散射、气溶胶消光、臭氧吸收、NO2吸收的叠加。我们要的是气溶胶那一项剩下三项必须从总量里扣出去。很多入门教程只提瑞利散射扣除把臭氧和NO2略过。短波段如果不管臭氧440nm处的AOD可能被高估0.01-0.03对干净大气来说这个偏差太明显了。所以完整流程是先算总光学厚度再逐项扣除分子散射和气体吸收剩余才是AOD。3.2 Rayleigh与臭氧、NO2扣除哪些通道必须做哪些可以忽略瑞利散射扣除的标准做法是用Bodhaine等人在1999年给出的经验公式压强归一化到观测站气压τ_R(λ) (P / P0) * 0.008735 * λ^(-4.08)其中λ以微米为单位P是观测站气压P0是标准大气压1013.25hPa。这个公式在340-1020nm范围内表现稳定比用单次散射近似更贴近真实大气。常用的几个波长下海平面高度的瑞利光学厚度参考值如下波长(nm)瑞利光学厚度(P01013.25hPa)3400.6603800.4234400.2495000.1486750.0438700.015410200.0068上表数值是直接用公式算的参考值适合做反演结果的合理性检查。如果某天440nm的AOD算出来是-0.05先把瑞利项复查一遍。臭氧修正主要在440nm附近。臭氧在Chappuis带的吸收会叠加在440-675nm区间其中500nm附近最明显。简化处理时用臭氧柱总量O3(DU)乘以吸收系数k(λ)臭氧吸收光学厚度τ_O3 k(λ) * O3。k的参考量级是440nm约0.003/DU500nm约0.03/DU675nm约0.001/DU。如果手里没有臭氧柱总量可以用全球气候态的月平均值或者从再分析资料里读取不要直接设零。NO2在城市和工业区影响明显但在清洁站通常可以忽略。反演时用440nm附近的非吸收通道做参照或者直接把NO2吸收项设为0。做城市观测时才有必要引入NO2柱总量修正这点在方案设计阶段就要想清楚否则后期追偏差异常困难。3.3 AOD与Angstrom指数计算代码输入输出和边界条件下面这段代码把前面所有公式串起来输入是清洗后的DataFrame和定标系数输出是逐通道AOD和Angstrom指数def calc_aod(df, v0_dict, wavelength_dict, station_pressure, o3_du): from scipy.interpolate import interp1d aod_out pd.DataFrame(indexdf.index) p_ratio station_pressure / 1013.25 for ch in v0_dict.keys(): lam wavelength_dict[ch] # 单位微米 v_raw df[fdn_{ch}].values d_factor df[distance_factor].values airmass df[airmass].values # 总光学厚度 tau_total (np.log(v0_dict[ch]) - np.log(v_raw) np.log(d_factor)) / airmass # 瑞利散射扣除 tau_rayleigh p_ratio * 0.008735 * (lam ** -4.08) # 臭氧修正中心波长越靠近 500nm 系数越大 k_o3 np.interp(lam * 1000, [340, 440, 500, 675], [0, 0.003, 0.03, 0.001]) tau_o3 k_o3 * o3_du # NO2 在城市观测时可加上这里先置零 aod_out[faod_{ch}] tau_total - tau_rayleigh - tau_o3 # Angstrom 指数440-870nm 双波长 aod_440 aod_out[aod_440].values aod_870 aod_out[aod_870].values valid (aod_440 0) (aod_870 0) aod_out[angstrom_440_870] np.nan aod_out.loc[valid, angstrom_440_870] ( -np.log(aod_440[valid] / aod_870[valid]) / np.log(440.0 / 870.0) ) return aod_out这段代码里v0_dict的每个通道必须和中心波长一一对应不要拿936nm的定标系数去算940nm的通道波长差几个纳米瑞利项和臭氧项的差别不大但后续水汽反演的误差会被放大。臭氧系数用np.interp线性插值如果只做科研级粗算也可以固定给0.003/DU。边界条件上AOD算出来小于0要标记而不是删除因为负值往往说明定标或者气体修正在某个通道有问题。Angstrom指数只在一对波长都有效时才计算负的或接近零的AOD会让指数变成无穷大或NaN。另外440和870的AOD对Angstrom指数的敏感性很高两个通道定标误差只要各有0.01指数就会偏0.05左右所以不要只看指数就说气溶胶类型。4. 反演大气水汽含量wv940nm通道的修正Beer-Lambert法4.1 为什么水汽通道不能照搬AOD公式饱和吸收与通道响应到了水汽反演很多照搬AOD公式的人会翻车。940nm通道测的不是简单消光水汽分子在这个波段有大量振动-转动吸收线而且CE318的940nm通道带宽约10nm通道内同时包含强线、弱线和连续吸收。Beer-Lambert定律在单色光条件下成立但宽通道里透过率与水汽总量的关系变成非线性直接用线性光学厚度公式会严重低估水汽。这一节是整个wv反演的关键。所以工程上普遍用修正的Beer-Lambert关系。940nm通道的信号可以写成V_w V0_w * D_factor * exp(-m * τ_a(940)) * exp(-m * a * wv^b)这里的a和b是经验参数描述通道有效透过率与水汽总量的幂律关系。wv的单位是cm表示整层大气可降水量。τ_a(940)是940nm处的气溶胶光学厚度但CE318在940nm没有单独的气溶胶通道必须从其他通道插值外推这是另一个容易踩坑的点。把上式取对数整理定义中间量LL (1 / m) * [ ln(V0_w) - ln(V_w) ln(D_factor) ] - τ_a(940) a * wv^b于是反演变成两步先算出L再解出wv (L / a)^(1/b)。整套流程核心不在算法多深而在参数匹配和中间量计算是否严格。4.2 从940nm透过率到可降水量wv指数参数a、b的标定参数a和b不是随便抄的实验室常数。它们取决于940nm通道的滤光片透过率函数、仪器响应和水汽吸收线强度分布。AERONET处理中有自己的一套标定参数一般a在0.65-0.70附近b在0.5-0.6附近但换一台仪器、换一个滤光片这两个值就要重新审视。我一般会做本地标定选若干个晴稳天把CE318反演结果与同站探空得到的水汽柱总量做回归。具体做法是固定b的值对(L, wv)数据做对数线性回归求a扫b在0.4-0.7范围找回归残差最小的组合。这样定出来的a、b对本地气候才靠谱。没有探空资料时用微波辐射计或GNSS水汽产品替代也可以但时间匹配要做准。参数敏感性上经验是这样的情况a偏大b偏大不经标定直接套用wv较小(1cm)结果偏低结果偏高偏差可能超过20%wv适中(1-3cm)结果偏低结果偏低偏差约10%-15%wv较大(3cm)结果偏低结果偏低偏差更大且不线性所以标定比反演公式本身重要得多。使用别人论文里的a、b值之前先看仪器型号和定标时间是否对得上。4.3 单次扫描的水汽反演代码气溶胶扣除与插值策略代码里最关键的是气溶胶光学厚度插值。我用440和870两个通道的AOD做Angstrom外推到940nm这种做法对粗模态气溶胶有误差但没有更直接的测量手段是通行做法def calc_wv(df, v0_wv, a_param, b_param): from scipy.interpolate import interp1d out pd.DataFrame(indexdf.index) # 先保证上一节算好的AOD在df里 aod_440 df[aod_440].values aod_870 df[aod_870].values lam_440, lam_870 440.0, 870.0 # Angstrom 外推 940nm 气溶胶光学厚度 valid (aod_440 0) (aod_870 0) tau_a_940 np.full_like(aod_440, np.nan) alpha np.full_like(aod_440, np.nan) alpha[valid] -np.log(aod_440[valid] / aod_870[valid]) / np.log(lam_440 / lam_870) tau_a_940[valid] aod_870[valid] * (940.0 / lam_870) ** (-alpha[valid]) v_w df[dn_940].values d_factor df[distance_factor].values airmass df[airmass].values # 计算 L 值 l_val (np.log(v0_wv) - np.log(v_w) np.log(d_factor)) / airmass - tau_a_940 # 反演可降水量单位 cm wv (l_val / a_param) ** (1.0 / b_param) out[wv_cm] wv out[tau_a_940_interp] tau_a_940 out[l_value] l_val return out这里有个前提前面的AOD计算函数要先运行并把aod_440和aod_870写回df。如果440nm或870nm的AOD出现负值插值结果就失真wv也会一起错。所以流程上必须先做AOD质量控制再算水汽顺序不能倒。水汽通道的V0_w也要单独定标。和普通气溶胶通道不同940nm通道Langley拟合的稳定性差因为大气水汽总量本身一直在变。一般做法是选多个晴天的数据段分段做Langley取截距的中位数或者直接依赖交叉定标传递。注意940nm通道的暗电流扣除比可见光通道更敏感因为水汽吸收强信号本身就低暗电流没扣干净会让wv整体系统性偏高。检查方式是用夜间数据看DN值是否归零。5. 数据处理避坑指南定标漂移、云判识与单位混乱的5个典型问题5.1 现象一AOD成片出现负值先查定标系数而不是查天气现象某个通道全天AOD都是负的早晨轻、中午重。原因几乎都是定标系数V0偏小或者滤光片窗口污染导致透过率下降。解决步骤是先看该通道的Langley定标时间超过半年就该怀疑再检查窗口镜片有没有结露、灰尘最后用晴稳天的数据重做Langley或与参考仪器交叉定标。负AOD不是一个数学问题是仪器状态问题。5.2 现象二wv整体偏低且湿度越大偏得越多现象wv反演结果比探空偏低夏天偏差能到20%。原因通常是a、b参数用的他人结果通道透过率函数不匹配。解决做法是把本地探空和光度计数据按时间匹配扫b参数做回归重新标定。另外检查τ_a(940)的插值是否用了被污染的气溶胶通道插值源本身有问题修正项就会引入系统偏差。5.3 现象三天顶角一大AOD和wv都开始跳变现象SZA超过60度后反演序列出现锯齿状波动。原因是大空气质量下平板大气近似失效且大气质量m的计算误差被放大。解决用Kasten-Young公式计算m并做SZA截断一般处理时不保留SZA大于70度的数据。如果必须保留要对m做球面大气修正否则宁可直接丢。5.4 现象四肉眼晴空但AOD偏高Angstrom指数异常偏小现象天空看起来很蓝但440nm AOD稳定在0.2以上且Angstrom指数小于0.8。原因往往是高空薄云或卷云进入了视场直射信号被削弱却没有触发跟踪异常。解决查看DN值时间序列若相邻两条记录信号波动超过3%-5%判定为可疑数据剔除再结合440nm和1020nm的AOD差值做一致性判断差值异常增大时标记为云影响。5.5 现象五两台CE318同站观测结果系统性偏差现象两台仪器并排架在同一个站房AOD差0.02-0.04wv差10%左右。原因多数是定标基准不同或者其中一台的滤光片透过率已经变化。解决做并排交叉定标用48小时连续直射数据传递V0确定偏差来源。处理时把两台仪器的V0统一到同一基准再重新反演。不要单独追算法差异算法一致时这种系统偏差九成在定标。5.6 数据处理技巧批量复算与验证指标反演过程越复杂越要保留中间变量。建议每条反演记录都输出airmass、distance_factor、tau_a_940、l_value这几个中间量后续排查会省很多时间。验证指标上用三个量平均偏差Bias、均方根误差RMSE和相关系数R。AOD与标准产品比对时440nm处偏差在±0.01以内是理想状态±0.02以内可接受wv与探空比对时相对偏差15%以内算合理。最后说句教训整套流程里最容易出错的不是公式而是定标系数和单位换算我把定标文件单独存档、每次反演前核对时间戳这习惯救过我好几次。希望帮到你。本文还有配套的精品资源点击获取