企业官网建设流程全解析

简介：面向大气科学、环境遥感领域的研究者与学习者，这份压缩包以太阳光度计原始观测数据为输入，提供了反演气溶胶光学厚度与水汽含量的完整编程实现。资源共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数据时默认的通道分工，新手可以直接照着建字段。

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'], keep=False).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 输入参数与常见误用的边界

实际处理中常见的误用有两种。一种是用太阳天顶角的余割函数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插值到936，440通道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反演问题的源头，管住定标就等于管住了后续所有产品的可靠性。希望帮到你。

本文还有配套的精品资源，点击获取