基于动态约束插值的PM10浓度空间分布模拟

Spatial Distribution Simulation of PM10 Concentrations Based on Dynamic Constraint Interpolation

摘要: PM10是一类极具代表性的大气污染物，其表面易吸附多环芳烃、重金属等有毒有害物质，进而对人体健康与生态环境构成多重威胁，相关研究因此受到广泛重视。对PM10污染的研究不仅能深入剖析其传输扩散规律，更为后续污染防控工作提供坚实的科学支撑。在PM10相关研究中，精准模拟其浓度空间分布是核心环节，对全面掌握PM10污染特征具有不可替代的作用。地面监测、数值模拟与空间插值是当前研究PM10浓度空间分布的主要技术手段，但数值模型对数据完整性要求较高，传统插值方法在观测稀疏区域难以达到理想精度。动态约束插值法(DCIM)已被证实可有效提升时空稀疏观测数据的利用率，为此，本研究提出一种基于伴随模型的动态约束插值方法，通过算法优化增强PM10浓度的空间重建能力。利用2015年3月2日至3月5日中国中南部85个地面监测站的观测数据，通过分步实验设计开展系统验证：首先验证DCIM在PM10数值模拟中的有效性，实验结果显示该方法显著提升了模拟精度、优化了模拟效果，充分证明其在该领域的适用性与可靠性；然后针对插值方法在稀疏数据场中的局限性，引入高阶守恒插值(PPM)算法，实验表明该算法在适配稀疏数据场景方面具有显著优势。

Abstract: PM10 is a highly representative atmospheric pollutant whose surface readily adsorbs toxic substances such as polycyclic aromatic hydrocarbons and heavy metals, posing multiple threats to human health and the ecological environment. Consequently, research in this area has garnered significant attention. Studies on PM10 pollution not only provide in-depth insights into its transport and dispersion patterns but also offer robust scientific support for subsequent pollution prevention and control efforts. Accurate simulation of PM10’s spatial concentration distribution is central to such research, playing an irreplaceable role in comprehensively understanding its pollution characteristics. Ground-based monitoring, numerical modeling, and spatial interpolation are the primary techniques currently used to study PM10’s spatial distribution. However, numerical models demand high data completeness, while traditional interpolation methods struggle to achieve ideal accuracy in sparsely observed regions. Dynamic Constrained Interpolation Method (DCIM) has been proven effective in enhancing the utilization of spatiotemporally sparse observational data. Therefore, this study proposes a dynamic constrained interpolation method based on an adjoint model, which enhances the spatial reconstruction capability of PM10 concentrations through algorithm optimization. Using observational data from 85 ground monitoring stations in central and southern China between March 2 and March 5, 2015, systematic validation was conducted through a stepwise experimental design: First, DCIM’s effectiveness in PM10 numerical simulation was validated. Results demonstrated significant improvements in simulation accuracy and optimized outcomes, fully confirming its applicability and reliability in this domain. Subsequently, addressing limitations of traditional interpolation methods in sparse data fields, the higher-order conservation interpolation (PPM) algorithm was introduced. Experiments revealed that this algorithm possesses distinct advantages in adapting to sparse data scenarios.

1. 引言

大气颗粒物PM10因其复杂的物理化学特性及对公众健康的显著影响，长期以来一直是全球环境与气候研究的重点方向。PM10可通过呼吸道深入肺泡，诱发哮喘、慢性阻塞性肺病等呼吸系统疾病[1]，其表面吸附的多环芳烃(PAHs)和重金属等有毒物质，还可通过氧化应激机制加剧心血管损伤与DNA突变风险[2]。流行病学调研数据显示，PM10浓度升高时呼吸系统疾病相关死亡率也升高[3]。在生态系统层面，PM10的长期沉降会改变土壤pH值与微生物群落结构，同时遮蔽阳光抑制植物光合作用，加速区域生态系统的脆弱化进程。

在区域环境与社会经济层面，中国中南部作为国家核心经济增长区，其主要工业城市与人口密集区域面临严重的PM10污染问题，长期浓度超出世界卫生组织标准2~3倍甚至更高[4]。鉴于PM10对人体健康与生态环境的多重危害，开展深入研究极具必要性。自2013年起，中国环境保护部构建了全国空气质量监测网络，该网络每小时发布PM10等空气污染物浓度数据，并启动专项科研项目鼓励相关研究探索。通过分析不同区域PM10浓度分布特征，可为空气污染状况的模拟预测及防控政策制定提供科学依据。

