简介:针对CE318太阳光度计观测数据的专用C++源码包,面向大气科学、遥感反演与气溶胶研究领域的科研人员和学生,解决从仪器原始数据到AOD与水汽含量WV反演的编程实现难题。资源共5个C++源文件,涵盖CE318数据读取、AOD计算、水汽含量反演、Angstrom参数拟合与定标处理等模块,压缩包整体仅9KB,代码精简且便于逐段分析。目前已有1080人学习,透过这份源码可了解完整的数据处理流程,包括算法函数用法、波长指数拟合和大气参数计算细节;对于需要快速上手CE318数据处理或开展相关算法验证的读者,是值得参考的轻量级代码示例。 做CE318数据处理这行有几年了,手头积累最多的工作就是把太阳光度计的原始计数值变成AOD和气柱水汽总量WV。很多刚接触遥感大气产品的人,拿到一台CE318,第一反应是“仪器不是直接给AOD吗?”,实际上仪器给的是电压或者DN值,AOD和WV全靠后期反演。这篇就把我从数据清洗到产品输出的完整流程写出来,尤其是WV那部分,很多文档语焉不详,我尽量用自己的实际经验补全。
先说清楚这东西能干什么。CE318(也叫POM-02或CE318N,不同型号略有差异)是法国Cimel公司生产的全自动太阳光度计,全球AERONET地基观测网的主力设备。它通过追踪太阳,在多个波段测量直射太阳辐射,进而反演大气气溶胶光学厚度AOD和大气水汽柱含量WV。AOD是气溶胶消光在整层大气的积分,单位无量纲,通常在0.05到1.0之间,沙尘或污染天能到2以上;WV单位是g/cm²或cm,代表单位截面积上整层大气水汽的质量。这两个参数是大气校正、气溶胶模式、辐射传输模拟的基础输入,做遥感反演、环境监测、气候评估的都绕不开。
这篇文章适合谁?如果你是刚拿到CE318数据的研一学生、准备做AERONET站点数据质量控制的观测员,或者想用太阳光度计做便携式大气校正的工程人员,这篇都能用得上。我尽量把算法原理讲透,同时给出可操作的处理路径。
1. 项目概述与技术路线
1.1 CE318观测与AOD/WV反演的基本逻辑
CE318的观测原理是老牌的太阳直射辐射测量。仪器内部有光电探测器,通过滤光片选择特定波段,测量到达地面的太阳直接辐射照度。原始输出是电压值(或数字DN值),与大气透过率之间存在确定的物理关系。
我把核心关系写成一个直观公式:
V(λ) = V0(λ) × R² × exp(-m × τ(λ)) × exp(-m_w × α_w × wv)
其中V是地面测量电压,V0是大气层顶的定标电压(也就是仪器在无大气情况下对同一太阳源的理论响应),R是日地距离修正因子,m是大气质量数,τ(λ)是波段总光学厚度,m_w是水汽有效大气质量数,α_w是水汽吸收系数,wv是水汽柱含量。
这个公式实际上就是把气溶胶和水汽的衰减从总大气衰减中分离开来。AOD反演走的是非水汽吸收通道,比如440nm、675nm、870nm、1020nm,在这些波段水汽吸收极弱,可以忽略;WV反演走的是936nm或940nm附近的水汽弱吸收带,利用水汽吸收对总透过率的影响来反推水汽柱含量。
1.2 我的技术路线总览
处理CE318数据,我一般按这么几条线走:
- 数据质量控制和预处理:检查原始电压信号是否有云遮挡、太阳跟踪是否正常、定标是否漂移。
- 计算太阳天顶角、日地距离修正、大气质量数。
- 对非吸收通道做Langley定标或者利用已有定标系数,反演AOD,并修正Rayleigh散射、臭氧吸收、NO₂吸收等。
- 对水汽通道,用改进型Langley法或比值法提取WV。
- 输出时间序列产品,做云检测和异常值剔除,生成最终Level 2数据。
这里我特别强调一点:不要把仪器自带软件给的结果当成最终产品。CE318配套软件(比如ASTPWIN)能输出AOD,但它的质量控制没有AERONET官方那么严,特别在云污染和数据连续性上容易出问题。自己动手处理一遍,既能理解参数意义,也能在论文里把误差来源写清楚。
2. 数据准备与质量控制
2.1 CE318原始数据的主要字段
CE318的数据文件常见有几种格式。以Level 1.0直射太阳数据为例,每个观测记录通常包括:
- 时间戳,一般是UTC。
- 太阳天顶角SZA和方位角AZM。
- 各通道的原始电压值或DN值,单位mV或counts。
- 仪器内部温度、湿度、跟踪误差等辅助信息。
拿到数据后我第一步不是直接算AOD,而是画一个全天的电压时间序列。太阳直射光强随时间的变化应当是平滑的,如果电压出现突然掉下去又恢复,那基本是云遮挡;如果电压整体漂移得很厉害,要检查是不是仪器内部湿度或温度异常。
2.2 云检测与无效数据剔除
云是太阳光度计最大的敌人。云遮挡时,测量到的直射辐射偏低,反演出来的AOD会异常偏高;如果云的边缘效应,甚至可能出现AOD跳变。
我现在常用的云检测方法有三种:
- 阈值法:单次观测的AOD如果与前一个观测值相对变化超过20%,且不是沙尘或火山灰事件,就判断为可疑。
- 稳定性法:同一站点在短时间内多次观测,AOD的标准差过大则剔除。
- DIMM法(差分图像运动监测):这个需要额外硬件,日常站点一般不做。
操作上,我建议先把原始电压数据按照每30秒或者每3分钟的观测频次画出曲线,用肉眼剔除明显云污染时段,然后再用代码做阈值筛查。AERONET的Level 2产品是经过严格云检测的,如果你只是用站点Level 1.5数据,等于是把质量控制的责任接了过来,不能偷懒。
2.3 定标系数的核对与校正
CE318的定标系数V0决定了AOD反演的绝对精度。V0的获取通常有三种途径:
- 实验室定标:在国家计量院或仪器厂商用标准灯标定,精度最高,但仪器需要返厂。
- 高山Langley定标:在海拔高、气溶胶极少的站点,利用太阳高度角连续变化做Langley图。
- 传递定标:与一台已定标的参考仪器在同一时间和地点进行对比观测。
我的经验是:每半年到一年要进行一次定标,定标系数变化超过5%就要特别注意。尤其是440nm和500nm通道,硅探测器老化敏感,最容易出现漂移。在实际数据处理时,如果发现多天的AOD基线整体抬升或下降,而气象上没有相应的污染或沙尘过程,十有八九是定标漂移,这时需要对V0做时间插值校正。
3. AOD反演原理与实操拆解
3.1 朗伯-比尔定律与大气质量数计算
对于非水汽吸收通道,朗伯-比尔定律可以写成:
τ_a(λ) = (1/m) × ln(V0(λ)×R² / V(λ)) - τ_R(λ) - τ_O3(λ) - τ_NO2(λ)
也就是说,总光学厚度中扣除分子Rayleigh散射、臭氧吸收、NO₂吸收后,剩下的就是气溶胶光学厚度AOD。
大气质量数m并不简单等于1/cos(SZA)。对平面平行大气,m=sec(θ),但考虑到地球曲率和大气折射,在太阳高度角很低时有偏差。我常用Kasten-Young公式:
m = 1 / (cos(SZA) + 0.50572 × (96.07995 - SZA)^(-1.6364))
SZA单位是度。这个公式在SZA达到85度时依然有较好效果。实际处理时我通常把SZA大于80度的数据直接剔除,因为此时大气质量数过大,反演误差迅速放大。
日地距离修正R也很关键。R = (rm/r)²,可以用近似公式计算:
R = 1 - 0.01674 × cos(0.0172×(day_of_year - 4))
其中day_of_year是年积日。如果想要更高精度,可以用NASA的ephemeris算法,但日常处理这个简单公式足够。
3.2 多通道AOD计算与分子散射修正
在我处理的CE318标准配置里,AOD核心通道是440nm、500nm、675nm、870nm、1020nm,有些型号还有1640nm。每个通道的Rayleigh散射光学厚度用公式计算:
τ_R(λ) = 0.008569 × λ^(-4) × (1+0.0113×λ^(-2)+0.00013×λ^(-4)) × (P/P0)
其中λ单位是微米,P是站点气压,P0是标准大气压1013.25hPa。这个公式在不同文献里略有差异,但一致性都很好。
臭氧吸收修正主要影响440nm和500nm通道,NO₂吸收影响440nm左右。这两项需要输入臭氧和NO₂柱总量,可以从卫星产品或再分析资料获取。我的经验是:在气溶胶反演精度要求不高的日常监测中,NO₂修正可以忽略,但臭氧修正不能省,尤其在热带地区臭氧柱浓度高,影响能达到0.01到0.02。
3.3 Angstrom指数与谱分布特征
得到多个波段的AOD后,我最常做的事情是拟合Angstrom公式:
τ_a(λ) = β × λ^(-α)
对两侧取对数,用440nm和870nm两个波段算:
α = -ln(τ_440 / τ_870) / ln(440/870)
这个α就是Angstrom指数,反映了气溶胶粒子谱分布的斜率。α越大说明细粒子占比越高,污染型气溶胶通常α在1.2到2.0;α越小说明粗粒子占比高,沙尘过程α常在0到0.5。这个参数对判断气溶胶类型和做大气校正很关键,我在处理数据时都会同步输出。
注意:如果反演出的某波段AOD是负值,不要直接丢掉,先检查是不是定标系数偏大或者云污染。负值往往意味着测量电压高于“大气层顶电压”,这在物理上不可能,除非定标失效。
4. 水汽柱含量WV反演详解
4.1 水汽通道的选择依据
CE318上用于水汽反演的通道一般是936nm,部分型号是940nm,带宽约10nm。这个波段位于水汽的弱吸收带,气溶胶消光和水汽吸收同时存在,但通过合适的数学变换可以把气溶胶影响剥离。
之所以不选择更强吸收的820nm或更强的水汽带,是因为如果水汽吸收太强,信号会迅速饱和,动态范围太小;如果太弱,反演灵敏度又不够。936nm这个位置是“弱吸收但灵敏”的折中。AERONET的WV反演用的就是这一通道。
4.2 改进型Langley法标定
水汽通道的Langley定标和普通通道不一样。普通通道的Langley图是ln(V)对m作图,截距是ln(V0)。但水汽通道的透过率同时受水汽吸收影响,直接做线性Langley会因为水汽变化而失败。
改进型方法的核心是:把水汽吸收项参数化为传输函数,通常用二参数的指数形式:
T_w = exp(-a × (m_w × wv)^b)
其中a和b是通道的带通特性决定的经验系数,m_w是水汽有效大气质量数。这样Langley图要处理的对象变成了:
ln(V0) - a×(m_w × wv)^b
实际上,AERONET的处理思路是把扣除气溶胶和Rayleigh散射后的水汽通道透过率,与一个已知水汽柱总量做比值标定。如果你没有条件做高山标定,最稳妥的做法是使用AERONET官方发布的定标系数,或找一台经过校准的参考仪器做传递标定。
我自己的具体做法是:从AERONET站点的Level 2数据里下载同站点同时间的WV值,利用它来反推自己仪器的等效吸收系数。这个做法在长期对比观测中很好用,能显著改善WV反演的绝对精度。
4.3 WV计算过程与结果解读
在实际反演时,我通常先对非吸收通道做AOD反演,然后通过Angstrom关系插值得到936nm处的气溶胶光学厚度τ_a,936。接着计算水汽通道的大气总透过率:
T_total = V(936) / (V0_936 × R²)
扣除Rayleigh散射、气溶胶消光后,得到水汽透过率:
T_w = T_total / exp(-m×τ_R,936) / exp(-m×τ_a,936)
然后利用水汽传输函数反算WV:
wv = [ -ln(T_w) / a ]^(1/b) / m_w
这里的m_w可取大致为1/cos(SZA)的修正形式,在SZA小于70度时与m接近。
我在实操中遇到过很多次WV反演值跳变的问题。常见原因是936nm通道的信号噪声大,尤其是雨后或高湿度环境下,水汽透过率很容易被云污染干扰。检测手段很简单:如果WV时间序列出现单点极端值,而相邻时刻正常,基本是云或视场内湿度不均一造成。这种情况直接剔除,不要用平滑滤波掩盖。
5. 数据处理实战:软件与脚本
5.1 官方软件ASTPWIN的使用要点
CE318配套的ASTPWIN软件能完成大部分基础处理,它把Langley定标、AOD计算和WV反演做成了图形界面。但说实话,我对这个软件又爱又恨。爱的是它开箱即用,适合现场快速查看;恨的是它的批量处理能力弱,质量控制也不透明。
如果你要用ASTPWIN,我建议关注几个关键设置:
- 定标文件路径要正确,V0和不确定度要同步更新。
- 数据过滤选项里,勾选太阳天顶角限制,通常设为小于80度。
- 输出格式选择ASCII,方便后续处理。
用软件算完,一定要把中间量比如各通道电压、太阳天顶角、气压一起导出。只导出最终AOD的话,后面排查问题会很痛苦。
5.2 Python批处理思路
我更推荐的路线是把数据处理写成Python脚本。这里给出一个简化版本的伪代码流程,主要是帮读者建立整体框架:
import pandas as pd import numpy as np # 读取CE318原始数据 df = pd.read_csv('ce318_raw.csv', parse_dates=['time']) # 计算太阳天顶角、日地距离、大气质量 df['sza'] = calculate_sza(df['time'], lon, lat) df['R'] = earth_sun_distance(df['time'].dt.dayofyear) df['m'] = kasten_young_airmass(df['sza']) # 通过朗伯比尔定律计算各波段AOD for wavelength in [440, 500, 675, 870, 1020]: df[f'aod_{wavelength}'] = ( np.log(df[f'v0_{wavelength}'] * df['R']**2 / df[f'v_{wavelength}']) / df['m'] - rayleigh_od(wavelength, pressure) - ozone_od(wavelength, ozone) - no2_od(wavelength, no2) ) # 云检测与异常剔除 df['flag'] = cloud_check(df, threshold=0.2) # 计算Angstrom指数 df['angstrom'] = -np.log(df['aod_440'] / df['aod_870']) / np.log(440/870) # 水汽通道反演WV df['wv'] = inversion_wv(df['v_936'], df['v0_936'], df['m'], df['aod_870'], pressure)实际脚本里,最需要注意的是把公式中的单位统一。波长录成纳米还是微米,气压录成hPa还是Pa,稍不留神就会差出好几个数量级。我的习惯是全部转成国际单位,波长用微米,气压用hPa,电压用mV。
5.3 精度验证与产品输出
处理完之后不能直接写结论,还要做精度验证。最常用的方法是与同站的AERONET Level 2产品对比。AERONET的AOD不确定度大约在0.01到0.02,WV在10%左右。如果你的反演结果和AERONET差得不多,说明仪器和流程靠谱;如果系统性偏差,优先怀疑定标。
输出产品时,我一般保存成NetCDF或CSV,包含时间、SZA、AOD各通道、Angstrom指数、WV、质量控制标志。这样不管后面是画图、做统计还是做辐射传输输入,都很方便。文件命名里建议加上站点名、日期、数据级别(L1.5或L2),这个习惯帮我省了很多事。
6. 常见问题与排查技巧
6.1 AOD出现负值或异常偏高
这是最频繁的问题。AOD为负,基本指向定标系数V0偏大或者大气质量数偏大。排查步骤:
- 检查定标文件是否过期,V0漂移了多少。
- 检查太阳天顶角计算是否正确,时区、经度纬度是否填对。
- 检查是不是混入了云污染,电压偏低会导致AOD偏高。
AOD异常偏高但序列相对平滑,通常是定标问题;AOD突然跳高又回落,大概率是云。把这两类区分开,问题就解决了一大半。
6.2 WV反演失败或结果离谱
WV反演结果离谱,先不要急着调模型,看看936nm通道电压值是否正常。如果电压波动剧烈,可能是水汽通道探测器故障或滤光片受潮。还有一个高频坑是:气溶胶通道和水汽通道的时间同步问题,CE318不同通道是依次测量的,如果此时太阳被云挡住又很快露出来,水汽通道的测量值可能对应到不同的云状态,导致WV异常。
我的处理经验是给WV反演增加一道约束:只有当936nm测量前后的1分钟内的非吸收通道AOD变化小于0.02时,才接受该WV值。这个条件在AERONET官方处理里也有类似逻辑。
6.3 仪器跟踪异常与数据断层
CE318靠太阳跟踪器工作,但阴天时跟踪器找不到太阳,数据会出现大段缺失。这是正常的,不需要补插。如果你发现数据在晴天依然断层,要先检查跟踪器是否被遮挡、电机是否卡顿、GPS授时是否漂移。GPS时间漂移会导致太阳位置计算错误,进而导致测量视场偏离太阳,数据看起来像云污染,其实是跟踪失败。
7. 一点个人体会
最后说点实在的。CE318的数据处理看似是“点击软件、出图”的简单事,真正难点在于把物理过程理解透。AOD和WV不是仪器直接测出来的物理量,而是通过模型从电压信号中反演出来的。定标、云检测、大气质量数、水汽吸收参数,每一步都带着误差,要清楚自己产品的误差来源,才能判断数据能不能用。
我在实际工作中最深的体会是,好数据不是处理出来的,而是观测出来的。日常仪器的清洁、干燥剂更换、窗口玻璃的检查,比后期算法更重要。再好的反演算法也救不回被污染的光学镜头数据。
如果后面要把这套流程做得更完善,我建议在云检测里引入更多天空辐射观测数据,或者直接对接全天空成像仪,彻底把云污染时段剔除干净。这样处理出来的AOD和WV时间序列,无论是做气候分析还是卫星验证,都会更让人放心。
本文还有配套的精品资源,点击获取