简介面向大气科学、环境遥感领域的研究者与学习者这份压缩包以太阳光度计原始观测数据为输入提供了反演气溶胶光学厚度与水汽含量的完整编程实现。资源共5个文件全部为C源码整体大小仅9KB包含观测数据读取、太阳位置角度拟合、气溶胶反演、水汽计算以及仪器定标等算法模块代码短小精悍便于逐段研读、修改与二次开发。已有1082人学习下载适合具备一定大气光学和编程基础、希望从实测辐射数据独立计算出大气参数的科研人员与学生使用。阅读这些源码可深入理解利用直接辐射与散射辐射推算气溶胶光学厚度的辐射传输思路掌握水汽吸收波段的反演处理方式同时体会太阳位置计算、大气折射修正等定标步骤对结果精度的重要影响。整套源码覆盖了从原始数据到最终参数的完整处理链路可作为相关课题的算法验证参考、教学演示样例或进一步开发的基础框架对开展气溶胶与水汽长期观测分析具有实用价值。1. 用CE318太阳光度计反演AOD和水汽wv数据链路比算法更值得花时间用CE318太阳光度计做AOD和水汽wv反演多数人以为核心是套公式实际动手后发现最大的成本在数据治理先要做Langley定标拿到仪器响应常数V0再做瑞利散射、臭氧吸收、质量数修正最后才轮到逐通道回归水汽通道还要单独处理吸收带。整个过程里90%的返工不是公式错而是定标漂移、云干扰、时间基准混乱这些不起眼的问题。这篇就按我从原始导出文件到最终AOD、wv产品的完整链路来写覆盖原理、可复现脚本、参数边界和常见坑适合正在维护地面太阳光度计、做气溶胶遥感校验或把地基数据拿去印证卫星产品的工程师。2. CE318原始数据预处理先把DN值、UTC和太阳天顶角整理干净2.1 CE318观测链路与通道用途CE318本质上是一台自动太阳跟踪光度计太阳辐射进入准直镜筒经滤光片轮分光后由探测器转成电信号最终输出为各通道的DN值。每完成一轮滤光片扫描就得到一条包含全部通道的观测记录。滤光片轮按固定波长序列工作常见六通道配置是440、500、675、870、936、1020 nm部分增强型号还会扩展340、380或1640 nm通道。六个通道在反演里的角色并不相同。440、500、675、870、1020以气溶胶散射贡献为主用于AOD反演936 nm位于水汽强吸收带用于反演整层水汽含量wv。这意味着一次完整的数据处理要跑两套链路非吸收通道走Langley定标加逐通道AOD回归936通道单独做水汽透过率换算。建议先把这两条链路的边界划清楚后面才不会被数据搞乱。下面这个表是我处理CE318数据时默认的通道分工新手可以直接照着建字段。通道(nm)主要大气过程反演产物440气溶胶散射AOD核心通道500气溶胶散射加轻微臭氧吸收AOD675气溶胶散射加微弱臭氧吸收AOD870气溶胶散射AOD、Angstrom指数锚点936水汽强吸收带水汽柱含量wv1020气溶胶散射AOD、水汽通道插值锚点2.2 读取导出文件并重建时间轴CE318配套上位机一般会把数据导出为带表头的文本文件字段通常包含UTC时间、各通道DN值、内部温度、跟踪偏差等。拿到文件后第一件事不是挑通道算AOD而是把时间轴和太阳几何参数重建好。很多用户直接拿导出的DN值开算却不检查时间是不是UTC漏掉时差方程和日地距离修正最后Langley回归斜率都是乱的。我一般用Python读文件后统一重建三个字段年积日、太阳天顶角、日地距离修正因子。太阳天顶角用经纬度和UTC时间算不直接信任仪器输出因为部分固件对时间的内部换算逻辑并不透明。下面这段代码可以直接复用。import numpy as np import pandas as pd raw pd.read_csv(ce318_export_001.csv, parse_dates[utc_time]) def sun_zenith_deg(lat_deg, lon_deg, dt): 由站点经纬度和UTC时间计算太阳天顶角度。 doy dt.dayofyear frac_hour dt.hour dt.minute / 60.0 dt.second / 3600.0 # 时差方程近似et单位是分钟 b 2.0 * np.pi * (doy - 81.0) / 364.0 et 9.87 * np.sin(2.0 * b) - 7.53 * np.cos(b) - 1.5 * np.sin(b) solar_time frac_hour lon_deg / 15.0 et / 60.0 hour_angle 15.0 * (solar_time - 12.0) decl np.radians(23.45 * np.sin(np.radians(360.0 * (284.0 doy) / 365.0))) lat_r, ha_r np.radians(lat_deg), np.radians(hour_angle) cos_zen (np.sin(lat_r) * np.sin(decl) np.cos(lat_r) * np.cos(decl) * np.cos(ha_r)) return np.degrees(np.arccos(np.clip(cos_zen, -1.0, 1.0))) def earth_distance_factor(doy): 日地距离修正因子d在平均距离附近小幅摆动。 return 1.0 - 0.01673 * np.cos(2.0 * np.pi * (doy - 3.0) / 365.25) raw[doy] raw[utc_time].dt.dayofyear raw[sza_deg] [sun_zenith_deg(LAT, LON, x) for x in raw[utc_time]] raw[d_fac] raw[doy].apply(earth_distance_factor)这段代码里LAT和LON要替换成站点真实经纬度东经为正、西经为负。solar_time是地方太阳时hour_angle以正午12点为0上午为负下午为正。时差方程的精度对日常反演足够太阳天顶角误差一般在0.1度以内换算到AOD影响小于0.001对水汽反演的影响会更敏感一些所以时间戳必须严格用UTC。2.3 反演前的三层质控原始记录不能直接全部进入AOD计算。我会先做三层预筛把明显无效的数据剔除云过滤放到后面算出AOD之后再做。原因很简单云判据依赖AOD时间序列的方差和相邻点跳变单靠DN值本身无法可靠识别薄云。第一层是太阳几何。太阳天顶角大于85度的数据直接标记无效大于80度的数据可以保留但不参与Langley定标。原因是气团数公式在80度以后开始变得敏感角度误差会被放大。第二层是跟踪偏差。CE318虽然自动跟踪太阳但薄云或阵风会让跟踪器短暂失锁导出的跟踪误差字段绝对值偏大时对应记录直接剔除。第三层是异常DN值。DN恒为常量、出现重复值或突然跳变到量级以外说明滤光片轮或电子学有异常也标记掉。qc raw.copy() qc[flag_geom] 0 qc.loc[qc.sza_deg 80.0, flag_geom] 1 qc.loc[qc.sza_deg 85.0, flag_geom] 2 qc[flag_track] (qc.track_err.abs() 0.05).astype(int) qc[flag_dup] qc.duplicated(subset[dn_440, dn_870], keepFalse).astype(int) qc[flag_invalid] qc[flag_geom].ge(1) | qc[flag_track].eq(1) | qc[flag_dup].eq(1)这里的track_err是导出文件里的跟踪偏差字段单位不同型号可能不同使用前先确认量纲。flag_geom等于2的数据我通常直接删flag_geom等于1的数据保留但在后续加权时不作为主样本。注意不要在这一步把疑似云的数据删光否则后续云过滤和Langley定标都会缺样本这是我反复体会到的一个教训。3. AOD反演链路Langley定标、质量数与瑞利扣除的最小实现3.1 反演公式怎么来的CE318直接测量的是到达仪器的太阳辐射通量在晴空条件下可以用比尔-朗伯定律的简化形式表达V V0 × d² × exp(-τ × m)其中V是实测DN值V0是大气质量数为零时仪器应有的响应常量d是日地距离修正因子τ是整层大气总光学厚度m是大气质量数。对非吸收通道总光学厚度拆成瑞利散射τ_R、气溶胶τ_a和臭氧τ_O3三部分于是AOD反演公式写为τ_a [ln(V0) 2×ln(d) - ln(V) - τ_R×m - τ_O3×m] / m这里面ln(V0)2ln(d)就是当前日地距离下的大气外响应ln(V)是实测信号减掉瑞利和臭氧后剩下的光学厚度归气溶胶。公式本身不难难点在V0怎么取、气压怎么修、臭氧要不要扣以及怎么判断反演结果是稳定可信的。3.2 用Langley定标求V0别只依赖出厂值V0按理可以由仪器出厂定标得到但滤光片老化、窗口污染、探测器温漂都会让它变化。现场最可靠的做法是Langley定标在气溶胶光学厚度稳定且太阳高度跨越较大的一段晴空时间里把ln(V)对大气质量数做线性回归截距就是ln(V0)斜率近似负的总光学厚度。Langley定标对大气稳定性要求极高。我一般选择清晨时段太阳天顶角从85度逐渐变化到55度左右这段时间气溶胶和瑞利散射变化最小。如果天顶角跨度太小比如只在70到75度之间回归斜率会被噪声主导V0截距就不可信。筛选条件要同时看天顶角跨度和回归残差。def airmass_kasten_young(sza_deg): Kasten-Young大气质量数公式80度内稳定85度以上不建议用。 z np.radians(sza_deg) if sza_deg 85.0: return np.nan return 1.0 / (np.cos(z) 0.50572 * (96.07995 - sza_deg) ** -1.6364) sel qc[(qc.flag_invalid False) (qc.sza_deg 55) (qc.sza_deg 85)] x sel[sza_deg].apply(airmass_kasten_young) y np.log(sel[dn_440]) k, intercept np.polyfit(x, y, 1) v0_440 np.exp(intercept)这段代码对每个通道要分别跑一遍得到各自的V0。判断定标是否可靠第一步看回归决定系数R²到0.999级别才说明整段大气稳定第二步看截距连续几天上午复测的V0相对偏差如果超过0.5%就要考虑清洁镜筒、检查滤光片污染不能强行用旧定标。还有一点Langley回归时如果云没有滤干净残差会明显偏大截距被云拉偏的幅度足以让全部AOD产生系统性平移。3.3 逐通道AOD计算把公式落成能跑的Python代码拿到V0之后AOD计算就是逐条记录逐通道套公式。我在代码里把瑞利光学厚度独立成一个函数气压作为输入参数而不是写死常量。瑞利光学厚度与波长成λ的负3.084次方关系所以440 nm通道的瑞利项明显大于870 nm这个波长依赖是AOD谱型正确的前提。def rayleigh_tau(wl_nm, pressure_hpa): 瑞利散射光学厚度波长单位nm气压单位hPa。 lam_um wl_nm / 1000.0 return (pressure_hpa / 1013.25) * 0.008735 * lam_um ** -3.084 def aod_from_dn(dn, v0, wl_nm, pressure_hpa, sza_deg, d_fac): 单通道AOD反演。 m airmass_kasten_young(sza_deg) if np.isnan(m) or dn 0: return np.nan tau_r rayleigh_tau(wl_nm, pressure_hpa) numerator (np.log(v0) 2.0 * np.log(d_fac) - np.log(dn) - tau_r * m) return numerator / m公式里2×ln(d)就是日地距离修正不要漏掉。如果不做这项修正北半球夏季和冬季的AOD会出现0.02到0.03量级的伪季节变化。气压参数强烈建议用站点实时气压而不是用海平面标准气压1013.25 hPa代替。高原站点气压可能只有700多百帕瑞利项直接偏差百分之二三十换算到AOD是0.005以上的系统性误差对做趋势分析来说完全不可接受。3.4 输入参数与常见误用的边界参数推荐取值影响气压站点实时气压瑞利光学厚度线性变化V0Langley定标结果整体平移所有AOD日地距离d按年积日计算产生季节伪信号必须修正大气质量数Kasten-Young公式80度以上误差快速放大臭氧气候值或卫生产品只影响500、675 nm约0.002~0.01实际处理中常见的误用有两种。一种是用太阳天顶角的余割函数sec(θz)代替Kasten-Young公式在60度以内两者差别不大但70度以后差距随角度迅速拉大Langley定标和AOD反演都会受影响。另一种是把出厂V0当作永久常量连续用几个月不重新定标。滤光片和窗口的污染是渐变的旧V0会造成AOD缓慢漂移这种漂移很难从单天数据里发现但时间序列一拉长就很刺眼。4. 水汽wv反演936nm吸收通道的透过率换算与质量数修正4.1 936通道为什么能反演水汽936 nm落在水汽的近红外强吸收带内水汽分子对936 nm辐射的吸收远大于气溶胶散射的影响。把瑞利散射和气溶胶的贡献扣除后剩下的透过率衰减主要就是水汽吸收造成的。水汽吸收通道的透过率T_w与垂直柱含水量wv之间存在强相关但又不是严格的比尔-朗伯特关系因为滤光片有一定带宽吸收带内包含多条强弱不等的吸收线。工程上最常用的反演方案是多项式法先由实测936通道DN值算出水汽斜程透过率T_w取对数后用一个三次多项式映射到水汽柱含量。多项式系数来自辐射传输模型或者仪器出厂水汽定标文件不是自己随便回归的这是整个水汽反演里最容易被轻视的点。4.2 气溶胶和瑞利在936nm怎么扣除用870和1020插值936通道没有独立的气溶胶测量必须从临近的非吸收通道推算气溶胶光学厚度。我一般用870和1020两个通道的AOD做Angstrom插值不用440。936与870只差66 nm与1020差84 nm插值误差很小如果把440拉进来细粒子模态会主导Angstrom指数反而把936的气溶胶扣除带偏。先由870、1020的AOD计算Angstrom指数α再外推到936。# 假设aod_870、aod_1020已经由第三章代码逐通道算好 alpha np.log(aod_870 / aod_1020) / np.log(1020.0 / 870.0) aod_936 aod_870 * (936.0 / 870.0) ** (-alpha) tau_r_936 rayleigh_tau(936.0, pressure_hpa) m_936 airmass_kasten_young(sza_deg) t_total dn_936 / (v0_936 * d_fac ** 2) t_w t_total * np.exp(m_936 * (tau_r_936 aod_936))代码里的t_total是936通道总的大气透过率exp因子把瑞利散射和气溶胶透过率除回去剩下的t_w就是水汽斜程透过率。这里v0_936必须来自936通道自己的Langley定标不能用其他通道的V0代替。936通道位于吸收带内Langley定标时不能直接使用需要在晴空干燥条件下用扣除水汽后的截距来定标这一点后面避坑章节还会展开。4.3 从透过率到水汽柱wv多项式回归与质量数修正得到t_w后水汽反演的关键是把ln(t_w)映射到wv。常见的做法是使用如下三次多项式wv a0 a1×ln(t_w) a2×ln(t_w)² a3×ln(t_w)³系数a0到a3由传输模型对一系列水汽柱含量和大气质量数组合模拟后拟合得到。模型模拟时使用的光谱响应函数、大气廓线和温度结构不同系数也会有差异。如果没有仪器自带的水汽定标文件就要用当地大气廓线重新模拟标定这点不能省。from numpy.polynomial import polynomial as P coefs [a0, a1, a2, a3] # 来自传输模型拟合单位对应cm ln_tw np.log(t_w) wv_slant P.polyval(ln_tw, coefs) # 斜程水汽柱含量cm wv_vertical wv_slant / m_936 # 垂直水汽柱含量cm这里为什么要除以m_936取决于系数拟合时用的透过率是怎么构造的。如果模型拟合时直接对斜程透过率与垂直水汽柱含量回归那么系数内部已经隐含了质量数关系输出就是垂直柱含量如果模型是先算出斜程水汽柱量再回归那么输出需要除以质量数才能得到垂直柱量。我的建议是处理历史数据时先按上面代码同时保留wv_slant和wv_vertical两个字段等确认定标文件形式后再决定用哪一个。正常大气条件下整层水汽柱含量在0.1到6 cm之间湿润地区夏季可能接近6 cm干燥高海拔地区冬季可能在0.2 cm以下。如果反演结果出现负值或超过8 cm不要怀疑大气先回去查936通道的V0定标和气溶胶扣除这两处是wv异常的最常见来源。5. CE318数据反演避坑五个让我来回返工的现场问题5.1 现象Langley定标截距每天都不稳定AOD整体平移连续几天做Langley定标每天求出的V0截距差异超过1%对应的AOD整体平移0.01到0.02。排查发现不是仪器故障而是定标时段里有薄卷云混入云没有被预筛剔除导致回归截距被拉偏。Langley回归对大气稳定极度敏感一段数据里只要混入少量变薄的云斜率就会失真。解决方法是把定标样本的筛选标准收紧只保留天顶角跨度足够大的连续晴空段回归决定系数要接近0.999同一时段用相邻通道的AOD时间序列做方差检查AOD标准差偏大的点直接剔除。如果连续几天截距仍然不稳定优先怀疑镜筒窗口污染清洁后重新定标而不是强行调整V0去拟合数据。5.2 现象上午AOD偏高湿度一高更明显有段时间上午的数据AOD系统性偏高0.02到0.03中午以后恢复正常。最初怀疑是Langley定标问题重新定标后依旧。后来发现是日出后镜筒窗口表面凝结了薄露水露水层对太阳辐射产生额外衰减DN值偏低反演出来的AOD就虚高。湿度大的站点这个问题尤其明显。解决思路不是删除整个上午的数据而是增加湿度判据相对湿度接近露点或者超过90%时给对应记录打上可疑标记不进入平均值和趋势统计。镜筒加装小型吹扫装置能明显缓解露水问题但不能完全消除数据质控里还是要保留湿度筛选这一层。5.3 现象936nm反演wv比探空偏大20%某次对比探空数据CE318反演的水汽柱含量系统性偏高。逐项排查后发现问题出在936通道的气溶胶扣除当时用的是440和870插值到936440通道AOD在细粒子污染时段偏高导致插值得到的aod_936偏大瑞利和气溶胶扣除过量t_w被高估最终wv偏大。解决方法是切换到870和1020两个近红外通道插值因为936靠近它们插值距离短、误差小。同时检查了936通道V0的定标日期发现V0文件已经用了半年重新做干燥晴空定标后偏差明显缩小。从此以后我的水汽反演链路默认使用870-1020插值不再把440拉进936的气溶胶扣除。5.4 现象500nm通道出现季节性偏高500 nm通道的AOD在春季总是比同期的周边站点高0.005左右而440和870通道正常。这个量级不大但做长期趋势分析时会出现明显的季节伪信号。原因是不正确处理500 nm附近的臭氧吸收。春季中高纬度臭氧柱浓度较高臭氧在500 nm附近有一定吸收如果不扣除臭氧光学厚度AOD就会被系统性高估。解决方法是把臭氧项加进反演公式臭氧柱量可以使用气候平均产品或者当日卫星臭氧柱数据。臭氧对500、675 nm的影响量级在0.002到0.01对440、870、1020几乎可以忽略。这不是一个漂亮的大问题但正因为它量级小最容易让人忽略最后反映在季节曲线上成为难以解释的假信号。5.5 现象数据出现连续重复值但时间戳正常导出的文件里某一条通道的DN值连续几十条完全一样时间戳却正常递增。这种现象不像云也不像露水用质控规则很难滤掉。排查发现是滤光片轮卡在某个位置每次扫描都重复测量同一个滤光片其他通道数据则正常。如果不检查重复值这些假数据会像真的观测一样进入统计结果。解决方法是给每一组扫描加唯一标识检查一轮扫描内六个通道的DN值是否都在合理范围内变化任一通道DN值与上一条完全相同时整组记录标记为不可信。这个过滤步骤放在预筛阶段代码简单但很有效。我通常用duplicated方法按440和870两列同时判重发现重复比例超过1%就要检查仪器硬件状态。6. 反演结果验证技巧不依赖外部数据的三个自检习惯拿到自己算出的AOD和wv序列后不要急着跟卫星或模式对比先做三个不依赖外部数据的自检。第一个是Angstrom谱型检查。用440和870的AOD计算Angstrom指数α正常气溶胶条件下α在0到2.5之间城市细粒子污染天可能接近1.5以上沙尘天接近0甚至为负。如果α标准化为负数或大于3说明某个通道的V0或云滤除出了问题。alpha_440_870 np.log(aod_870 / aod_440) / np.log(440.0 / 870.0) if not (0.0 alpha_440_870 2.5): print(谱型异常检查V0和云过滤)第二个是相邻通道差值检查。AOD在440到870 nm之间应该是平滑变化的任意两个相邻通道的AOD差不应出现剧烈跳变。特别是936通道所在位置反演出的aod_936应该与870、1020通道保持同一谱型如果aod_936明显偏离水汽反演的结果也不可信。第三个是时间连续性检查。晴空大气AOD在分钟时间尺度上不应有锯齿状跳变把一天的AOD序列画出来看到离散的尖峰基本就是云没滤干净看到整段台阶就是定标变化。与周边同类型站点比对时注意比较同一时刻、同一波长的数据。两台经过独立定标的仪器AOD 440 nm日均值差异在0.01以内属于正常范围wv相差10%以内可以接受。如果差异明显偏大先查两边的V0定标日期和气压输入再查海拔修正。这些自检都过了反演结果才算真正能用。这些年我养成的习惯是每次Langley定标后都存一张回归图截距变化超过0.5%就停下来查原因而不是继续往下跑数据。定标漂移是所有AOD和wv反演问题的源头管住定标就等于管住了后续所有产品的可靠性。希望帮到你。本文还有配套的精品资源点击获取