1. 从数据到决策:为什么我们需要亲手计算AQI?

大家好,我是老张,一个在数据分析和环境科学交叉领域摸爬滚打了十来年的“老码农”。这些年,我见过太多朋友,无论是做数据分析的同行,还是环境专业的学生,一提到空气质量指数(AQI),第一反应就是去查天气App或者环保部门的官网。这当然没错,但如果你想真正理解数据背后的故事,甚至想自己动手做一些预测、分析或者参加个数学建模比赛,那光会看现成的数字是远远不够的。

这就好比你想成为一个美食家,不能只满足于品尝成品,还得知道这道菜是怎么从食材一步步做出来的。AQI就是这个“成品”,而六种污染物(PM2.5、PM10、SO₂、NO₂、CO、O₃)的浓度就是“原始食材”。官方的计算逻辑,就藏在国家发布的《环境空气质量指数(AQI)技术规定(试行)》(HJ 633-2012)这份文档里。这份文档就是我们的“菜谱”。

自己动手算AQI,好处太多了。首先,你能彻底搞懂“轻度污染”和“中度污染”之间那个关键的数值界限是怎么来的,理解不同污染物对最终指数的“贡献”权重。其次,当你拿到一批原始的监测站小时数据时,你可以独立验证数据的可靠性,甚至能自己搭建一个简单的空气质量预警系统。最后,对于参加数学建模、数据分析竞赛的同学来说,这更是一个展示你数据处理和建模能力的绝佳机会——别人还在用现成的AQI数据,你已经能从原始浓度开始构建整个分析链条了,这差距一下子就拉开了。

所以,这篇文章,我就打算用最接地气的Python代码,手把手带你走一遍这个“从食材到菜肴”的全过程。我们不谈空泛的理论,就对着那份“国标菜谱”,一行行代码实现它。我会分享我实际写代码时踩过的坑,比如数据边界怎么处理、小数点位怎么四舍五入才符合规范,以及算出来的结果怎么用图表变得一目了然。目标就一个:让你看完就能自己动手,把一坨枯燥的浓度数据,变成有决策价值的空气质量洞察。

2. 拆解“国标菜谱”:AQI的计算公式到底在说什么?

在动手写代码之前,我们必须先把那份“国标菜谱”——HJ 633-2012看懂。别被技术规定的名头吓到,它的核心逻辑其实非常直观,我们可以把它想象成一个“分段评分器”。

这个评分器针对每一种污染物(比如PM2.5)都有一张对应的“评分卡”。这张卡把污染物的浓度范围划分成了几个区间,每个区间对应一个固定的AQI子指数(IAQI)范围。举个例子,PM2.5的24小时平均浓度,0到35微克/立方米对应IAQI 0-50(优),35到75对应IAQI 50-100(良),以此类推,直到浓度大于500对应IAQI 300-500(严重污染)。这个对应关系是标准里明确给定的表格,我们的代码本质上就是在查这个表。

那么,如果某个PM2.5浓度值正好落在两个区间中间,比如50微克/立方米,该怎么评分呢?这里就用到了线性插值法。公式长这样:

IAQI_p = [(IAQI_high - IAQI_low) / (C_high - C_low)] * (C_p - C_low) + IAQI_low

看起来有点复杂?我们用人话翻译一下:假设PM2.5浓度50,落在了“35-75”这个浓度区间,这个区间对应的IAQI是“50-100”。那么:

  • C_low=35, C_high=75
  • IAQI_low=50, IAQI_high=100
  • 我们测到的浓度 C_p=50

套进公式就是:IAQI_p = [(100-50)/(75-35)] * (50-35) + 50 = (50/40)*15 + 50 = 1.25*15 + 50 = 68.75

最后,对这个结果进行四舍五入取整,得到这个PM2.5浓度对应的IAQI就是69。看,是不是很像初中数学里的比例计算?这就是AQI计算最核心的数学原理。

