大气污染溯源算法及其技术实现

污染溯源基础概念知识

大气污染溯源是指识别并追踪污染物的来源及其传输过程，以确定造成大气污染的根本原因和污染物传播路径的技术和方法。这对于制定有效的控制和减轻污染策略至关重要。大气污染的溯源主要涉及以下几个方面：

污染源识别：首先要确定污染源的类型，大气污染主要分为自然源（如火山爆发、沙尘暴、森林火灾产生的烟尘）和人为源（工业排放、汽车尾气、农业燃烧、燃煤电厂、建筑活动）。
化学与物理分析 - 通过对大气样本（如降水、颗粒物、气溶胶体、气体）进行化学成分分析，识别污染物的种类（如二氧化硫化物、氮氧化物、臭氧族化合物、重金属、挥发性有机化合物）。
同位素标记与分子标志物 - 使用稳定同位素和放射性同位素标记、分子标志物追踪特定污染物，比如13H、13C、14N等，可以揭示污染物的来源和传输路径。
模型模拟 - 利用大气传输模型模拟污染物的扩散过程，如高斯羽烟、WRF模型、CMAQS模型，结合气象数据（风向、温湿度、压力、风速、逆温层）模拟污染物的传播轨迹。
遥感技术 - 利用卫星遥感、无人机、航空遥感图像分析大气污染状况，监测大范围的污染分布，如PM2.5、NOx、O3浓度。
时空分析 - 结合时间序列分析和空间统计方法，如GIS分析，追踪污染变化趋势和空间分布，识别热点区域。
数据融合与大数据 - 利用多源数据（环境监测站、移动监测、社会经济数据）融合分析，通过AI、机器学习算法提高溯源精度。大气污染溯源是一项复杂的跨学科工作，需要环境科学、化学、气象学、计算机科学、地理信息科学等多领域知识的综合应用。通过这些技术手段，可以准确找到污染源头，为政策制定和环保措施提供科学依据，减少大气污染对人体健康和环境的负面影响。

大气污染溯源流程

2.1站点污染物超标异常报警监测

大气污染物超标监测是环境监测中一个关键环节，它涉及到对空气中存在的各种有害物质浓度进行定期或连续监测，以确保其浓度保持在国家或地方规定的安全限值内。当监测结果显示污染物浓度超过这些限值时，即认为是超标。以下是对大气污染物监测超标处理的一般流程：

实时监测与数据收集：使用在线监测设备（如自动监测站、传感器网络）连续监测大气中的主要污染物，如二氧化硫氧化物（SOx）、氮氧化物（NOx）、颗粒物（PM2.5、PM10）、臭氧、异味等，实时收集数据。
数据验证与分析：对收集到的数据进行质量控制和校验，确保数据的准确性和可靠性。运用统计分析、趋势分析、时空分布等方法，识别超标原因。

3. 超标预警：设定阈值预警系统，一旦监测数据超过设定的阈值，立即触发预警，可视化展示实时高值预警站点。

2.2融合气象条件溯源

气象条件下的污染物高斯扩散原理：

环境科学中的高斯扩散模型在环境科学领域，高斯模型常用于模拟污染物在大气、水体或土壤中扩散的过程。高斯模型假设污染物的浓度分布遵循高斯分布，其中污染物的浓度随着距离源的距离增加而呈指数下降。这通常用于预测污染物的浓度场、评估环境影响、规划应急响应或制定污染控制策略。

统计学中的高斯分布在统计学中，高斯分布是一种连续概率分布，常用于描述自然界中许多随机变量的分布，如误差或测量误差。高斯分布由均值（μ）、方差（σ²）决定，公式为: [ f(x)=\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}]在环境模型中，高斯分布用来描述污染物浓度的概率分布，即在某一特定点的污染物浓度值的概率，随着距离污染源的距离增加而按高斯分布规律衰减小。

