一种复合受体模型污染源解析方法与流程

本发明属于环境保护技术领域，具体涉及一种复合受体模型污染源解析方法。

背景技术：

在环境保护领域，当进行环境质量评估、环境风险评价、污染物削减时，需要通过定量方法筛选出特征污染物并解析出污染源及其贡献率，从而切实有效地控制污染，保障环境安全。所谓污染源解析，是研究污染源对周围环境污染影响和作用的一系列方法，常用的解析技术有清单分析法、扩散模型分析法和受体模型分析法。其中，清单分析法是通过观测和模拟污染物的源排放量、排放特征及排放地理分布等,建立列表模型进行污染源解析；扩散模型分析法属于预测式模型，它是通过输入各个污染源的排放数据和相关参数信息来预测污染物的时空变化情况；受体模型分析法则通过对受体样品的化学和显微分析，确定各污染源贡献率的一类技术，其最终目的是识别对受体有贡献的污染源，并且定量计算各污染源的分担率。

图1为现有技术中扩散模型分析过程示意图。如图1所示，利用扩散模型进行污染源解析，不仅需要完整的排放源清单，还要知道排放源的排污特征，排放率，排放时间，受体位置等排放参数。可以看出，利用扩散模型进行源解析成本较高，难度较大。

图2为现有技术中受体模型分析过程示意图。如图2所示，受体模型的优势在于它利用真实的受体数据来驱动模型，能够对每一个样品进行定量源解析，也可以根据适应度分析来识别模型源与真实源的误差。其劣势在于模型要有受体样品浓度，甚至还要求有与受体浓度相对应的不确定性值，要求源线独立，不能有效区分外来源和本地源。受体模型开展污染源解析，需要满足几个基本的假设条件：(1)源线固定；(2)受体数据获取过程及条件没有明显的变化；(3)受体数据具有代表性；(4)污染源存在时间或空间变异。与扩散模型相比，受体模型的输入数据来源较为容易，不需要完整的源谱清单，模型运转成本较低，同时，作为扩散模型的有效补充，可以使人们通过对偏差的认识理解源和受体的关系。

相对而言，受体模型不依赖于排放源排放条件及水文水质等数据,不用追踪污染物的迁移过程，且与预测式的扩散模型相比，受体模型算是一种诊断式模型，它解释过去而不是未来，比较适合用来实现污染源的源解析。根据是否需要事先掌握详尽的污染源排放特征谱，受体模型可分为化学质量平衡模型(CMB)、多元统计模型、UNMIX模型等。随着环境监测网络的日益完善，监测手段、方法和技术的不断进步，研究者可以取得大量而精准的观测样品数据，这为开展污染源解析的定量研究提供了保障，受体模型源解析技术得到了较快的发展和应用，被广泛用来开展大气、水体、沉积物和土壤等环境介质中各种污染物的污染源解析研究。

其中，化学质量平衡模型是一类发展较为成熟的源解析模型，它不依赖于排放源条件、气象、地形等数据，不用追踪颗粒物的迁移过程，模型原理清楚，易于接受，是应用最为广泛的一种受体模型。它对于无机污染物的源解析已经较为成功，被美国环保局定为无机物源解析的首选方法。该模型依据质量守恒定律建立污染源定量解析模型，通过测量源与受体样品的物理、化学性质，定性识别对受体有贡献的污染源并定量计算各污染源的分担率，其理论依据是各个污染源的指纹谱是有一定差别的，从而可以通过检测受体中的各种物质的含量(组成)来确定各污染源的贡献率。但化学质量平衡模型应用需要满足多项假设：(1)可以识别出对环境受体中有明显贡献的所有污染源类，并且各源类所排放的污染物化学组成有明显的差别；(2)各源类所排放的污染物化学组成相对稳定，化学组分之问无明显影响；(3)各源类所排放的污染物之间没有相互作用，在传输过程中的变化可以被忽略；(4)所有污染源成分谱是线性无关的；(5)污染源种类少于或等于化学组分种类；(6)测量的不确定度是随机的、符合正态分布。更为重要的是，污染源指纹谱的获取难度较大，成本较高，需要根据环境的不断变化，追踪排放源清单，更新排放特征指纹库，特别是在当前监测网络并不是覆盖所有污染排放点的现实情况下，污染源指纹谱已知的条件常常难以满足。