当前，地面监测站能够提供高精度、高时间分辨率的观测数据，已广泛应用于城市空气质量实时警报与健康暴露评估，例如中国国家环境监测中心(CNEMC)已通过全国1400余个站点构建覆盖339个地级市的监测网络，为《大气污染防治行动计划》的实施提供了关键数据支撑[5]。数值模型在PM10浓度空间分布研究中发挥着关键作用，但这类模型对数据完整性的高度依赖限制了其应用场景[6]。传统空间插值技术如Cressman插值、Kriging插值等虽在站点数据空间化中应用广泛，但在观测数据稀缺区域易出现过度平滑的问题，难以精准反映污染物浓度的真实空间分布特征。

动态约束插值(DCIM)通过观测数据与数值模型，为上述问题提供了创新解决方案，显著提升了稀疏数据的利用率。在海洋学领域，Zheng等人利用DCIM方法重建渤海湾潮汐场，相比传统数学插值方法，均方根误差(RMSE)实现显著降低[7]；在大气科学领域，Li等人通过该方法将全国PM2.5监测数据同化的空间分辨率提升至5千米，展现出其在大规模污染物分布重建中的潜力[8]。空间插值方法的准确性与精度一直是研究关注的重点，Cressman插值虽具备简单高效的优势，但在处理远离观测站的格点时易产生较大误差，影响半径的不合理设定会进一步削弱插值效果。相比之下，高阶守恒插值(PPM)能够保证局部质量守恒，具备更高精度，适用于大气学、流体力学等领域。Li等人(2017)应用PPM有效弥补了Cressman插值在大网格间距下的精度缺陷，其能够通过少量独立点构建更贴合实际的初始场，即便在大网格间距条件下仍能保持较高精度，为优化动态约束插值方法的模拟效果提供了新思路，对PM10浓度空间重建研究具有重要参考价值[9]。因此，本研究旨在结合地面监测数据与高阶守恒插值(PPM)优化算法，提升中国中南部PM10浓度空间重建的精度。

本文后续内容安排如下：第2节详细阐述所用材料与方法；第3节呈现实际实验与分析结果；第4节总结讨论并提出结论。

2. 材料与方法

2.1. 研究区域与数据

近年来，中国持续推进基础监测网络的科学规划与部署，截至2022年，已建成超过1400个地级监测站，覆盖339个地级市，形成了较为广泛的监测网络。本研究聚焦中国中南部地区，该区域人口密集、经济发展迅速，在全国城市体系中占据关键地位。研究选取该地区85个地面监测站，通过网络爬虫技术获取2015年主要城市的PM10观测数据，数据覆盖2015年3月2日0点至3月5日0点的72小时观测期。将每小时浓度值取平均得到日平均浓度，对于拥有多个监测站的城市，取所有可用测量值的平均值85个监测站的空间分布如图1所示。

Figure 1. Spatial distribution of 85 basic monitoring stations in the central and southern regions

图1. 中南部地区85个基础监测站的空间分布

2.2. 动态约束插值方法

伴随模型通过求解伴随方程计算目标函数梯度，结合梯度下降法等优化算法逐步逼近最优解，该方法广泛应用于海洋与大气研究中的参数优化及动态约束插值。例如，Wu等人(2021)利用该方法优化Ekman 模型参数[10]，Xu等人(2017)通过其高效模拟2014年APEC期间中国区域PM2.5浓度分布[11]，Li等人(2021)成功模拟全国PM2.5浓度分布[8]。这些研究证实，该方法能有效提升时空稀疏数据的利用率与模拟精度，尤其适用于大气污染物研究，在全国PM10观测数据重建中具有显著优势。

2.2.1. 正向模型

PM10浓度是空气质量评估的核心指标，其空间分布特征与变化趋势的精准刻画可为空气质量管理提供重要依据。大气中的PM10受到传输、扩散、光化学转化(产生二次PM10)及湿干沉降等多重耦合过程影响，其传输主要由风向与湍流扩散驱动。本研究采用二维对流扩散方程描述其动力学机制，将一次与二次PM10合并为“源汇项”(S)，同时忽略垂直扩散与复杂化学过程，建立简化质量守恒方程：

