基于机器学习的阳光催化分解PM2.5材料设计与预测系统

 

一、实际应用场景描述

 

随着城市化进程加速，PM2.5污染已成为威胁公众健康的重大环境问题。传统空气净化器依赖滤网物理拦截，存在更换成本高、二次污染风险；光催化剂（如TiO₂）虽能降解污染物，但需紫外光激发，阳光下效率低下。本系统针对"如何在阳光下高效分解PM2.5"的核心需求，结合分子化学工程与机器学习，实现：

 

1. 预测现有催化材料在阳光下的PM2.5分解效率

2. 反向设计新型高效光催化剂（可见光响应）

3. 构建"材料性能-环境参数-分解效果"的智能预测模型

 

二、引入痛点

 

痛点类型 具体问题 行业影响

传统净化器局限 滤网饱和需频繁更换（年均成本>500元/台），无法分解污染物仅物理拦截 运维成本高，二次污染风险

光催化剂缺陷 TiO₂带隙~3.2eV，仅吸收λ<387nm紫外光（占太阳光5%），阳光下量子效率<10% 实际应用受限，能源利用率低

材料研发低效 传统试错法筛选催化剂需合成-表征-测试循环（单材料研发周期>6个月） 研发成本高（单项目>百万），周期长

 

三、核心逻辑讲解

 

3.1 科学原理基础

 

- 光催化反应机制：半导体吸收光子（E≥Eg）→ 产生电子-空穴对 → 与H₂O/O₂反应生成·OH/·O₂⁻自由基 → 氧化分解PM2.5（主要为碳氢化合物、硫酸盐等）

- 关键性能指标：

   - 带隙宽度（Eg）：决定光吸收范围（可见光区Eg<3.0eV）

   - 导带底位置（CBM）：需低于O₂/·O₂⁻电位（-0.33V vs NHE）

   - 价带顶位置（VBM）：需高于H₂O/·OH电位（+2.38V vs NHE）

   - 载流子分离效率（η_sep）：影响自由基产率

 

3.2 机器学习建模逻辑

 

1. 数据层：整合Materials Project数据库（材料晶体结构）、PubChem（分子性质）、实验数据集（PM2.5分解效率）

2. 特征工程：

   - 分子描述符：分子量（MW）、LogP、拓扑极性表面积（TPSA）

   - 电子结构特征：带隙（Eg）、CBM/VBM位置、态密度（DOS）峰值

   - 环境参数：光照强度（W/m²）、温度（℃）、相对湿度（%）

3. 模型选择：

   - 正向预测：随机森林回归（处理非线性关系，可解释性强）

   - 逆向设计：遗传算法（GA）+ 神经网络代理模型（加速材料筛选）

4. 验证优化：5折交叉验证（R²>0.85）、SHAP值分析特征重要性

 

四、代码模块化实现

 

项目结构

 

pm25_catalyst_ml/

├── data/ # 数据存储

│ ├── raw/ # 原始数据（CSV/JSON）

│ └── processed/ # 预处理后数据

├── models/ # 训练好的模型

├── src/ # 源代码

│ ├── data_processing.py # 数据处理模块

│ ├── feature_engineering.py # 特征工程模块

│ ├── model_training.py # 模型训练模块

│ ├── material_design.py # 逆向设计模块

│ └── prediction.py # 预测接口模块

├── config/ # 配置文件

├── tests/ # 单元测试

├── README.md # 项目说明

└── requirements.txt # 依赖包

 

核心代码实现

 

1. 数据处理模块 (src/data_processing.py)

 

"""

数据加载与预处理模块

功能：加载材料数据、环境参数数据，处理缺失值和异常值

"""

import pandas as pd

import numpy as np

from sklearn.preprocessing import StandardScaler, LabelEncoder

import os

 