与化学质量平衡模型不同的是，多元受体模型不需要事先已知具体的污染源成分谱，而是通过因子分解来提取污染源指纹谱，并获得各污染因子的源贡献率。多元受体模型可以解析数据集合，压缩数据维数，分析多个变量之间的关系，这个特性对于大量观测数据的分析极为有用，可以使用较少的有代表性的因子来说明众多变量的主要信息。然而，多元受体源解析模型也存在一些缺陷，例如，因子分解不能得到“唯一”的因子载荷和因子得分，通过因子分析提取的因子载荷和因子得分可以通过正交旋转得到多种形式的“不适定解”，更重要的利用因子分析得到的因子载荷和因子得分中常常出现负值，这与实际应用不相符，因为负值的源成分谱或源贡献率没有实际的物理意义；正定矩阵分解模型尽管考虑了受体数据的不确定性问题，但通过引入惩罚因子和辅助因子来保证提取的污染源贡献率也没有考虑污染源排放的不确定性问题。

技术实现要素：

本发明要解决的技术问题是提供一种复合受体模型污染源解析方法，用以解决现有技术中污染源解析方法的不足，提高源解析结果的准确性和可靠性。

根据本发明的一个方面，提供了一种复合受体模型污染源解析方法，所述方法包括：

步骤S1，受体样本数据预处理；

步骤S2，提取主成分因子数；

步骤S3，分解非负约束矩阵；

步骤S4，求解污染源贡献率。

上述方案中，所述步骤S1进一步包括：

步骤S101，从数据处理的角度，审核受体样本数据；

步骤S102，从污染物的角度，选择复合受体模型变量；

步骤S103，对所述模型变量进行标准化处理。

上述方案中，所述从污染物的角度选择复合受体模型变量，具体包括如下过程：

第一，不能同时选择活性强，或存在反应的污染物；

第二，受污染事故影响而检出的污染物不能用于模型计算；

第三，不能同时选择含有相同元素，且迁移过程中会发生相互转换的污染物；

第四，浓度未检出或数据缺失的污染物，通过数据替换并设置合适的不确定性值；

第五，污染物均有浓度值但数据质量未知，通过再分析以确定是否排查该污染物。

上述方案中，所述步骤S2中主成分因子特征为：

特征值大于1；

累积方差贡献率值大于85％；

决定系数大于0.9；

Exter函数值小于0.1。

上述方案中，所述步骤S2进一步包括：

步骤S201，求解协方差矩阵Z：Z＝DD^t，D为受体数据集；

步骤S202，求解特征值E及特征向量Q：Z＝QEQ^t；

步骤S203，求解无旋转的因子载荷矩阵S：S＝QE^1/2；

步骤S204，求解无旋转的因子得分矩阵R：R＝(S^tS)^-1S^tD；

步骤S205，生成特征值、累积方差、决定系数和Exner函数判据矩阵，提取主成分因子数。

上述方案中，所述步骤S3中分解非负约束矩阵，通过改进的交替回归法进行分解。

上述方案中，所述改进的交替回归法，进一步包括：

步骤S301，k＝0，设置两个未知矩阵的初始值S(k)＝S₀和R(k)＝R₀；

步骤S302，保持R＝R(k)不变，寻求合适的ΔS，使得S＝S(k)+ΔS最小化Q值(D-SR)；

步骤S303，保持S＝S(k)+ΔS不变，寻求合适的ΔR，使得R＝R(k)+ΔR最小化Q值；

步骤S304，寻求合适的扩展系数α，使得R＝R(k)+α·ΔR，S＝S(k)+α·ΔS最小化Q值；

步骤S305，k＝k+1，重复步骤(b)～(c)直至迭代收敛。

上述方案中，所述步骤S4，进一步包括：

首先假定某污染源中污染物的标准偏差线性依赖于该污染源中指定污染物的贡献比例，即有：

其中：σ_i表示提取出的源成分谱中第i类污染物的标准偏差；δ_i是受体样本中第i类污染物的标准偏差；表示提取出的源成分谱中第i类污染物的排放浓度；D_i表示受体样本中第i类污染物的浓度；