这里有一个非常关键的细节,也是我早期写代码时栽过跟头的地方:AQI最终的值,不是各种污染物IAQI的平均值,而是它们的最大值! 也就是说,我们会分别计算出PM2.5、PM10、SO₂等六种污染物各自的IAQI,然后取其中最大的那个数,作为最终的AQI值。同时,这个最大值对应的污染物,就被称为“首要污染物”。这很好理解,空气质量是由最差的那个“短板”决定的。你的PM2.5可能是“良”,但如果臭氧突然爆表到了“重度污染”,那么整体的AQI就是“重度污染”,首要污染物就是臭氧。

理解了这些,我们再看标准里的两个主要指数:1小时AQI和24小时AQI。它们的区别主要在于两点:一是 averaging time(平均时间),1小时AQI用的是污染物最近1小时的平均浓度,反应瞬时状况;24小时AQI用的是过去24小时滑动平均浓度,反应长期暴露水平。二是部分污染物的浓度限值表不同,比如O₃(臭氧),它有1小时和8小时两种评价标准。在实际计算中,1小时AQI通常直接借用24小时AQI的IAQI限值表来计算PM2.5和PM10(这是标准中明确说明的),而SO₂、NO₂、CO、O₃则使用1小时特定的限值表。我们在代码里准备两套不同的“评分卡”(专业叫法是Idata),根据需求切换使用就行了。

3. 实战准备:搭建你的Python数据分析环境

工欲善其事,必先利其器。咱们这个实战项目不需要多么复杂的环境,但几个核心的Python库必不可少。我强烈建议你使用Anaconda来管理环境,它能避免很多包依赖的麻烦。

首先,打开你的终端(Windows用CMD或PowerShell,Mac/Linux用Terminal),创建一个新的虚拟环境,这样不会和你电脑上其他项目冲突。命令很简单:

conda create -n aqi_analysis python=3.9
conda activate aqi_analysis

接下来,安装我们需要的“四大金刚”:

pip install pandas numpy matplotlib seaborn
  • pandas:数据处理的瑞士军刀。我们读取JSON、CSV数据,整理表格全靠它。
  • numpy:提供高效的数组运算和数学函数。后面那个线性插值公式用numpy写起来会非常简洁。
  • matplotlib & seaborn:图表绘制的黄金组合。Matplotlib是基础,Seaborn基于它,能让我们用更少的代码画出更美观的统计图表。

数据从哪里来?我们可以从中国环境监测总站等官方平台的公开数据接口获取,或者直接使用一些竞赛、研究公开的数据集。为了演示,我这里准备了一个模拟的air_data.json文件,里面包含了某个监测点连续几天、每小时记录的六项污染物浓度。文件结构大致如下,你可以用文本编辑器自己创建一个:

{
  "site_name": "示范监测点",
  "data": [
    {
      "time": "2023-10-01 00:00:00",
      "PM2.5": 28.5,
      "PM10": 45.2,
      "SO2": 8.1,
      "NO2": 32.4,
      "CO": 0.8,
      "O3": 102.3
    },
    // ... 更多时间点的数据
  ]
}

有了环境和数据,我们的厨房就算准备好了,接下来可以开始按照“菜谱”正式烹饪了。

4. 核心代码逐行解析:手把手实现AQI计算器

现在,我们进入最核心的环节——写代码。我会把代码分成几个逻辑清晰的函数,并逐块讲解。你可以跟着我一起写,也可以先通读理解。

第一步:定义“评分卡” 这是整个计算的地基,必须严格按照国标来。我们定义两个字典,分别存放24小时和1小时AQI计算时,各污染物浓度区间对应的IAQI上限。注意,列表中的浓度值是每个IAQI区间的上限(比如IAQI=50时,PM2.5的浓度上限是35)。