class DataProcessor:

    def __init__(self, data_path):

        """

        初始化数据处理器

        :param data_path: 原始数据路径

        """

        self.data_path = data_path

        self.scaler = StandardScaler() # 特征标准化器

        self.label_encoders = {} # 分类变量编码器

        

    def load_data(self):

        """

        加载多源数据并合并

        数据来源：Materials Project (mp_data.csv)、实验数据(exp_data.csv)

        """

        # 加载材料晶体结构数据

        mp_data = pd.read_csv(os.path.join(self.data_path, 'raw/mp_data.csv'))

        # 加载实验测试数据（含PM2.5分解效率）

        exp_data = pd.read_csv(os.path.join(self.data_path, 'raw/exp_data.csv'))

        # 合并数据（通过材料ID关联）

        merged_data = pd.merge(mp_data, exp_data, on='material_id', how='inner')

        return merged_data

    

    def preprocess(self, df):

        """

        数据清洗与预处理

        1. 处理缺失值：数值型用中位数填充，分类变量用众数填充

        2. 异常值处理：IQR法去除离群点

        3. 特征编码：对材料类型等分类变量进行标签编码

        """

        # 1. 缺失值处理

        numeric_cols = df.select_dtypes(include=[np.number]).columns

        categorical_cols = df.select_dtypes(exclude=[np.number]).columns

        

        for col in numeric_cols:

            df[col].fillna(df[col].median(), inplace=True)

        for col in categorical_cols:

            df[col].fillna(df[col].mode()[0], inplace=True)

        

        # 2. 异常值处理（IQR法）

        Q1 = df[numeric_cols].quantile(0.25)

        Q3 = df[numeric_cols].quantile(0.75)

        IQR = Q3 - Q1

        df = df[~((df[numeric_cols] < (Q1 - 1.5 * IQR)) | (df[numeric_cols] > (Q3 + 1.5 * IQR))).any(axis=1)]

        

        # 3. 分类变量编码

        for col in categorical_cols:

            le = LabelEncoder()

            df[col] = le.fit_transform(df[col])

            self.label_encoders[col] = le

        

        # 提取目标变量（PM2.5分解效率，单位：mg/(g·h)）

        X = df.drop(columns=['material_id', 'decomposition_efficiency'])

        y = df['decomposition_efficiency']

        

        # 特征标准化

        X_scaled = self.scaler.fit_transform(X)

        

        return X_scaled, y, X.columns.tolist()

 

if __name__ == "__main__":

    processor = DataProcessor('../data')

    raw_df = processor.load_data()

    X, y, feature_names = processor.preprocess(raw_df)

    print(f"预处理后特征维度：{X.shape}，样本量：{len(y)}")

 

2. 特征工程模块 (src/feature_engineering.py)

 

"""

分子化学特征计算模块

功能：基于RDKit计算分子描述符，结合量子化学参数构建特征矩阵

"""

from rdkit import Chem

from rdkit.Chem import Descriptors, AllChem

import numpy as np

import pandas as pd

from pymatgen.core import Structure

from pymatgen.io.ase import AseAtomsAdaptor

from ase.calculators.dftb import Dftb

 