高斯扩散方程在数学上，高斯扩散方程描述了这种浓度随时间变化的扩散过程，它是一个偏微分方程，通常形式为： [ \frac{\partial c}{\partial t} = D \nabla \Delta c] 其中 (c) 是浓度，(t) 是时间，(D) 是扩散系数，(\Delta) 是拉普朗克算子，代表空间二阶导数。编程语言（如.NET、Java等）中实现高斯模型，或是在GIS、气象、大气科学、医疗、警务等领域的应用；

示例应用验证

2.3污染溯源查找

风速小于1.5m/s查找周边25公里范围内企业

风速1.5m/s以上查找上下风口22.5度范围内企业

2.4相关性分析

相关性分析是统计学中一种重要的方法，用于研究两个或多个变量之间相互关系的性质，主要包括关系的方向（正相关或负相关）、强度以及是否存在线性。相关性分析不意味着因果关系，而是表明变量之间的伴随关系或关联度量度。以下是相关性分析的一些核心概念和方法：

相关系类型

•正相关：**当一个变量增加时，另一个变量也倾向于增加，如身高和体重。

•负相关：**一个变量增加时，另一个变量倾向于减少，如冰淇淋销量与气温。

•无相关：**两个变量变化无明显关联，或关系随机。

相关系强度

相关系数值（通常用r表示）衡量相关强度：

•|r ≈ 0：几乎无相关

•0.3 ≤ |r| < .5：低度相关

•.5 ≤ |r| < .8：中度相关

•|r| ≈ 1：高度相关

见方法

•Pearson相关系数：用于度量连续变量，假设数据近似正态分布。

•Spearman秩相关：连续或分类变量，不关心变量分布，关注排序关系。

•Kendall相关：非参数，适用于小样本量数据，排序数据。

•偏相关：控制其他变量影响，分析两变量间净相关。

应用

•气象：影响分析。

分析过程

数据收集：确保数据质量，代表性。
预处理：缺失值处理，异常值处理。
选择模型：根据数据类型选相关系数。
计算：软件执行分析。
解读：图表，理解r值。

相关性分析是理解和预测模型、决策支持科学假设的重要工具。

关联企业和企业站点进行溯源相关性分析（推最近26小时，相关性分析连续24小时）

2.4集成AI精准溯源

后台代码算法实现