def get_aqi_breakpoints(aqi_type='24h'):
    """
    根据AQI类型(24h或1h)返回污染物浓度限值表。
    参数 aqi_type: '24h' 或 '1h'
    返回: 一个字典,键为污染物名称,值为浓度限值列表。
    """
    # IAQI 对应的标准分段点
    iaqi_points = [0, 50, 100, 150, 200, 300, 400, 500]

    # 24小时平均AQI的浓度限值 (单位: PM2.5, PM10, SO2, NO2为 ug/m3; CO为 mg/m3; O3为 ug/m3)
    breakpoints_24h = {
        'PM2.5': [0, 35, 75, 115, 150, 250, 350, 500],
        'PM10': [0, 50, 150, 250, 350, 420, 500, 600],
        'SO2': [0, 150, 500, 650, 800, 1600, 2100, 2620],
        'NO2': [0, 100, 200, 700, 1200, 2340, 3090, 3840],
        'CO': [0, 5, 10, 35, 60, 90, 120, 150], # 单位 mg/m3
        'O3': [0, 160, 200, 300, 400, 800, 1000, 1200] # 这是8小时平均的O3限值,用于24h AQI
    }

    # 1小时平均AQI的浓度限值 (注意:PM2.5和PM10通常沿用24h标准,O3使用1小时标准)
    breakpoints_1h = {
        'PM2.5': breakpoints_24h['PM2.5'], # 沿用
        'PM10': breakpoints_24h['PM10'],   # 沿用
        'SO2': [0, 150, 500, 650, 800], # 1h SO2限值分段更少
        'NO2': [0, 200, 700, 1200, 2340, 3090, 3840],
        'CO': [0, 10, 35, 60, 90, 120, 150],
        'O3': [0, 100, 160, 215, 265, 800] # 1小时O3限值
    }

    if aqi_type == '24h':
        return iaqi_points, breakpoints_24h
    else: # '1h'
        return iaqi_points, breakpoints_1h

第二步:实现单个污染物的IAQI计算函数 这个函数就是前面讲的“分段评分器”的逻辑实现。它接收一个污染物浓度值、该污染物的浓度限值列表和IAQI分段点列表,返回计算出的IAQI。

def calculate_iaqi(concentration, bp_list, iaqi_list):
    """
    计算单个污染物的IAQI。
    参数 concentration: 污染物浓度值
           bp_list: 该污染物的浓度限值列表
           iaqi_list: IAQI分段点列表
    返回: 计算得到的IAQI值 (整数)
    """
    # 处理浓度低于最低限值或高于最高限值的情况
    if concentration <= bp_list[1]: # 浓度在第一个区间(0-bp_list[1])
        # 按比例计算,因为第一个区间起点是0
        ratio = (concentration - bp_list[0]) / (bp_list[1] - bp_list[0])
        iaqi = ratio * (iaqi_list[1] - iaqi_list[0])
    elif concentration > bp_list[-1]: # 浓度超过表格最大值,直接赋最高值
        iaqi = iaqi_list[-1]
    else:
        # 找到浓度落在哪个区间
        for i in range(1, len(bp_list)):
            if concentration <= bp_list[i]:
                low_bp, high_bp = bp_list[i-1], bp_list[i]
                low_iaqi, high_iaqi = iaqi_list[i-1], iaqi_list[i]
                # 应用线性插值公式
                iaqi = (high_iaqi - low_iaqi) / (high_bp - low_bp) * (concentration - low_bp) + low_iaqi
                break

    # 四舍五入取整,这是国标要求
    return int(round(iaqi))

注意:这里我特别处理了浓度在最低区间和超出最高区间的情况。国标规定,浓度低于最低限值时,按比例计算(从0开始);超过最高限值时,IAQI直接取500。很多网上的示例代码会忽略这个边界处理,导致计算结果不准确。

第三步:计算整体AQI和首要污染物 现在,我们可以对一行数据(一个时间点的六种浓度)计算整体AQI了。

def calculate_aqi_for_row(row, aqi_type='24h'):
    """
    计算单行数据(一个时间点)的AQI和首要污染物。
    参数 row: 一个字典或Pandas Series,包含六项污染物浓度
           aqi_type: '24h' 或 '1h'
    返回: (aqi_value, primary_pollutant) 元组
    """
    iaqi_points, breakpoints = get_aqi_breakpoints(aqi_type)
    pollutants = ['PM2.5', 'PM10', 'SO2', 'NO2', 'CO', 'O3']

    iaqi_dict = {}
    for poll in pollutants:
        conc = row[poll]
        bp = breakpoints[poll]
        iaqi = calculate_iaqi(conc, bp, iaqi_points)
        iaqi_dict[poll] = iaqi

    # AQI是各IAQI中的最大值
    aqi_value = max(iaqi_dict.values())
    # 首要污染物是达到最大IAQI的污染物(可能多个)
    primary = [p for p, i in iaqi_dict.items() if i == aqi_value]

    return aqi_value, primary