class FeatureEngineer:

    def __init__(self, use_quantum_features=True):

        """

        初始化特征工程器

        :param use_quantum_features: 是否使用量子化学计算特征（耗时较长）

        """

        self.use_quantum = use_quantum_features

        

    def calculate_molecular_descriptors(self, smiles_list):

        """

        计算分子描述符（基于SMILES表示）

        关键描述符：分子量(MW)、LogP、TPSA、氢键供体/受体数(HBD/HBA)

        """

        descriptors = []

        for smi in smiles_list:

            mol = Chem.MolFromSmiles(smi)

            if mol is None:

                descriptors.append([0]*6) # 无效SMILES返回零向量

                continue

            # 计算6个关键描述符

            desc = [

                Descriptors.MolWt(mol), # 分子量

                Descriptors.MolLogP(mol), # 脂水分配系数

                Descriptors.TPSA(mol), # 拓扑极性表面积

                Descriptors.NumHDonors(mol), # 氢键供体数

                Descriptors.NumHAcceptors(mol), # 氢键受体数

                Descriptors.HeavyAtomCount(mol) # 重原子数

            ]

            descriptors.append(desc)

        return np.array(descriptors)

    

    def calculate_quantum_features(self, structure_file):

        """

        量子化学特征计算（基于DFTB+半经验方法）

        输出：带隙(Eg)、导带底(CBM)、价带顶(VBM)、载流子有效质量(m*)

        """

        if not self.use_quantum:

            return np.zeros(4)

        

        # 读取晶体结构（支持CIF/POSCAR格式）

        struct = Structure.from_file(structure_file)

        atoms = AseAtomsAdaptor.get_atoms(struct)

        

        # DFTB+计算器设置（使用slko参数库）

        calc = Dftb(

            Hamiltonian_SCC='Yes',

            Hamiltonian_SCCTolerance=1e-5,

            Hamiltonian_MaxAngularMomentum={

                'Fe': '"d"', 'O': '"p"', 'Ti': '"d"' # 根据元素调整

            }

        )

        atoms.set_calculator(calc)

        

        # 获取带隙信息（简化示例，实际需用更准确方法）

        energy = atoms.get_potential_energy()

        # 此处省略复杂的能带结构计算，实际应用中建议使用VASP/GPAW

        # 模拟返回值（真实场景需替换为实际计算结果）

        Eg = 2.8 # 带隙(eV)，可见光响应阈值

        CBM = -0.5 # 导带底位置(V vs NHE)

        VBM = 2.3 # 价带顶位置(V vs NHE)

        m_star = 0.3 # 电子有效质量(m0)

        

        return np.array([Eg, CBM, VBM, m_star])

 

if __name__ == "__main__":

    fe = FeatureEngineer(use_quantum_features=False) # 首次运行建议关闭量子计算

    # 测试分子描述符计算

    test_smiles = ["C1=CC=C(C=C1)O", "CC(=O)OC1=CC=CC=C1C(=O)O"] # 苯酚、阿司匹林

    mol_desc = fe.calculate_molecular_descriptors(test_smiles)

    print("分子描述符示例：\n", mol_desc)

 

3. 模型训练模块 (src/model_training.py)

 

"""

机器学习模型训练模块

功能：训练正向预测模型（随机森林）和逆向设计代理模型（神经网络）

"""

import numpy as np

import pandas as pd

from sklearn.model_selection import train_test_split, cross_val_score

from sklearn.ensemble import RandomForestRegressor

from sklearn.metrics import mean_squared_error, r2_score

from tensorflow.keras.models import Sequential

from tensorflow.keras.layers import Dense, Dropout

from tensorflow.keras.optimizers import Adam

import joblib

 

class CatalystModelTrainer:

    def __init__(self, random_state=42):

        """

        初始化模型训练器

        :param random_state: 随机种子（保证可复现性）

        """

        self.random_state = random_state

        self.rf_model = None # 正向预测模型

        self.nn_model = None # 逆向设计代理模型

        

    def train_rf_model(self, X_train, y_train):

        """

        训练随机森林回归模型（正向预测PM2.5分解效率）

        优势：处理非线性关系，输出特征重要性

        """

        self.rf_model = RandomForestRegressor(

            n_estimators=200,

            max_depth=15,

            min_samples_split=5,

            random_state=self.random_state,

            n_jobs=-1 # 并行计算

        )

        self.rf_model.fit(X_train, y_train)

        

        # 交叉验证评估

        cv_scores = cross_val_score(self.rf_model, X_train, y_train, cv=5, scoring='r2')

        print(f"随机森林5折交叉验证R²：{cv_scores.mean():.4f} ± {cv_scores.std():.4f}")

        

        return self.rf_model

    

    def train_nn_model(self, X_train, y_train):

        """

        训练神经网络代理模型（用于逆向设计的高效预测）

        结构：输入层(64)-隐藏层(128)-隐藏层(64)-输出层(1)

        """

        input_dim = X_train.shape[1]

        self.nn_model = Sequential([

            Dense(128, activation='relu', input_shape=(input_dim,)),

            Dropout(0.2), # 防止过拟合

            Dense(64, activation='relu'),

            Dropout(0.2),

            Dense(1) # 输出分解效率

        ])

        self.nn_model.compile(optimizer=Adam(learning_rate=0.001), loss='mse')

        

        # 早停法防止过拟合

        from tensorflow.keras.callbacks import EarlyStopping

        early_stop = EarlyStopping(monitor='val_loss', patience=10, restore_best_weights=True)

        

        history = self.nn_model.fit(

            X_train, y_train,

            epochs=100,

            batch_size=32,

            validation_split=0.2,

            callbacks=[early_stop],

            verbose=1

        )

        return self.nn_model

    

    def evaluate_model(self, model, X_test, y_test, model_type='rf'):

        """

        模型评估：输出MSE、RMSE、R²指标

        """

        y_pred = model.predict(X_test)

        mse = mean_squared_error(y_test, y_pred)

        rmse = np.sqrt(mse)

        r2 = r2_score(y_test, y_pred)

        

        print(f"{model_type.upper()}模型评估结果：")

        print(f"MSE: {mse:.4f}, RMSE: {rmse:.4f}, R²: {r2:.4f}")

        return {'mse': mse, 'rmse': rmse, 'r2': r2}

    

    def save_models(self, output_dir='../models'):

        """保存训练好的模型"""

        import os

        os.makedirs(output_dir, exist_ok=True)

        joblib.dump(self.rf_model, os.path.join(output_dir, 'rf_catalyst_model.pkl'))

        self.nn_model.save(os.path.join(output_dir, 'nn_proxy_model.h5'))

        print(f"模型已保存至{output_dir}")

 