$\frac{\partial C}{\partial t}+u\frac{\partial C}{\partial x}+v\frac{\partial C}{\partial y}=\frac{\partial }{\partial x}\left(A_H\frac{\partial C}{\partial x}\right)+\frac{\partial }{\partial y}\left(A_H\frac{\partial C}{\partial y}\right)+S$

其中，C表示PM10浓度，t表示时间，x和y为水平坐标，u和v分别为x和y坐标处的水平风速，$A_H$ 为水平粘度系数，S表示源与汇。

$C\left(t^0,x,y\right)=C^0\left(x,y\right)$

$\frac{\partial C}{\partial t}=0\text{on}{\Gamma }^{\text{IN}}$

$\frac{\partial C}{\partial \stackrel{\rightharpoonup }{n}}=0\text{on}{\Gamma }^{\text{OUT}}$

在流入边界处施加常数边界条件，而在流出边界处施加零梯度边界条件。该模型的数值方案与文献[12]中描述的方案一致。

2.2.2. 伴随模型

伴随模型是解决参数优化的强大工具。通过计算目标函数对参数的梯度，并结合优化算法(如梯度下降法)，它们能逐步逼近最优解。在大气科学领域，伴随模型常用于参数估计。该过程通过优化正向模型参数，使模拟结果与观测数据的偏差最小化，从而获得更精确的时空分布。

基于此，我们构建了伴随模型并定义了成本函数，以量化观测与模拟之间的差异：

$J=\frac{1}{2}{\displaystyle {\int }_{\sum }K{\left(C-C_{OBS}\right)}^2\text{d}\sigma }$

拉格朗日函数：

$L=J+\left\{p\times \left[\frac{\partial C}{\partial t}+u\frac{\partial C}{\partial t}+v\frac{\partial C}{\partial y}-\frac{\partial }{\partial x}\left(A_H\frac{\partial C}{\partial x}\right)-\frac{\partial }{\partial y}\left(A_H\frac{\partial C}{\partial y}\right)-S\right]\right\}\text{d}\sigma$

其中，C和$C_{OBS}$ 分别表示模拟与观测的PM10浓度；p是C的伴随变量；K是权重矩阵，理论上为观测误差协方差矩阵的逆矩阵。假设误差不相关且权重相等，则可简化为1 (含观测值)或0 (无观测值)。为了最小化成本函数，如下所示：

$\frac{\partial L}{\partial p}=0$

$\frac{\partial L}{\partial C}=0$

$\frac{\partial L}{\partial A_H}=0$

$\frac{\partial L}{\partial S}=0$

2.2.3. 独立点方案

在空间上选取部分网格点作为“独立点”。其他网格点的PM10浓度通过独立点浓度值进行线性插值计算。该方法既减少了需调整的网格点数量，又提升了反演结果的平滑度。设 $\left\{C_{ii,jj},ii=1,2,\cdots ,in;jj=1,2,\cdots ,jn\right\}$ 为独立点的PM10浓度，$in\times jn$ 为独立点总数，则网格点$\left(i,j\right)$ 的PM10浓度值为：

$C_{i,j}={\displaystyle \sum }_{ii,jj}\left({\varphi }_{i,j,ii,jj}\times C_{ii,jj}\right)$

${\varphi }_{i,j,ii,jj}$ 是线性插值系数：

${\varphi }_{i,j,ii,jj}=W_{i,j,ii,jj}/{\displaystyle \sum _{ii,jj}{W}_{i,j,ii,jj}}$

$W_{i,j,ii,jj}$ 是Cressman插值法的权重系数：

$W_{i,j,ii,jj}=\left(R^2-r_{i,j,ii,jj}^2\right)/\left(R^2+r_{i,j,ii,jj}^2\right)$

$r_{i,j,ii,jj}$ 是一个独立点。该点$\left(ii,jj\right)$ 到网格点$\left(i,j\right)$ 的距离为d，影响半径为R。独立点的梯度公式可直接推导。在初始场反演中，我们采用独立点方案。空间上，每四个网格点选取一个独立点，影响半径设定为4˚。独立点的初始PM10浓度设定为50 μg/m³。

2.2.4. 插值方法

1. Cressman插值