try {Map<String, Object> kv = new HashMap<>();kv.put("labelMap", new HashMap<>());Map<String, String> headerMap = new HashMap<>();if(authorization==null)authorization="eyJhbGciOiJIUzUxMiJ9.eyJ1c2VyX2lkIjoxLCJ1c2VyX2tleSI6IjlkNGZjN2VjLWRkMDEtNDE1MC1iMjRjLTZlMDFhM2EzOWMzYSIsInVzZXJuYW1lIjoiYWRtaW4ifQ.Yoxhd_UCCMiBPHHsFaFAZp2q0kQZ7guXR6oTA4wz3CNEYKIOZ1U-IOE2pxDzEzVG2N-CWI8S4arj3ffNXF5V8g";headerMap.put("Authorization", authorization);String[] confs= FileUtils.ReadTxt("peConfig.txt").split(";");JSONArray datalist = null;String url = confs[0] + "/stationExceed/changingTrends?contaminant=" + factor + "&end=" + end + "&start=" + start + "&stationIds=" + id;JSONObject obj = JSON.parseObject(HttpUtils.get(url,kv, headerMap));logger.warn("请求的地址为："+url+"返回的数据为"+obj);datalist = (JSONArray)((JSONObject)((JSONArray)(((JSONObject)obj.get("data")).get("series"))).get(0)).get("values");double[] x = new double[datalist.size()-2];for (int i = datalist.size()-1,j=0; i > 1; i--,j++) {x[j] = TransferDouble(datalist.get(i));}String entIds = "";JSONArray datalist1 = null;String url1 = confs[0] + "/station/list?stationTypeDictIds="+stationTypeDictIds;JSONObject obj1 = JSON.parseObject(HttpUtils.get(url1,kv, headerMap));
//            logger.warn("请求的地址为："+url1+"返回的数据为"+obj1);datalist1 = (JSONArray)obj1.get("data");JSONArray eArr = new JSONArray();String thetas = "";if(windSpeed>=1.5){for(int i = 0; i < datalist1.size(); i++){double distancel=InterpolationUtils.GetDistance(lon,lat,TransferDouble(((com.alibaba.fastjson.JSONObject)datalist1.get(i)).get("lon")),TransferDouble(((com.alibaba.fastjson.JSONObject)datalist1.get(i)).get("lat")));if(((com.alibaba.fastjson.JSONObject)datalist1.get(i)).get("id").equals(id)||distancel>100000){continue;}double dgr=InterpolationUtils.getDegrees(lon,lat,TransferDouble(((com.alibaba.fastjson.JSONObject)datalist1.get(i)).get("lon")),TransferDouble(((com.alibaba.fastjson.JSONObject)datalist1.get(i)).get("lat")));if(dgr <0 ){dgr += 360.0;}if(wd>=11.25&wd<=349.75){if(Math.abs(dgr-wd)<=11.25){eArr.add(datalist1.get(i));entIds = entIds + ((JSONObject)datalist1.get(i)).get("id").toString() + ",";thetas=thetas+String.valueOf(Math.abs(dgr-wd))+",";}}else {if(dgr<11.25){dgr += 360.0;}if(wd<11.25){wd += 360.0;}if(Math.abs(dgr-wd)<=11.25){eArr.add(datalist1.get(i));entIds = entIds + ((JSONObject)datalist1.get(i)).get("id").toString() + ",";thetas=thetas+String.valueOf(Math.abs(dgr-wd))+",";}}if(eArr.size()==5){break;}}int ecount = eArr.size()+5;for(int i = 0; i < datalist1.size(); i++){double distancel=InterpolationUtils.GetDistance(lon,lat,TransferDouble(((com.alibaba.fastjson.JSONObject)datalist1.get(i)).get("lon")),TransferDouble(((com.alibaba.fastjson.JSONObject)datalist1.get(i)).get("lat")));if(((com.alibaba.fastjson.JSONObject)datalist1.get(i)).get("id").equals(id)||distancel>100000){continue;}double dgr=InterpolationUtils.getDegrees(lon,lat,TransferDouble(((com.alibaba.fastjson.JSONObject)datalist1.get(i)).get("lon")),TransferDouble(((com.alibaba.fastjson.JSONObject)datalist1.get(i)).get("lat")));if(dgr <0 ){dgr += 360.0;}double wd1=wd;if(wd>180){wd1=wd-180;}else {wd1=wd+180;}if(wd1>=11.25&wd1<=349.75){if(Math.abs(dgr-wd1)<=11.25){eArr.add(datalist1.get(i));entIds = entIds + ((JSONObject)datalist1.get(i)).get("id").toString() + ",";thetas=thetas+String.valueOf(Math.abs(dgr-wd))+",";}}else {if(dgr<11.25){dgr += 360.0;}if(wd1<11.25){wd1 += 360.0;}if(Math.abs(dgr-wd1)<=11.25){eArr.add(datalist1.get(i));entIds = entIds + ((JSONObject)datalist1.get(i)).get("id").toString() + ",";thetas=thetas+String.valueOf(Math.abs(dgr-wd1))+",";}}if(eArr.size()==ecount){break;}}}else {for(int i = 0; i < datalist1.size(); i++){if(((com.alibaba.fastjson.JSONObject)datalist1.get(i)).get("id").equals(id)){continue;}double dgr=InterpolationUtils.getDegrees(lon,lat,TransferDouble(((com.alibaba.fastjson.JSONObject)datalist1.get(i)).get("lon")),TransferDouble(((com.alibaba.fastjson.JSONObject)datalist1.get(i)).get("lat")));double ds=InterpolationUtils.GetDistance(lon,lat,TransferDouble(((com.alibaba.fastjson.JSONObject)datalist1.get(i)).get("lon")),TransferDouble(((com.alibaba.fastjson.JSONObject)datalist1.get(i)).get("lat")));if(ds <25000 ){eArr.add(datalist1.get(i));entIds = entIds + ((JSONObject)datalist1.get(i)).get("id").toString() + ",";thetas=thetas+String.valueOf(Math.abs(dgr-wd))+",";}if(eArr.size()==10){break;}}}for(int i = 0; i < datalist1.size(); i++){if(((com.alibaba.fastjson.JSONObject)datalist1.get(i)).get("id").equals(id)){eArr.add(datalist1.get(i));entIds = entIds + ((JSONObject)datalist1.get(i)).get("id").toString() + ",";thetas=thetas+"0,";break;}}if(entIds.equals("")){return "";}JSONArray datalist2 = null;String url2 = confs[0] + "/stationExceed/changingTrends?contaminant=" + factor + "&end=" + end + "&stationIds=" + entIds.substring(0, entIds.length() - 2) + "&start=" + start;JSONObject obj2 = JSON.parseObject(HttpUtils.get(url2,kv, headerMap));logger.warn("请求的地址为："+url2+"返回的数据为"+obj2);datalist2 = (JSONArray)((JSONObject)obj2.get("data")).get("series");double[] pearsonCorrelations=new double[datalist2.size()];double[] emissionCorrelations=new double[datalist2.size()];for(int i=0;i<datalist2.size();i++){JSONArray datalistL=(JSONArray)(((JSONObject)datalist2.get(i)).get("values"));double[] y1 = new double[datalistL.size()-2];double[] y2 = new double[datalistL.size()-2];double[] y3 = new double[datalistL.size()-2];for (int k = datalist.size()-1,j=0; k > 1; k--,j++) {y1[j] = TransferDouble(datalistL.get(k));}for (int k = datalist.size()-2,j=0; k > 0; k--,j++) {y2[j] = TransferDouble(datalistL.get(k));}for (int k = datalist.size()-3,j=0; k > -1; k--,j++) {y3[j] = TransferDouble(datalistL.get(k));}double correlation1 = StatisticsUtil.correlation(x, y1, PEARSON_ID);double correlation2 = StatisticsUtil.correlation(x, y2, PEARSON_ID);double correlation3 = StatisticsUtil.correlation(x, y3, PEARSON_ID);if(Double.isNaN(correlation1)){correlation1 = 0;}if(Double.isNaN(correlation2)){correlation2 = 0;}if(Double.isNaN(correlation3)){correlation3 = 0;}pearsonCorrelations[i]=correlation1;double q = y1[0];if(correlation2>pearsonCorrelations[i]){pearsonCorrelations[i]=correlation2;q = y2[0];}if(correlation3>pearsonCorrelations[i]){pearsonCorrelations[i]=correlation3;q = y3[0];}double lz=Math.sqrt(Math.pow(((JSONObject)eArr.get(i)).getDouble("lon")-lon,2)+Math.pow(((JSONObject)eArr.get(i)).getDouble("lat")-lat,2));double lX=lz*Math.cos(TransferDouble(thetas.split(",")[i])/180);double lY=lz*Math.sin(TransferDouble(thetas.split(",")[i])/180);emissionCorrelations[i]=InterpolationUtils.interpolation_GaussPlumeP(lX, lY,1, windSpeed, q, 0,"A");((JSONObject) eArr.get(i)).put("emissionCorrelation",emissionCorrelations[i]*1000);((JSONObject) eArr.get(i)).put("pearsonCorrelation",pearsonCorrelations[i]);((JSONObject) eArr.get(i)).put("lineGeoField","[["+((JSONObject)eArr.get(i)).getString("lon")+","+((JSONObject)eArr.get(i)).getString("lat")+"],["+lon+","+lat+"]]");}return eArr.toJSONString();}catch (Exception e){return "";}