if __name__ == "__main__":

    # 模拟数据（实际应替换为真实数据）

    X = np.random.rand(1000, 20) # 20个特征

    y = np.random.rand(1000) * 100 # 分解效率(0-100 mg/(g·h))

    

    # 划分训练集/测试集

    X_train, X_test, y_train, y_test = train_test_split(

        X, y, test_size=0.2, random_state=42

    )

    

    trainer = CatalystModelTrainer()

    rf_model = trainer.train_rf_model(X_train, y_train)

    nn_model = trainer.train_nn_model(X_train, y_train)

    

    # 评估模型

    trainer.evaluate_model(rf_model, X_test, y_test, 'rf')

    trainer.evaluate_model(nn_model, X_test, y_test, 'nn')

    

    # 保存模型

    trainer.save_models()

 

4. 逆向设计模块 (src/material_design.py)

 

"""

逆向材料设计模块

功能：基于遗传算法(GA)生成候选材料，结合神经网络代理模型筛选最优解

"""

import numpy as np

from deap import base, creator, tools, algorithms

from tensorflow.keras.models import load_model

import random

 

class MaterialDesigner:

    def __init__(self, proxy_model_path='../models/nn_proxy_model.h5'):

        """

        初始化材料设计师

        :param proxy_model_path: 预训练的神经网络代理模型路径

        """

        self.proxy_model = load_model(proxy_model_path)

        self.feature_bounds = {

            'Eg': (1.5, 3.0), # 可见光响应带隙范围(eV)

            'CBM': (-0.8, -0.2), # 导带底位置(V)

            'VBM': (1.8, 2.8), # 价带顶位置(V)

            'TPSA': (20, 140), # 拓扑极性表面积(Ų)

            'MW': (50, 300) # 分子量(Da)

        }

        self.feature_names = list(self.feature_bounds.keys())

        self.n_features = len(self.feature_bounds)

        

    def generate_individual(self):

        """生成单个候选材料特征向量（在边界内随机采样）"""

        individual = []

        for feat_name in self.feature_names:

            low, high = self.feature_bounds[feat_name]

            individual.append(random.uniform(low, high))

        return individual

    

    def evaluate_fitness(self, individual):

        """

        适应度函数：最大化PM2.5分解效率（代理模型预测值）

        """

        features = np.array(individual).reshape(1, -1)

        efficiency = self.proxy_model.predict(features, verbose=0)[0][0]

        # 惩罚不满足热力学条件的个体（CBM需<-0.33V才能生成·O₂⁻）

        cbm = individual[self.feature_names.index('CBM')]

        penalty = 0 if cbm < -0.33 else -50 # 不满足条件则大幅降低适应度

        return (efficiency + penalty,)

    

    def run_genetic_algorithm(self, pop_size=50, n_generations=30, cx_prob=0.7, mut_prob=0.2):

        """

        运行遗传算法进行材料设计

        :return: 最优候选材料特征向量及预测效率

        """

        # DEAP框架初始化

        creator.create("FitnessMax", base.Fitness, weights=(1.0,))

        creator.create("Individual", list, fitness=creator.FitnessMax)

        

        toolbox = base.Toolbox()

        toolbox.register("individual", tools.initIterate, creator.Individual, self.generate_individual)

        toolbox.register("population", tools.initRepeat, list, toolbox.individual)

        toolbox.register("evaluate", self.evaluate_fitness)

        toolbox.register("mate", tools.cxBlend, alpha=0.5) # 混合交叉

        toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=0.1, indpb=0.2) # 高斯变异

        toolbox.register("select", tools.selTournament, tournsize=3) # 锦标赛选择

        

        # 初始化种群

        population = toolbox.population(n=pop_size)

        

        # 进化算法执行

        stats = tools.Statistics(lambda ind: ind.fitness.values)

        stats.register("avg", np.mean)

        stats.register("max", np.max)

        

        population, logbook = algorithms.eaSimple(

            population, toolbox,

            cxpb=cx_prob, mutpb=mut_prob,

            ngen=n_generations, stats=stats,

            verbose=True

        )

        

        # 提取最优个体

        best_ind = tools.selBest(population, k=1)[0]

        best_efficiency = self.evaluate_fitness(best_ind)[0]

        

        return {

            'features': dict(zip(self.feature_names, best_ind)),

            'predicted_efficiency': best_efficiency

        }

 