模型运行过程中插值方法的选择对模拟结果影响显著，直接关系到空间数据处理能力和模拟结果的准确性。不同插值方法在原理、算法及适用性方面存在差异。例如：Cressman插值法基于距离与浓度变化采用加权分布法估算未知点浓度；Kriging插值法运用空间统计原理，考虑空间自相关性与变异性以生成更精确连续的浓度分布；正交多项式拟合通过构建多项式函数拟合数据点，适用于分布均匀且趋势明显的数据，可生成平滑的浓度分布曲线。

这些插值方法可归纳为以下公式：

$\widehat{Z}\left(x_0\right)={\displaystyle \sum _{i=1}^n{\lambda }_{i}z\left(x_i\right)}$

Cressman插值法是一种距离加权插值方法。因其模拟精度高、操作简便且适用于网格化数据，常被应用于气象研究。将Cressman插值法应用于观测数据处理，其定义为：

$C\left(x,y\right)=\frac{{\displaystyle {\sum }_{i=1}^n{W}_i{z}_i}}{{\displaystyle {\sum }_{i=1}^n{W}_i}}$

$W_i$ 是第$i$ 个邻近观测点的权重，通常根据$x_0$ 与$x_i$ 的距离使用以下公式计算：

$W_i=\frac{R-d_i}{R}$

此处，$R$ 表示影响半径，$d_i$ 表示从$x_0$ 到$x_i$ 的距离。

2. 高阶守恒插值法(PPM)

PPM方法应用于动态约束插值，二维计算域划分为规则欧拉单元，单元边界与中心定义如下：

${\Omega }_{i,j}=\left[x_{i-1/2},x_{i+1/2}\right]\times \left[y_{i-1/2},y_{i+1/2}\right]$

$x_i=\frac{1}{2}\left(x_{i-1/2}+x_{i+1/2}\right)$

$y_i=\frac{1}{2}\left(y_{j-1/2}+y_{j+1/2}\right)$

设$\left(x_i,y_i\right)$ 表示一个独立点，$C_{i,j}$ 表示该独立点处的参数值。对于二维空间，计算采用分步策略：先沿x方向进行，再沿y方向进行。本节中，鉴于x方向计算原理与y方向类似，我们将重点阐述x方向保守插值过程以突出关键点并避免冗余。通过这一系列计算步骤，最终获得$t=t_n$ 时间处${\Omega }_{i,j}$ 的值。

本应用中PPM ${\left[ℛC^n\right]}_i\left(x\right)$ 的具体应用如下：

${\left[ℛC^n\right]}_i\left(x\right)=C_{i-1/2}^n+\frac{x-x_{i-1/2}}{h}\left(\Delta C_i^n+C_{6,i}^n\frac{x_{i+1/2}-x}{h}\right)$

$\Delta C_i^n={\langle C^n\rangle }_{i+1/2}-{\langle C^n\rangle }_{i-1/2}$

$C_{6,i}^n=6\left(C_i^n-\frac{{\langle C^n\rangle }_{i+1/2}-{\langle C^n\rangle }_{i-1/2}}{2}\right)$

插值分布${\left[ℛC^n\right]}_i\left(x\right)$ 定义为，其中${\Omega }_i=\left(x_{i-1/2},x_{i+1/2}\right)$ ，${\langle C^n\rangle }_{i+1/2}$ 表示$x_{i+1/2}$ 点处的数值。

为获得PPM的二阶近似表达式，给出${\langle C^n\rangle }_{i+1/2}$ 如下公式：

${\langle C^n\rangle }_{i+1/2}=\frac{C_i^n+C_{i+1}^n}{2}$

然后，我们可以使用以下方程获得点上所有值：

$\begin{array}{c}{\left[ℛC^n\right]}_i\left(x\right)=\left[\frac{1}{2}-\frac{x-x_{i-1/2}}{2h}-\frac{3\left(x-x_{i-1/2}\right)\left(x_{i+1/2}-x\right)}{2h^2}\right]C_{i-1}^n\\ +\left[\frac{1}{2}+\frac{3\left(x-x_{i-1/2}\right)\left(x_{i+1/2}-x\right)}{h^2}\right]C_i^n\\ +\left[\frac{x-x_{i-1/2}}{2h}-\frac{3\left(x-x_{i-1/2}\right)\left(x_{i+1/2}-x\right)}{2h^2}\right]C_{i+1}^n\end{array}$

2.3. 模型设置