第四步:批量处理与数据整合 最后,我们写一个主函数,读取整个数据集,循环计算每个时间点的AQI,并把结果保存到新的DataFrame中。

import pandas as pd

def calculate_aqi_from_file(filepath, aqi_type='24h'):
    """
    从文件读取数据,批量计算AQI。
    参数 filepath: 数据文件路径(支持JSON或CSV)
           aqi_type: '24h' 或 '1h'
    返回: 包含原始数据和计算结果的DataFrame
    """
    # 读取数据
    if filepath.endswith('.json'):
        df = pd.read_json(filepath)
        # 假设数据在嵌套结构中,根据实际情况调整
        # 这里假设数据直接是列表形式
        if 'data' in df.columns:
            df = pd.DataFrame(df['data'].tolist())
    elif filepath.endswith('.csv'):
        df = pd.read_csv(filepath)
    else:
        raise ValueError("仅支持JSON或CSV文件")

    # 确保时间列是datetime类型
    if 'time' in df.columns:
        df['time'] = pd.to_datetime(df['time'])

    # 应用计算函数到每一行
    results = df.apply(lambda row: calculate_aqi_for_row(row, aqi_type), axis=1, result_type='expand')
    df['AQI'] = results[0]
    df['Primary_Pollutants'] = results[1]

    # 根据AQI值添加空气质量等级
    def get_aqi_level(aqi):
        if aqi <= 50:
            return '优'
        elif aqi <= 100:
            return '良'
        elif aqi <= 150:
            return '轻度污染'
        elif aqi <= 200:
            return '中度污染'
        elif aqi <= 300:
            return '重度污染'
        else:
            return '严重污染'

    df['AQI_Level'] = df['AQI'].apply(get_aqi_level)

    return df

# 使用示例
if __name__ == '__main__':
    result_df = calculate_aqi_from_file('air_data.json', aqi_type='1h')
    print(result_df[['time', 'PM2.5', 'AQI', 'AQI_Level', 'Primary_Pollutants']].head())
    # 保存结果
    result_df.to_csv('calculated_aqi_results.csv', index=False)

跑完这段代码,你的原始数据旁边就会多出三列:AQI(具体的指数值)、Primary_Pollutants(首要污染物列表)、AQI_Level(空气质量等级)。至此,从原始数据到AQI决策信息的核心转换就完成了。你可以打开生成的CSV文件检查一下,看看计算出的AQI和首要污染物是否符合你的直观预期。

5. 让数据说话:用可视化洞察空气质量趋势

算出AQI只是第一步,如何让这些数字变得直观、有洞察力,才是数据分析的精华所在。图表是最好的语言。下面我用最常用的matplotlibseaborn来展示几种实用的可视化方法。

第一种:AQI时间序列折线图 这是最基础也最有效的图表,能一眼看出空气质量随时间的变化趋势、污染峰值出现在什么时候。

import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

def plot_aqi_timeline(df, title='AQI时间序列变化'):
    """
    绘制AQI随时间变化的折线图。
    """
    plt.figure(figsize=(14, 6))
    plt.plot(df['time'], df['AQI'], marker='o', linestyle='-', linewidth=1.5, markersize=4, color='steelblue')
    # 添加等级分区背景色
    plt.axhspan(0, 50, facecolor='green', alpha=0.1, label='优')
    plt.axhspan(50, 100, facecolor='yellow', alpha=0.1, label='良')
    plt.axhspan(100, 150, facecolor='orange', alpha=0.1, label='轻度污染')
    plt.axhspan(150, 200, facecolor='red', alpha=0.1, label='中度污染')
    plt.axhspan(200, 300, facecolor='purple', alpha=0.1, label='重度污染')
    plt.axhspan(300, 500, facecolor='maroon', alpha=0.1, label='严重污染')

    plt.title(title, fontsize=16)
    plt.xlabel('时间', fontsize=12)
    plt.ylabel('AQI', fontsize=12)
    plt.grid(True, linestyle='--', alpha=0.6)
    plt.legend(loc='upper right')
    plt.xticks(rotation=45)
    plt.tight_layout()
    plt.show()