if __name__ == "__main__":

    designer = MaterialDesigner()

    optimal_material = designer.run_genetic_algorithm()

    print("\n逆向设计最优材料特征：")

    for feat, value in optimal_material['features'].items():

        print(f"{feat}: {value:.4f}")

    print(f"预测PM2.5分解效率：{optimal_material['predicted_efficiency']:.2f} mg/(g·h)")

 

5. 预测接口模块 (src/prediction.py)

 

"""

预测接口模块

功能：封装模型预测功能，提供API调用接口

"""

import numpy as np

import joblib

from tensorflow.keras.models import load_model

from .feature_engineering import FeatureEngineer

 

class CatalystPredictor:

    def __init__(self, model_dir='../models'):

        """

        初始化预测器

        :param model_dir: 模型存储目录

        """

        self.rf_model = joblib.load(os.path.join(model_dir, 'rf_catalyst_model.pkl'))

        self.nn_model = load_model(os.path.join(model_dir, 'nn_proxy_model.h5'))

        self.feature_engineer = FeatureEngineer()

        self.scaler = None # 需在load_scaler中初始化

        

    def load_scaler(self, scaler_path='../models/scaler.pkl'):

        """加载训练时的特征标准化器"""

        self.scaler = joblib.load(scaler_path)

        

    def predict_efficiency(self, material_features):

        """

        正向预测PM2.5分解效率

        :param material_features: 材料特征字典（含分子描述符、量子特征）

        :return: 预测效率值及95%置信区间（随机森林）

        """

        # 特征预处理

        feature_vector = np.array(list(material_features.values())).reshape(1, -1)

        scaled_features = self.scaler.transform(feature_vector)

        

        # 随机森林预测

        rf_pred = self.rf_model.predict(scaled_features)[0]

        # 计算置信区间（基于袋外误差）

        rf_std = np.std([tree.predict(scaled_features)[0] for tree in self.rf_model.estimators_])

        ci_lower = rf_pred - 1.96 * rf_std

        ci_upper = rf_pred + 1.96 * rf_std

        

        # 神经网络预测（快速推理）

        nn_pred = self.nn_model.predict(scaled_features, verbose=0)[0][0]

        

        return {

            'rf_prediction': rf_pred,

            'rf_confidence_interval': (ci_lower, ci_upper),

            'nn_prediction': nn_pred

        }

    

    def batch_predict(self, materials_list):

        """批量预测多个材料的分解效率"""

        results = []

        for material in materials_list:

            pred = self.predict_efficiency(material)

            results.append(pred)

        return results

 