其次，计算源贡献率过程如下：

步骤S401，假定初始源贡献率其中n表示污染源数；

步骤S402，计算有效方差对角阵其元素为：

其中：为第i种污染物的不确定性标准偏差；为第j类污染源的第i种污染物浓度的不确定性标准偏差；为上次求得的源贡献率；

步骤S403，源贡献率计算：

其中：S为污染源指纹图谱；ST为其转置矩阵；为步骤(b)得到的有效方差对角阵；D为样品浓度；

步骤S404，源贡献率迭代偏差计算：

如果偏差大于某设定误差精度(如0.001)，则重复步骤S401～403；否则，终止迭代。

上述方案中，所述计算源贡献率还可以包括：

步骤S405，通过

计算源贡献率不确定性偏差。

本实施例基于复合受体模型的污染源解析方法，将两类受体模型进行复合应用，提出了一种基于非负约束矩阵分解化学质量平衡复合模型的污染源解析方法。本发明利用非负约束矩阵分解反推出线性无关的污染源指纹图谱，以满足化学质量平衡模型对污染源成分谱的线性独立要求，并利用CMB模型计算出污染源贡献率，从而充分发挥两种模型各自的优势，着力提高源解析结果的准确性和可靠性。本发明可以为环境质量评估、环境风险评价、水源保护、总量减排、生态修复、污染事故调查、损害赔偿提供必要的科技支撑，对于有效地控制环境污染，保障生态安全具有重要的现实意义。本发明成果提出的源解析方法耦合应用的思路，对污染源解析理论与方法的完善也具有一定的科学价值。

附图说明

图1是现有技术中扩散模型污染源解析过程示意图；

图2是现有技术中受体模型污染源解析过程示意图；

图3是本发明实施例复合受体模型污染源解析方法过程示意图；

图4是本发明实施例非负约束矩阵分解化学质量平衡复合受体模型算法逻辑图。

具体实施方式

为使本发明要解决的技术问题、技术方案和优点更加清楚，下面将结合附图及具体实施例进行详细描述。

通过对现有技术的分析可知，基于单一的受体模型解析技术进行污染源解析，存在很多不足。例如，化学质量平衡模型是应用比较成熟、解析比较精准的污染源识别方法，但其应用条件比较严格，除了需要预先知道污染源指纹图谱外，各污染源之间还不能出现“混源”现象，这往往与实际应用条件不相符合。而多元受体模型可以反推出污染源指纹图谱，求得因子得分，但完全基于因子分析的源成分谱和源贡献率往往存在解的“不适定性”问题，需要通过增加附加条件(如非负约束、惩罚函数、R项等等)以压缩解空间，即便如此，利用因子分解不能同时考虑受体浓度和污染源排放浓度的不确定性问题，弱化了源贡献率计算的精度与可靠性。

针对单一受体模型的污染源解析方法的局限性，本发明提出了一种复合受体模型的污染源解析方法，即非负约束矩阵分解化学质量平衡复合模型(Nonnegative Matrix Factorization Chemical Mass Balance,简称为NMFCMB)。该模型利用非负约束矩阵分解得到线性无关的污染源排放成分谱及对应的不确定性值，以满足化学质量平衡模型对污染源成分谱的线性独立要求，再利用化学质量平衡模型计算出相应的污染源贡献率，从而充分发挥CMB模型和因子分析模型的各自优势，达到污染源解析的最优化处理。

图3为本发明具体实施例的复合受体模型污染源解析方法的过程示意图。如图3所示，本实施例的复合受体模型污染源解析方法，包括如下步骤：

步骤S1，受体样本数据预处理；

步骤S2，提取主成分因子数；

步骤S3，分解非负约束矩阵；

步骤S4，求解污染源贡献率。

图4为本实施例非负约束矩阵分解化学质量平衡复合受体模型算法逻辑图。如图4所示：

所述步骤S1具体包括：

步骤S101，从数据处理的角度，审核受体样本数据。