计算域覆盖东经105˚至125˚、北纬20˚至40˚，水平分辨率是0.5˚ × 0.5˚，网格为x方向43个单元、y方向41个单元。时间步长是3600秒。模拟周期自2015年3月2日00:00起至2015年3月5日00:00止，共计72小时。背景风数据覆盖2015年3月2日至5日，空间分辨率为0.5˚ × 0.5˚，时间分辨率是1小时。迭代步数是100次。模型在流入边界采用常数边界条件，在流出边界采用无梯度边界条件。依据2013年环境保护部颁布的中国新环境空气质量标准，入流边界PM10浓度限值设定为35.0 μg/m³。水平扩散系数取值100.0 m²/s，背景值设定为50 μg/m³。为提高模拟精度，每四个网格选取一个独立网格进行优化，其余网格的PM10浓度通过Cressman法计算。初始场通过对初始时刻观测PM10数据进行Cressman插值获得，以确保动态模型模拟结果的准确性。

3. 实际实验与分析

3.1. 实验设计

为开展动态约束插值性能的计算研究，我们实施了两组实际实验。第一组实验验证动态约束插值法的有效性，为进一步探索该方法与不同插值技术结合应用于PM10数值模拟研究时的效果，我们设计了第二组实验，旨在验证优化算法的动态约束插值法在PM10数值模拟中的优势。所有实验均采用交叉验证设计以确保结果可靠性，因交叉验证是评估插值结果的有效方法[13] [14]。

本实验采用四折交叉验证方案评估动态约束插值法的有效性。具体而言，随机抽取四分之一观测数据排除在同化处理之外，该部分数据作为独立观测数据用于验证同化结果，称为“检验数据”。其余四分之三观测数据作为模型输入数据，称为“同化数据”。

采用平均绝对误差(MAE)和均方根误差(RMSE)作为误差指标。MAE能直观反映模拟值与观测值的偏差幅度，是评估模拟精度的常用指标，但仅考虑误差绝对值而忽略正负方向，可能在正负误差相互抵消时无法准确反映真实误差；RMSE通过对误差进行平方处理，全面考量正负误差，且对较大偏差更敏感，能更综合地评估模拟精度。定义如下：

$\text{MAE}={\displaystyle \sum }_{i=1}^N\left|M_i-O_i\right|/N$

$\text{RMSE}={\left[\frac{1}{N}{\displaystyle \sum }_{i=1}^N{\left(M_i-O_i\right)}^2\right]}^{\frac{1}{2}}$

此处，N表示检验数据点的总数$\left(i=1,2,\cdots ,N\right)$ ，而$M_i$ 和$O_i$ 分别代表PM10浓度的模拟值与观测值。

3.2. 第一组结果与讨论

本研究采用基于Cressman插值的动态约束插值方法，模拟2015年3月2日0:00至3月5日0:00的PM10浓度数据，通过比较模型应用前后成本函数和误差的变化(实验PE11)评估方法性能；设计实验PE12和PE13，PE12计算监测站日均浓度后采用Cressman方法直接插值，得到日均PM10浓度空间分布；PE13以地面监测站实际观测数据为模型输入，通过伴随模型优化初始场和源汇项等关键参数运行模型得到PM10浓度空间分布。通过对比两种方法的浓度空间分布与85个地面监测站的日均浓度分布，评估动态约束插值方法的模拟准确性。

2015年3月2日至3月5日的四折交叉验证实验结果如表1所示，四组实验的约化成本函数下降趋势如图2所示。表1显示，同化后验证点观测值与模拟值的MAE从44.53 μg/m³降至34.34 μg/m³，降幅达22.89%；RMSE显著降低30.41%，表明动态约束插值方法有效减小了模拟结果与实际观测值的误差，使模拟值更接近真实PM10浓度。图2显示成本函数随迭代过程持续下降，进一步证实该方法在提高模拟精度方面的有效性。