# 调用函数
plot_aqi_timeline(result_df, '监测点近三日AQI小时变化趋势')

这张图能立刻告诉你,污染通常发生在一天中的哪个时段(比如早晚高峰),以及是否有持续的污染过程。

第二种:首要污染物堆叠面积图 AQI只告诉我们总体情况,但究竟是哪种污染物在“作祟”?堆叠面积图可以清晰展示不同污染物对IAQI的贡献度随时间的变化。

def plot_primary_pollutant_stacked(df):
    """
    绘制首要污染物的堆叠面积图(简化版,展示各污染物IAQI)。
    注意:这里需要先计算各污染物的IAQI并存储。我们在主流程中稍作修改。
    """
    # 假设我们在计算时,将各污染物的IAQI也存成了单独的列,例如‘IAQI_PM2.5’
    # 这里为了演示,我们重新快速计算一下(实际应用中应在主流程保存)
    pollutants = ['PM2.5', 'PM10', 'SO2', 'NO2', 'CO', 'O3']
    iaqi_points, breakpoints = get_aqi_breakpoints('1h')
    for poll in pollutants:
        df[f'IAQI_{poll}'] = df[poll].apply(lambda x: calculate_iaqi(x, breakpoints[poll], iaqi_points))

    plt.figure(figsize=(14, 7))
    # 准备堆叠数据
    stack_data = df[[f'IAQI_{p}' for p in pollutants]]
    colors = ['lightblue', 'lightgreen', 'gold', 'salmon', 'violet', 'grey']
    plt.stackplot(df['time'], stack_data.T, labels=pollutants, colors=colors, alpha=0.8)

    plt.title('各污染物IAQI贡献度(堆叠面积图)', fontsize=16)
    plt.xlabel('时间', fontsize=12)
    plt.ylabel('IAQI', fontsize=12)
    plt.legend(loc='upper left')
    plt.grid(True, linestyle='--', alpha=0.3)
    plt.xticks(rotation=45)
    plt.tight_layout()
    plt.show()

从这张图,你可以分析出污染类型。比如,如果红色(代表SO2)和金色(代表NO2)的区域在冬季供暖期显著增大,可能意味着燃煤排放是主要来源;如果夏季灰色(O3)区域突出,则可能与光化学污染有关。

第三种:空气质量等级分布饼图 对于一份长时间的数据,我们可能想宏观了解空气质量等级的分布情况:有多少天是优,多少天是良?

def plot_aqi_level_distribution(df):
    """
    绘制空气质量等级分布饼图。
    """
    level_counts = df['AQI_Level'].value_counts()
    # 按优、良、轻度污染...的顺序排序
    order = ['优', '良', '轻度污染', '中度污染', '重度污染', '严重污染']
    level_counts = level_counts.reindex(order, fill_value=0)

    plt.figure(figsize=(8, 8))
    colors = ['green', 'yellowgreen', 'gold', 'orange', 'red', 'darkred']
    patches, texts, autotexts = plt.pie(level_counts.values, labels=level_counts.index, colors=colors,
                                         autopct='%1.1f%%', startangle=90)
    # 美化文本
    for text in texts:
        text.set_fontsize(12)
    for autotext in autotexts:
        autotext.set_color('white')
        autotext.set_fontweight('bold')
    plt.title('空气质量等级分布', fontsize=16)
    plt.axis('equal')  # 保证饼图是圆形
    plt.show()

这张图非常适合用在报告里,给人一个整体印象。你可以轻松地告诉别人:“在过去一个月里,我们监测点有70%的时间空气质量是优良。”

第四种:污染物浓度与AQI的散点图矩阵 想深入探究污染物浓度和最终AQI之间的关系吗?散点图矩阵(Pair Plot)是探索多变量关系的利器。

def plot_pollutant_aqi_correlation(df):
    """
    绘制污染物浓度与AQI的散点图矩阵。
    """
    # 选取我们关心的变量
    plot_df = df[['PM2.5', 'PM10', 'O3', 'AQI', 'AQI_Level']].copy()
    # 使用seaborn的pairplot,按等级着色
    g = sns.pairplot(plot_df, hue='AQI_Level', palette={'优':'green','良':'yellow','轻度污染':'orange','中度污染':'red','重度污染':'purple','严重污染':'brown'},
                     diag_kind='kde', plot_kws={'alpha':0.6, 's':30})
    g.fig.suptitle('污染物浓度与AQI关系矩阵图', y=1.02, fontsize=16)
    plt.tight_layout()
    plt.show()