首先需要说明的是，这里的受体样品数据，主要指污染物样品的浓度测量数据，包括未检出项、缺失项、浓度测量值等。因此，本实施例还可以包括步骤S0，收集样本浓度数据。这里的样本浓度数据可以表达为D＝m×n，D为样本数据矩阵。

不确定性分析包括模型本身的不确定性和模型输入的不确定性。本发明中的复合受体模型，利用足够的受体样本集，利用多元统计技术，提取污染源成分谱，计算源贡献。很明显，作为模型输入项的受体样本数据对于模型计算结果的精度有着重要的意义。应用中，受体样本数据集可以是时间和/或空间序列的污染物浓度矩阵。而无论是时间还是空间序列的污染物浓度集，其数据质量总是会受采样、测量、实验室分析、数据记录(计算机录入)、数据拷贝等过程的影响。例如，不同采样时间采集的样品，可能受环境因素、排放周期特征影响出现较大的变异。另外，当前的监测网络(能力)存在较大的区域差异性，自动监测能力还不是很强，大部分的常规监测还是依靠人工方式开展。因此，测量仪器，分析方法，数据记录，数据拷贝也都可能带来受体浓度值的不确定性误差。再有，受体样本中的污染物来自于不同污染源的贡献，而污染源排放的污染物从源到受体过程也存在一定的不确定性，如物理、化学、生物降解等。作为受体模型输入项的受体样本，其数值存在三种情况：测量浓度值、未检出(低于检出限)和空值(数据丢失或没有测量)。

因此，这里的数据审核，包括如下三个方面：

对于有浓度测量值的情况，根据环境调查和已有经验对污染物浓度测量值进行识别，判断是否存在异常现象，其不确定性值考虑分析不确定性和方法检出限；所述数据审核，是对未检出项进行处理；对无测量值项的处理。

对于未检出项，用检出限的一半替代，而对应的不确定性则取值为检出限的5/6倍；

对于空值，处理方式可以通过下述三种方式实现：第一，剔除空值样本；第二，包含空值的污染物样品数据变量不纳入模型计算；第三，以污染物浓度测量值算术平均或几何平均代替空值，并赋予平均浓度三到四倍的不确定性值。

步骤S102，从污染物的角度，选择复合受体模型变量。

本发明中基于复合受体模型进行污染源解析时，其模型变量指的是不同污染物的浓度值，因此选择合适污染物的浓度值作为模型变量对于模型计算结果有着重要的影响。

所述选择合适的污染物，包括：

第一，根据复合受体模型的基本假设，污染物从源到受体的过程不能相互之间发生反应，因此不能同时选择那些活性较强，或存在反应的污染物；

第二，由于本发明的复合受体模型假定，不同时间采集的受体样本，其污染浓度受相同污染源的影响，因此，受污染事故影响而检出的污染物不能用于模型计算；

第三，对于含有相同元素，且迁移过程中会发生相互转换的污染物(即污染指标)不能同时选择，如氨氮、总氮的浓度值不能同时选择作为模型变量；

第四，浓度未检出或数据缺失的污染物，通过数据替换并设置合适的不确定性值减少其对模型结果的影响；

第五，污染物均有浓度值(测量值或预处理值)，但数据质量未知，需要通过深入分析以确定是否排查该污染物。

步骤S103，对模型变量进行标准化处理。

为解决不同污染物浓度的量纲不匹配问题，消除分析过程的偏差，需要对所选择的样品数据即模型变量进行标准化处理。本发明数据标准化可以通过平均值标准化方法或指数标准化方法进行，并得到标准数据。如：

所述步骤S2，提取主成分因子数，具体包括：

本实施例主成分因子数可根据特征值、累积方差贡献率、决定系数、Exner函数、自由度来确定。首先确定主成分因子。本实施例选取的主成分因子特征包括：特征值大于1的所有因子作为主成分因子；如果累积方差贡献率值大于85％，就可以认为因子分析的结果解释了足够多的变量信息；决定系数大于0.9的可以确定为主成分因子。特别的，对于Exner函数值，理论上可以从0到无限大，其值为0.5为可接受的最高值，小于0.1则认为结果很好。

提取主成分因子数的过程如下：

步骤S201，求解协方差矩阵Z：Z＝DD^t；