图3展示了85个地面监测站的3天平均观测结果，图4展示了传统Cressman插值法和基于Cressman的动态约束插值法(DCIM)的PM10浓度分布图。对比可见，两种方法的平均浓度分布图与平均观测图高度相似，均能有效模拟PM10浓度的空间分布特征。其中，传统Cressman插值法作为经典方法已得到广泛应用[15]，而DCIM方法产生的分布特征更平滑、空间连续性更优，这得益于2015年3月中国中南部处于春季过渡阶段，受东亚冬季风残余影响，区域盛行偏北风，同时气温回升导致边界层高度抬升，湍流扩散能力增强，这些气象条件直接影响PM10的传输与扩散过程。从物理机制来看，动态约束插值法(DCIM)通过耦合二维对流扩散方程与伴随模型，将水平风速(u、v)、水平扩散系数($A_H$ )等参数纳入优化过程，能够动态响应区域气象场的空间差异，例如在风速较大的边界区域，DCIM通过梯度优化捕捉到PM10的远距离传输特征，而在湍流扩散较弱的盆地或河谷区域，则通过源汇项调整还原浓度累积效应，从而更精准地捕捉PM10在复杂气象条件下的传输扩散规律。此外，DCIM结果与2015年春季中国环境空气中PM10浓度的空间分布特征具有一致性，进一步验证了DCIM在空间分布模拟中的有效性。

需要注意的是，本研究的监测站点在空间分布上存在明显稀疏性，尤其是陕西、山西和内蒙古边境地区及内蒙古东部，85个监测站在该区域的平均间距超过50公里，远大于中南部平原城市群的站点密度。这种稀疏性导致观测数据无法充分反映区域内PM10的微尺度空间变异，进而影响插值方法的优化效果：传统Cressman插值法模拟的PM10平均浓度约为106.017 μg/m³，虽低于观测浓度(约118.607 μg/m³)但数值更接近观测值；而DCIM方法模拟的平均浓度约为146.982 μg/m³，显著偏高，推测该区域DCIM方法存在过拟合现象，其核心原因是稀疏站点无法提供足够的约束信息，导致模型对局部观测值的敏感度异常升高，进而放大了局部浓度特征。

Table 1. MAE and RMSE of observed versus simulated values at assimilation and verification points (unit: μg/m³)

表1. 同化点与检验点观测值与模拟值的MAE和RMSE (单位：μg/m³)

	同化点		检验点
	Before	After	Before	After
MAE	44.74	24.55	44.52	34.34
RMSE	619.99	39.47	80.81	56.23

Figure 2. Changes in the cost function across iterations for the four experimental groups in the PE11 quadruple cross experiment

图2. PE11四折交叉实验中四个实验组迭代过程中代价函数的变化

Figure 3. Daily average PM10 concentrations at 85 monitoring stations in central and southern China (00:00 on March 2, 2015 to 00:00 on March 5, 2015). The intensity of each point’s color indicates its daily average PM10 concentration level

图3. 中国中部和南部85个监测站的日均PM10浓度(2015年3月2日00:00至2015年3月5日00:00)。每个点的颜色深浅表示其日均PM10浓度的大小

Figure 4. Average fields obtained from experiments PE12 and PE13 ((a): PE12, (b): PE13)

图4. 实验PE12和PE13所得的平均场((a): PE12, (b): PE13)

3.3. 第二组结果与讨论

本研究采用基于Cressman插值和高阶守恒插值(PPM)的两种动态约束插值方法，模拟分析2015年3月2日0:00至3月5日0:00的PM10浓度数据，通过四折交叉验证实验评估结果并重建PM10浓度空间分布。

表2展示了两种方法在同化前后，同化点和检验点观测值与模拟值之间的MAE和RMSE变化，图5展示了实验PE21和PE22的误差折线图。图表显示，与基于Cressman的插值方法相比，基于PPM的动态约束插值方法在同化数据中产生的误差更小，表明其对同化数据实现了更优的优化。表2显示，基于PPM插值的动态约束插值方法在同化点的性能提升显著：MAE降低了59.6%，RMSE降低了95.36%；而基于Cressman插值的方法同化点MAE降低了45.1%，RMSE降低了93.63%。尽管PPM方法在检验点误差略高于Cressman方法，但其同化点数据误差更小，综合来看在兼顾同化点与检验点模拟精度方面更具优势。

Table 2. Comparison of MAE and RMSE between the two methods (unit: μg/m³)

表2. 两种方法的MAE与RMSE对比(单位：μg/m³)

	同化点MAE	检验点MAE	同化点RMSE	检验点RMSE
	Before	After	Before	After
Cressman	44.74	24.55	44.53	34.34
RMSE	44.74	17.66	44.53	41.43

Figure 5. Mean Absolute Error (MAE) and Root Mean Square Error (RMSE) of data assimilated using corresponding methods in the four-fold cross-validation experiment (blue dashed line: dynamic constrained interpolation based on Cressman interpolation; red solid line: dynamic constrained interpolation based on PPM interpolation)