在这张矩阵图中,对角线是每个变量的分布曲线,非对角线是两两之间的散点图。你可以清晰地看到,PM2.5和AQI的散点图呈明显的正相关,且随着颜色从绿变红(污染加重),点群向右上方移动。这直观地验证了PM2.5是影响AQI的关键因子。同时,你可能会发现O3和PM2.5在某些情况下呈现负相关,这反映了夏季臭氧污染和颗粒物污染往往不同时出现的复杂大气化学过程。

把这些图表组合起来,一份专业的空气质量分析报告就有了坚实的数据可视化基础。你可以根据分析目的,选择最合适的图表来讲述你的数据故事。

6. 避坑指南与进阶思考:我踩过的那些“坑”

代码跑通了,图也画出来了,是不是就万事大吉了?别急,在实际项目中,有几个坑我几乎每次都会提醒自己和团队注意,这也是新手最容易出错的地方。

第一个大坑:数据源的时区与时间格式。 很多公开数据接口返回的时间是UTC时间或带时区信息,而我们的分析往往需要本地时间。直接用pd.to_datetime()转换有时会出错。我的经验是,在读取数据后,立刻检查时间列的前几行,明确时区。如果是UTC,就用df['time'] = pd.to_datetime(df['time']).dt.tz_convert('Asia/Shanghai')进行转换。确保所有时间都在同一时区下,否则你的“日变化规律”可能会出现诡异的偏移。

第二个坑:浓度单位的统一。 国标中,CO的单位是毫克/立方米(mg/m³),而其他污染物是微克/立方米(μg/m³)。有些数据集可能把CO也记录成了μg/m³,如果不加转换直接计算,结果会错得离谱。所以,在计算前,一定要确认你的数据框中CO浓度的单位。一个简单的检查:print(df['CO'].mean()),如果数值在0.5-5之间,单位很可能是mg/m³;如果是几百到几千,那可能就是μg/m³,需要除以1000进行转换。

第三个细节:四舍五入的规则。 国标规定IAQI和AQI都取整数。Python内置的round()函数在处理.5时,采用的是“银行家舍入法”(round half to even),这有时会导致与预期不符。比如2.5会被舍入为2,3.5会被舍入为4。为了严格符合常规的“四舍五入”,我通常使用int(round(x + 1e-9))这个小技巧,或者使用Decimal模块进行更精确的控制。在数学建模等严谨场合,这个细节必须注意。

第四个进阶点:缺失值处理。 真实数据常有缺失。如果某个时间点缺少一种污染物的数据,怎么算AQI?简单的做法是跳过该污染物,用剩余的有效污染物计算最大值。但更严谨的做法是,根据数据缺失的严重程度,决定是插值(如用前后时刻的平均值)、删除该时间点,还是标记为“数据不足”。这需要结合业务需求来判断。

第五个思考:1小时AQI vs 24小时AQI的应用场景。 我们代码里实现了两种。那什么时候用哪种呢?1小时AQI更像“实时快照”,适用于健康出行提示(比如现在出门跑步合不合适)、污染事件的实时追踪。24小时AQI则反映长期暴露水平,更适合用于每日空气质量评价、健康影响评估、以及向公众发布的“今日空气质量”。在你的分析报告中,一定要明确说明你使用的是哪种指数,为什么选用它。

最后,我想说,亲手实现一遍AQI计算,最大的收获不是学会了几个函数,而是建立了一种“数据透视”能力。下次当你再看到App上那个简单的AQI数字时,你脑子里能立刻浮现出背后六条浓度曲线如何交织、竞争,最终由那个“短板”决定全局的图景。这种从底层理解数据生成逻辑的能力,会让你在面临更复杂的环境数据分析任务,比如预测建模、溯源分析时,拥有更扎实的起点和更清晰的思路。希望这份详细的实战指南能真正帮到你,如果在复现代码时遇到任何问题,随时可以再来交流。

更多推荐