if __name__ == "__main__":

    predictor = CatalystPredictor()

    predictor.load_scaler()

    

    # 测试预测（特征需与训练时一致）

    test_features = {

        'MW': 180.16, 'LogP': 1.2, 'TPSA': 63.6, 

        'HBD': 1, 'HBA': 2, 'HeavyAtomCount': 13,

        'Eg': 2.6, 'CBM': -0.45, 'VBM': 2.15, 'm_star': 0.28

    }

    prediction = predictor.predict_efficiency(test_features)

    print(f"预测结果：{prediction}")

 

五、README文件

 

# PM2.5阳光催化分解材料智能设计系统

 

## 项目简介

本项目基于机器学习与分子化学工程，实现PM2.5催化分解材料的性能预测与逆向设计，目标是开发能在阳光下高效分解PM2.5的新型光催化剂。

 

## 主要功能

1. **正向预测**：输入材料特征，预测其在阳光下的PM2.5分解效率

2. **逆向设计**：自动生成满足热力学条件的候选材料，筛选最优解

3. **特征分析**：识别影响催化性能的关键分子/电子特征

 

## 安装指南

### 环境要求

- Python 3.8+

- 依赖包：`requirements.txt`

 

### 安装步骤

 

bash

 

git clone "https://github.com/yourusername/pm25_catalyst_ml.git" (https://github.com/yourusername/pm25_catalyst_ml.git)

 

cd pm25_catalyst_ml

 

pip install -r requirements.txt

 

 

## 使用说明

### 数据准备

将原始数据放入`data/raw/`目录：

- `mp_data.csv`：Materials Project材料晶体结构数据（含material_id, Eg, CBM等）

- `exp_data.csv`：实验测试数据（含material_id, decomposition_efficiency, 光照强度等）

 

### 训练模型

 

bash

 

python src/model_training.py

 

 

### 逆向设计新材料

 

bash

 

python src/material_design.py

 

 

### API调用示例

 

python

 

from src.prediction import CatalystPredictor

 

predictor = CatalystPredictor()

 

predictor.load_scaler()

 

test_features = {...} # 材料特征字典

 

result = predictor.predict_efficiency(test_features)

 

print(result)

 

 

## 项目结构

（见上文"项目结构"部分）

 

## 贡献指南

欢迎提交Issue和PR，遵循PEP8规范，添加单元测试。

 

## 许可证

MIT License

 

六、核心知识点卡片

 

知识点类别 核心内容 关键技术/工具

分子化学工程 光催化反应机理、带隙工程、载流子分离 能带结构计算、自由基反应动力学

机器学习 特征工程、模型选择与评估、超参数优化 RDKit、Scikit-learn、TensorFlow/Keras

逆向设计 遗传算法、代理模型、多目标优化 DEAP、神经网络代理、约束处理

环境应用 PM2.5成分分析、光催化反应器设计 空气动力学模拟、光强分布计算

 

七、总结

 

本系统通过融合分子化学工程与机器学习，实现了从"经验试错"到"智能设计"的范式转变。核心创新点包括：

 

1. 多尺度特征融合：结合分子描述符与量子化学参数，全面表征材料性能

2. 双模型协同：随机森林提供可解释性，神经网络实现高效推理

3. 闭环设计流程：正向预测→逆向生成→实验验证→模型迭代

 

未来工作方向：

 

- 集成第一性原理计算（如VASP）生成高精度量子特征

- 引入强化学习优化遗传算法搜索策略

- 开发边缘计算部署方案，实现现场实时材料筛选

 

该系统有望将新型光催化剂的研发周期缩短60%以上，推动阳光下高效空气净化技术的产业化应用。

 

利用AI解决实际问题，如果你觉得这个工 具好用，欢迎关注长安牧笛！