图5. 四折交叉实验中，采用相应方法同化的数据的平均绝对误差(MAE)和均方根误差(RMSE) (蓝色虚线：基于Cressman插值的动态约束插值；红色实线：基于PPM插值的动态约束插值)

图6展示了实验PE21和PE22所得平均浓度的空间分布：PE21为基于Cressman插值的结果，PE22为基于PPM插值的结果。对比可见，两种方法的平均浓度分布与观测分布高度相似，且均表现出平滑的分布模式，整体一致性良好。值得注意的是，PE22在陕西、山西和内蒙古边境地区及内蒙古东部的过拟合现象得到明显缓解，其模拟的PM10平均浓度约为134.196 μg/m³，相较于PE21的146.982 μg/m³，更接近观测值(118.607 μg/m³)。这一改善的核心原因是：PPM的高阶守恒特性降低了模型对稀疏站点的敏感度，即使在站点间距较大的区域，也能通过多项式重构维持浓度分布的物理合理性，避免了Cressman插值和传统DCIM方法不足导致的极值偏离。从算法机理来看，PPM方法之所以能缓解稀疏数据区域的过拟合问题，核心在于其高阶守恒特性与分段多项式重构逻辑：与Cressman插值基于距离加权的局部拟合不同，PPM通过在每个网格单元内构建二次多项式分布，同时满足单元间的质量守恒和梯度连续性，能够在有限观测点的约束下，合理推断未观测区域的浓度变化趋势，避免对单一观测值的过度依赖。具体而言，PPM通过计算单元边界的通量守恒(公式中${\langle C^n\rangle }_{i+1/2}$ 项)和内部梯度($C_{6,i}^n$ 项)修正，有效抑制了稀疏数据场景下的“局部极值放大”现象，这与2015年3月研究区域的气象背景高度契合：该时段区域扩散条件整体均匀，PM10浓度空间梯度以平缓变化为主，PPM的高阶平滑特性能够更好地匹配这种气象驱动下的浓度分布规律，而Cressman插值的距离加权逻辑则易受稀疏站点的局部异常值影响。

同时需要客观指出，本研究的站点分布稀疏性仍对PPM方法的性能产生了一定限制：在内蒙古东部部分无站点覆盖的区域，PE22的模拟浓度仍存在约15%的偏差，这表明尽管PPM的算法机理更适配稀疏数据场景，但当站点间距超过其有效插值尺度时，仍难以完全消除数据稀缺带来的不确定性。

Figure 6. PM10 concentration distribution obtained from experiments PE21 and PE22 ((a): PE21; (b): PE22)

图6. 实验PE21和PE22得到的PM10浓度分布((a): PE21; (b): PE22)

4. 结论与讨论

本研究提出的基于伴随模型的动态约束插值法(DCIM)，显著提高了中国中南部复杂地理区域PM10浓度的时空重建精度。该方法通过地面监测数据以及优化的高阶守恒插值(PPM)算法实现。借助伴随模型求解伴随方程以计算目标函数梯度，并结合迭代优化算法逼近最优解，该方法显著克服了传统插值法在观测稀疏区域的精度局限性。

从模拟精度和效果的角度出发，本研究将传统的Cressman插值法与基于Cressman的动态约束插值法进行了比较。结果表明，模拟值与观测值的对比显示，动态约束插值法有效降低了模拟误差，提高了精度。更重要的是，它生成了更平滑的空间分布，其核心优势在于通过物理模型约束减少了纯数学插值的主观性。这表明该方法能够有效地重建PM10的空间浓度分布，尽管在陕西、山西和内蒙古交界地区以及内蒙古东部出现了过拟合现象。进一步分析表明，过拟合可能源于观测数据不足。当Kriging插值法应用于大网格间距的数据集时，其插值误差往往无法满足要求。观测数据不足导致反演问题的不适定性，这是一个需要解决的关键问题。