步骤S202，求解特征值E及特征向量Q：Z＝QEQ^t；

步骤S203，求解无旋转的因子载荷矩阵S：S＝QE^1/2；

步骤S204，求解无旋转的因子得分矩阵R：R＝(S^tS)^-1S^tD；

步骤S205，生成特征值、累积方差、决定系数和Exner函数判据矩阵，提取主成分因子数。

所述步骤S3，分解非负约束矩阵，进一步包括：

现有技术中，因子分解最普遍的算法就是交替回归最小二乘法，该算法首先假定其中某一个目标矩阵已知，寻求另一个未知矩阵使得目标函数值最小；然后再假定另一个目标矩阵已知，重复上述过程；如此反复直到趋于收敛。很明显，该算法简单，但效率低，计算量较大。为了改进交替回归法的计算效率，本实施例采取一种改进的交替回归法，所述改进的交替回归算法，通过实现两个未知矩阵的同步递增，从而达到提升收敛效率的目标，其基本过程如下：

步骤S301，k＝0，设置两个未知矩阵的初始值S(k)＝S₀和R(k)＝R₀；

步骤S302，保持R＝R(k)不变，寻求合适的ΔS，使得S＝S(k)+ΔS最小化Q值；

步骤S303，保持S＝S(k)+ΔS不变，寻求合适的ΔR，使得R＝R(k)+ΔR最小化Q值；

步骤S304，寻求合适的扩展系数α，使得R＝R(k)+α·ΔR，S＝S(k)+α·ΔS最小化Q值；

步骤S305，k＝k+1，重复步骤(b)～(c)直至迭代收敛。

所述步骤S4，求解污染源贡献率，进一步包括：

步骤S41，利用化学质量平衡模型计算污染源贡献率需要提供受体样品浓度及其不确定性，以及污染源排放成分谱(即指纹谱)及其不确定性。在本实施例的NMFCMB模型中，利用非负约束矩阵分解反推出的污染源成分谱的不确定性值计算方法如下：假定某污染源中污染物的标准偏差线性依赖于该污染源中指定污染物的贡献比例，即有：

其中：σ_i表示提取出的源成分谱中第i类污染物的标准偏差；δ_i是受体样本中第i类污染物的标准偏差；表示提取出的源成分谱中第i类污染物的排放浓度；D_i表示受体样本中第i类污染物的浓度。

步骤S42，源贡献率计算采用有效方差最小二乘法计算源贡献率，其步骤如下：

步骤S401，假定初始源贡献率其中n表示污染源数。

步骤S402，计算有效方差对角阵其元素为：

其中：为第i种污染物的不确定性标准偏差；为第j类污染源的第i种污染物浓度的不确定性标准偏差；为上次求得的源贡献率。

步骤S403，源贡献率计算：

其中：S为污染源指纹图谱；ST为其转置矩阵；为步骤(b)得到的有效方差对角阵；D为样品浓度。

步骤S404，源贡献率迭代偏差计算：

如果偏差大于某设定误差精度(如0.001)，则重复步骤S401～403；否则，终止迭代。

步骤S405，计算源贡献率不确定性偏差：

本实施例还可以包括步骤S5，进行可靠性检验。

具体的，本实施例中可以通过残差方差、方差均值、质量百分比、T统计值中的一种或多种进行可靠性检验，从而获得最终的源贡献率。

由以上技术方案可以看出，本实施例基于复合受体模型的污染源解析方法，将两类受体模型进行复合应用，提出了一种基于非负约束矩阵分解化学质量平衡复合模型的污染源解析方法。本发明利用非负约束矩阵分解反推出线性无关的污染源指纹图谱，以满足化学质量平衡模型对污染源成分谱的线性独立要求，并利用CMB模型计算出污染源贡献率，从而充分发挥两种模型各自的优势，着力提高源解析结果的准确性和可靠性。本发明既可以为环境质量评估、环境风险评价、水源保护、总量减排、生态修复、污染事故调查、损害赔偿提供必要的科技支撑，对于有效地控制环境污染，保障生态安全具有重要的现实意义。本发明成果提出的源解析方法耦合应用的思路，对污染源解析理论与方法的完善也具有一定的科学价值。

以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明所述原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。