为解决观测数据不足导致的不适定性问题，本研究进一步引入了高阶守恒插值(PPM)算法。将该算法应用于独立点方案，有效缓解了反演问题的不适定性。将优化后的模型与未优化的模型进行比较，结果表明前者在陕西、山西和内蒙古的边界地区以及内蒙古东部部分缓解了过拟合现象，这得益于PPM的高阶守恒特性与分段多项式重构逻辑，能够在有限观测约束下维持浓度分布的物理合理性，避免对局部观测值的过度依赖。这一结果验证了PPM算法的独特优势，即在观测数据稀少的情况下能提供更优的插值质量和稳健的性能。因此，通过地面监测数据以及优化的高阶守恒插值(PPM)算法可作为模拟PM10分布的有效研究方法，并为制定中国中南部地区的空气污染控制政策提供参考。并且本研究开发的DCIM方法不仅为PM10模拟提供了一种新颖的技术范例，还为优化算法同化提供了一个可扩展的理论框架。

值得注意的是，当前研究区域的站点分布稀疏性对本研究结论的可靠性与普适性造成了不可忽视的影响：

从空间约束有效性来看，研究区域内85个监测站的空间分布不均衡，中南部平原城市群(如长江中下游、珠江三角洲)站点间距约为20~30公里，而内蒙古东部、晋陕蒙边境等区域的站点间距普遍超过50公里，部分区域甚至达到80公里以上，远超本研究模型的有效插值尺度(约40公里)。这种不均衡的稀疏分布导致观测数据无法完整捕捉PM10浓度的空间变异特征，尤其是在无站点覆盖的“数据空白区”，模型只能依赖数学假设与邻域站点的间接推断，使得插值结果的物理真实性难以保证。

对方法性能评估的干扰方面，站点稀疏性直接影响了交叉验证结果的客观性。在四折交叉验证中，检验数据与同化数据的空间重叠度较低，部分检验点周围缺乏足够的同化站点作为约束，导致MAE、RMSE等误差指标可能被低估或高估，例如晋陕蒙边境区域的检验点，因周边仅1~2个同化站点，其模拟误差无法真实反映方法的实际性能，进而可能导致对DCIM结合PPM方法优势的判断存在偏差。

对结论普适性的限制尤为明显：本研究得出“PPM算法可有效缓解过拟合”的结论，主要基于春季短期数据与当前稀疏站点条件。但在站点更稀疏或稀疏模式不同的场景下(如西北干旱区、西南山区)，PPM算法的多项式重构逻辑可能因缺乏足够观测约束而失效，无法维持浓度分布的物理合理性；同时，站点稀疏导致本研究无法验证方法在不同污染梯度、不同地形条件下的适应性，使得“DCIM结合PPM方法适用于中国中南部PM10浓度模拟”的结论仅能局限于当前站点覆盖范围内，难以推广至更广泛的区域。

本研究存在以下局限性，需在未来研究中进一步完善：

1. 模型简化带来的不确定性：本研究采用二维对流扩散方程，忽略了垂直扩散过程与复杂光化学转化机制，而实际PM10的传输扩散是三维耦合过程，且二次PM10的生成机制会显著影响浓度分布，模型简化可能导致高海拔或工业排放集中区域的模拟偏差；2. 研究时段的局限性：实验数据仅涵盖2015年3月2日至5日的72小时，属于春季短期观测，未涵盖冬、夏、秋等不同季节的污染特征(如冬季供暖期PM10浓度峰值、夏季降水湿沉降影响)，导致方法的季节普适性尚未得到充分验证；3. 站点分布稀疏的影响：尽管PPM算法缓解了过拟合，但研究区域内部分区域(如内蒙古东部、晋陕蒙边境)的监测站点间距过大，仍导致插值结果存在一定偏差，稀疏站点无法提供足够的空间约束信息，是影响模拟精度的关键因素；4. 数据来源的单一性：仅利用地面监测数据进行同化与验证，缺乏卫星遥感(如MODIS AOD)、数值模型输出等多源数据的融合，可能导致观测误差被放大，且无法有效弥补地面站点的空间覆盖缺口。

为提升研究的鲁棒性与普适性，未来可从以下方面开展拓展研究：

1. 扩展研究时段至典型季节(如冬季1个月、夏季1个月)，涵盖不同污染水平与气象条件(如沙尘、暴雨、静稳天气)，系统验证DCIM结合PPM方法在不同场景下的适应性；2. 优化模型结构，引入三维对流扩散方程与二次PM10生成模块，提升模型的物理完整性；3. 融合多源数据(地面监测、卫星遥感、数值模型输出)构建多维数据集，通过数据同化技术减少单一数据源的误差，同时弥补地面站点的空间稀疏性。