两年前我第一次接到烟气扩散模拟的需求时甲方要求一周内给出污染物的空间分布图。我当时的思路还是传统 Java 后端那套写个循环、算个表格、塞到数据库再让前端去画。结果所有数据交上去以后对方 GIS 工程师问了一句你给的坐标是 WGS84 还是 CGCS2000我当场愣住了。那次之后我才意识到把高斯羽烟模型用 Java 落地难的不是公式而是从污染气象学参数到 GIS 空间网格的整条链路如何打通。这篇文章就把这条链路拆开讲清楚高斯羽烟公式怎么用 Java 实现、浓度场怎么转成等值面格网点、最后怎么输出成 GIS 能直接加载的图层。1. 高斯羽烟模型不只是公式GIS 集成前先搞懂它的适用边界1.1 标准浓度公式里每个符号到底代表什么先放下代码把模型本身说透。高斯羽烟扩散模型Gaussian Plume Model是描述连续点源排放的污染物在下风向空间浓度分布的一种半经验模型它不求解复杂的湍流方程而是假设污染物浓度在垂直和水平方向都服从正态分布。对绝大多数项目来说这个精度已经够用了。最常用的无障碍连续点源地面反射公式长这样C(x, y, z) Q / (2 * π * u * σ_y * σ_z) * exp(-y² / (2 * σ_y²)) * [exp(-(z - H)² / (2 * σ_z²)) exp(-(z H)² / (2 * σ_z²))]逐个拆开看C(x, y, z)空间点 (x, y, z) 处的污染物浓度单位是 g/m³ 或者 μg/m³Q污染源排放速率单位是 g/s这是模型里最重要的一个输入参数u烟囱出口高度处的平均风速单位 m/sσ_y、σ_z水平方向和垂直方向的扩散参数单位 m它们随下风向距离 x 增大而增大H有效源高单位 m等于烟囱几何高度加烟气抬升高度很多项目直接拿烟囱高度当有效源高会导致地面浓度被严重低估y计算点离烟羽中心线的水平横向距离风轴方向为 xz计算点离地面的高度m。公式里有两个 exp 项第二个 exp 就是所谓的“地面反射项”。污染物往下扩散碰到地面后不会消失而是会被反射回大气等效于在 z -H 的位置有一个虚像源在同时排放。如果你只计算地面浓度z 0这两个指数项就变成 2 * exp(-H² / (2 * σ_z²))公式会简化不少。GIS 项目里犯得最多的错误就是把 (x, y) 直接当成经纬度去算。x 是下风向距离y 是侧风向距离单位是米不是经纬度。建模第一步是先在以污染源为原点的笛卡尔坐标系里算出格网点相对坐标再用 GIS 的坐标转换把它换成实际地图坐标这个顺序不能反。1.2 扩散参数 σ_y、σ_z 怎么取查表还是公式拟合扩散参数是整个模型里最看经验的地方。国内外项目最常用的是 Pasquill-Gifford 曲线和 Briggs 公式。Pasquill 模型把大气稳定度分成 A、B、C、D、E、F 六类A 类最不稳定对应强太阳辐射、小风速F 类最稳定通常对应晴朗夜间。如果你没有本地实测数据可以直接用 Briggs 开放乡村型公式这套参数在平坦地形下通用性不错σ_y a * x^bσ_z c * x^d其中 x 是下风向距离a、b、c、d 是跟稳定度类型有关的系数。我给一个常用的系数表方便你直接塞进 Java 代码里稳定度abcdA0.220.8940.200.801B0.160.8940.120.801C0.110.8940.080.801D0.080.8940.060.801E0.060.8940.030.801F0.040.8940.020.801你说我记不住这张表怎么办没关系写项目时直接在配置里用枚举维护运行时按稳定度读取。但要知道 Briggs 公式的适用范围是 100 m 到 10 km 左右超过 10 km 扩散参数会失真尤其是 σ_z 在远距离下被高估导致浓度偏低。1.3 模型失效的几种现场别硬算高斯羽烟模型建立在“均匀稳定流场”的假设上现实里一旦出现下面四种情况算出来的等值面格网点基本就是数字垃圾一是静风或微风。u 接近 0 时公式直接除以 0此时湍流扩散主导高斯模型不再适用。行业里一般把 u 0.5 m/s 视为静风需要切换到其他模型。二是复杂地形。山谷、城市建筑物密集区会改变风向和湍流高斯羽烟假设平坦均匀地形在这种场景下要谨慎至少要把有效源高和初始扩散参数做修正。三是远距离传输。超过 10~20 km 后化学转化、干湿沉降会明显削减浓度纯高斯模型只适合中等距离的扩散估算。四是逆温层和混合层顶约束。公式里的垂直反射项只考虑地面如果混合层高度很低顶层也会反射需要加多次反射项。很多工程简化模型会完全忽略顶层反射导致地面浓度偏高。在开始写 Java 代码之前先把这个边界划好省得后面做完网格和等值面发现结果自洽但跟实测完全对不上那才是最尴尬的。2. Java 化建模从大气扩散公式到可运行的核心代码2.1 核心类怎么设计别把参数全塞进一个工具方法里我在第一个版本里偷懒把所有计算都写在一个 static 方法里参数一排排传进去结果后期加 GIS 坐标转换时恨不得重写。后来重构成了“参数对象 计算器 结果快照”的结构维护成本一下降下来了。基础配置类用于装载污染源的各项参数public class PlumeParams { private double q; // 排放速率 g/s private double u; // 烟囱出口风速 m/s private double h; // 有效源高 m private double yRef; // 源点经度十进制度 private double xRef; // 源点纬度十进制度 private double sourceY; // 源点的 GIS 坐标E/N private double sourceX; private Stability stability; // A-F 稳定度 private double zCalc; // 计算高度地面通常是 0 }计算器类只负责算单点浓度和格网浓度场不负责任何地理坐标转换。这样设计的好处是模型逻辑和 GIS 逻辑解耦后面不管你是接 WGS84 还是 CGCS2000计算器都不用动。2.2 单位换算和坐标基准最容易出错的一步所有计算都必须用米制。你从 GIS 里拿到的是经纬度或者投影坐标第一步要做的是确定一个平面坐标基准否则算出来的等都是乱的。我的做法是用源点作为局部坐标系原点东西方向为 y侧风向北方向为 x下风向。比如源点经纬度是 (E113.5, N34.2)用 UTM 投影转成投影坐标后源点坐标是 (460000, 3785000)。那么计算域里任意一点的相对偏移为deltaE pointE - sourceEdeltaN pointN - sourceN其中 deltaN 就是下风向距离 xdeltaE 就是侧风向距离 y。反过来生成等值面格网点时把每个格点的 (x, y) 相对偏移加上源点投影坐标再转回经纬度就能写出 GeoJSON。2.3 单点浓度计算实现核心方法我按公式直译没有做花哨的优化关键是让公式和代码一一对应方便排查public double calcConcentration(double downwindX, double crossY, double z) { if (this.u 0.5) { throw new IllegalStateException(风速过低高斯羽烟模型不适用); } double sigmaY getSigmaY(downwindX); double sigmaZ getSigmaZ(downwindX); double expY Math.exp(-crossY * crossY / (2.0 * sigmaY * sigmaY)); double expZ1 Math.exp(-(z - this.h) * (z - this.h) / (2.0 * sigmaZ * sigmaZ)); double expZ2 Math.exp(-(z this.h) * (z this.h) / (2.0 * sigmaZ * sigmaZ)); return this.q / (2.0 * Math.PI * this.u * sigmaY * sigmaZ) * expY * (expZ1 expZ2); }这里 sigmaY 和 sigmaZ 用 Briggs 公式x 用的是从污染源算起的下风向距离 downwindX。要特别注意这个距离不是网格点到源点的直线斜距而是沿风向的距离。如果你的场景里风向和正北方向有夹角那就得先把风轴旋转到正北方向再做坐标映射。最简单的做法是把所有点先绕源点旋转用风向角使得风轴和 x 轴重合算完浓度再逆旋转回去。3. 等值面格网点生成流程从零构建浓度场数据3.1 网格范围怎么取拍脑袋还是算一下网格范围和分辨率决定了计算量和等值面精度。我会先做一个快速判断地面浓度最大值通常出现在下风向 x 约等于有效源高 H 到 10 倍 H 的位置具体取决于稳定度。为了保证等值面覆盖完整下风向最大距离建议取 50 倍 H 到 100 倍 H侧风向覆盖到 5 倍 σ_y 以上。举例有效源高 50 m下风向取 5 km侧风向取 ±500 m网格间距 25 m那么网格规模是 200 × 40 8000 个点Java 单线程算起来很快。如果间距取 10 m规模变成 500 × 100 50000 个点也完全可接受。真正卡性能的是后面做等值线光滑和 GIS 导出。我给出一个自适应的估算逻辑int gridX (int) (downwindMax / gridStep) 1; int gridY (int) (crossSpread * 2 / gridStep) 1;这里 downwindMax 可以设为 max(100 * h, 5000)但不要超过混合层高度的量级否则模型失真。3.2 打网格并计算浓度场二维数组还是内存快照计算浓度场我通常维护一个 double[][] 矩阵行索引对应下风向位置列索引对应侧风向位置。每个格点都调用 calcConcentration 得到浓度值。有几点经验第一浓度跨度可能非常大烟囱近中心线可达几百 μg/m³边缘可能只有 0.001 μg/m³。等值面要想画出漂亮效果建议先对浓度取对数再用对数尺度生成等值线。这样即使浓差跨度很大制图时也不会出现边缘一片空白。第二对于浓度低于某个阈值比如 0.01 μg/m³的格点直接设成 NoData不参与等值线生成能避免后期异常等值线乱飞。第三步网格是规则矩阵但实际 GIS 坐标不是规整的因此要把矩阵和地理坐标映射关系保存下来。最简单的是一组双线性参数起始投影坐标、网格间距、旋转角任何格点都能还原成投影坐标。如果风向和正北方向存在夹角要额外保存旋转角我这里先按风轴与正北方向一致处理。3.3 从浓度格网提取等值线Marching Squares 的思路等值面格网点是散点浓度数据光是格网点没法直接当矢量面用必须提取等值线。GIS 里常见的等值线图本质上就是一系列浓度阈值对应的线。纯 Java 实现里最轻量的算法是 Marching Squares。它对网格中每个四边形单元取四个顶点浓度值和目标阈值 L 比较确定该单元与等值线的相交方式。一个四边形有 4 个顶点每个顶点两种状态大于 L 或小于 L组合起来共 16 种情况。等值线会穿过边界上那些一个顶点过阈值、另一个顶点未过阈值的边交点位置用线性插值得到。比如一条水平边两端浓度是 C1 和 C2那么插值因子 t (L - C1) / (C2 - C1)交点的横向坐标就在这条边的相对位置 t 处。把所有交点段连接起来就得到了该阈值下的等值线。Marching Squares 只生成线要素。要生成等值面填充面还需要判断每个网格单元是落在超过阈值的区域还是低于阈值的区域再把相邻单元合并成多边形。实际项目里如果只是想看图我一般生成等值线就够了如果要分析“超标面积”就得生成面。Java 侧我用 JTSJava Topology Suite做多边形合并和缓冲比手写几何运算稳定太多。4. GIS 集成输出把计算网格变成可分析的地图图层4.1 为什么不直接导出一张 PNG 图片很多初学 GISJava 的人会问我已经算出了浓度场直接画一张热力图不行吗答案是可以但你一旦把结果导出成 PNG就丢失了空间坐标信息后续叠加行政区划、计算超标人口、统计受影响面积都没法做。GIS 分析的核心对象是带坐标的矢量或者栅格数据不是图片。所以我们的最终输出应该是至少下面三种中的一种GeoJSON适合小型数据量前端 Leaflet、GitHub 直接能渲染Shapefile / GeoPackage适合进 ArcGIS、QGIS 做深度分析GeoTIFF适合做栅格分析比如重分类、叠加建模。从实践角度看Java 生成 GeoJSON 的回归成本最低因为 JTS 原生支持 GeoJSON 读写不需要额外转化库。如果项目需要发布成地图服务再考虑把 GeoJSON 导入 PostgreSQL PostGIS 后发布。4.2 栅格还是矢量CRS 和坐标基准决定数据命运栅格适合表达连续浓度场矢量适合表达等值线和等值面。两者可以共存我用网格算浓度场抽等值线输出矢量同时把原始浓度矩阵输出为 GeoTIFF 栅格这样 GIS 端既能做图又能做数值分析。需要极度重视坐标系。栅格 GeoTIFF 写入时必须带上投影信息不然 GIS 打开会乱套。Java 里我用 GeoTools 写 GeoTIFF坐标参考系统设为源数据的投影坐标比如 EPSG:32650UTM 50N。如果你的源数据是 WGS84 经纬度也可以直接把网格设置成等经纬度格网但要注意低纬度地区经纬度网格会导致距离变形模板计算不建议这么干。4.3 用 JTS 生成 GeoJSON 等值面核心代码片段先单独提一句JTS 的坐标顺序是 (x, y)对应 GIS 里的 (E, N)而经纬度是 (lon, lat)也就是 x经度、y纬度。很多哥们儿在写 GeoJSON 时把经纬度写反了生成的地物跑到非洲去了。下面这段代码展示了如何把一条等值线折线转成 GeoJSON FeatureGeometryFactory gf new GeometryFactory(); Coordinate[] coords contourLine.stream() .map(p - new Coordinate(p.x, p.y)) .toArray(Coordinate[]::new); LineString line gf.createLineString(coords); Feature feature new Feature(); feature.setGeometry(line); feature.setProperty(level, thresholdValue);当然实际项目里我更推荐直接用 JTS 的 GeoJSONWriter它能把 Geometry 直接序列化成标准 GeoJSON 字符串省去手写 JSON 结构的时间。5. 完整案例实测烟囱排放模拟从参数到图层5.1 场景设定让数据尽量贴近现实我在测试项目里设定一个简化场景某工厂烟囱高度 80 m有效源高按 100 m 处理SO₂ 排放速率 12 g/s出口处平均风速 3.5 m/s大气稳定度取 D 类中性风向为西北风也就是烟羽向东南方向扩散。源点经纬度取东经 113.5°、北纬 34.2°投影用 UTM 50N。网格设置为下风向 6000 m、侧风向 ±600 m、步长 30 mz 取 0地面浓度。这样生成 201 × 41 8241 个格点计算量非常小。5.2 代码流程计算、定阈值、抽等值线第一步用 UTM 转换把源点转成投影坐标。第二步遍历网格计算每个格点相对源点的东向偏移 deltaE 和北向偏移 deltaN。因为风向是西北风相当于风轴与正北方向有一个夹角我把 deltaE 和 deltaN 做旋转变换使得旋转后的坐标符合公式里的 x下风向、y侧风向。第三步计算每个格点地面浓度 c(x, y, 0)并把浓度矩阵保存起来。第四步设定一组等值线阈值比如 10、30、50、100、200 μg/m³用 Marching Squares 提取每个阈值下的等值线。第五步把等值线折线的点从相对坐标转回 UTM 投影坐标再转回经纬度写成 GeoJSON。5.3 结果合理性检查算完不能直接交差。我会做三个快速校验一是浓度峰值量级是否合理。经验公式是地面最大浓度约为 Q / (π * u * σ_y * σ_z * e) 附近的数量级如果你算出来峰值几百克每立方米那一定是单位错了。单位换算最容易出错的点在 Q 是 g/s但网格坐标是 m最后浓度应该落在一个合理的量级。二是等值线是否沿下风向拉长。侧风方向等值线应该比下风向方向窄。三是源点位置浓度是否过高。源点处由于 x 接近 0σ_y 和 σ_z 接近 0公式会出现异常高值所以源点附近几百米范围内要谨慎解读。用这次案例数据峰值大致出现在下风向约 500~800 m 处峰值浓度在 150 μg/m³ 上下整体形态符合 D 类稳定度的预期。等值线从源点开始呈锥形向下风向展开在 6 km 处已经衰减到 1 μg/m³ 以下。6. 项目落地中的常见坑与调优心得6.1 静风除零和源点发散静风处理我在代码里直接抛异常但真实项目里你不可能看天气不好就让系统崩掉。更合理的做法是当 u 低于 0.5 m/s 时切换到漫散型烟团模型或者直接提示“不满足高斯模型条件”。相比抛异常返回一个带有状态码的计算结果对象更友好前端还能源源不断渲染。源点发散问题尤其隐蔽。当 downwindX 很小且接近 0 时σ_y、σ_z 也趋向 0公式会出现巨大的浓度值。如果你对源点周边做了网格点计算很可能出现一个数值高出其他区域好几个数量级的点等值线图直接崩掉。我的经验是设一个最小下风向距离比如 10 m小于这个距离的格点浓度直接取源点近邻可用值或者标记为“源区不参与评估”。6.2 大网格性能优化从单线程到并行流如果网格扩大到 1000 × 1000 100 万个点单个线程循环计算可能耗时秒级。对 GIS 在线服务来说这个速度不能忍。Java 里最简单的优化就是并行流或者线程池因为每个格点计算之间没有数据依赖ConcurrentMapGridPoint, Double result gridPoints.parallelStream() .collect(Collectors.toConcurrentMap( p - p, p - model.calcConcentration(p.x, p.y, 0.0) ));还有一个很多人忽略的优化点σ_y 和 σ_z 每格都计算一次指数函数和幂函数成本高。如果固定风向下风向坐标相同的整列格点它们的 σ_y、σ_z 是一样的完全可以缓存。我把每行的 σ 值缓存到数组里性能提升非常明显尤其在大网格场景下。6.3 坐标基准混淆是最贵的坑我在开头说的 WGS84 和 CGCS2000 的问题不是段子。国内很多数据的底图是 CGCS2000但第三方接口提供的是 WGS84两者在高纬度地区差几十厘米在低纬度也有亚米级差异对等值面边界的位置会造成偏移。用 Java 做坐标转换时推荐用 GeoTools 自带的 CRS 转换能力不要自己写七参数公式。还要强调一个空间参考细节模型计算里的 x、y 应该是平面直角坐标系的米制单位而不是经纬度。一旦你把经纬度直接代进公式因为 1 度纬度和 1 度经度的米数不一样σ_y 和 σ_z 的单位全乱出来的等值面形状完全是错的。6.4 大结果数据集在 GIS 里的渲染问题如果你生成的是面状等值面而不是线状等值线相邻多边形之间容易出现裂缝或者重叠。原因在于 Marching Squares 生成的相邻单元公共边插值点位置可能会有 10⁻⁹ 级别的浮点误差。解决办法一是做节点捕捉二是用 JTS 的 Snap 工具把相邻多边形的节点合并到容差范围内。另外GeoJSON 里的坐标精度如果保留太多位小数文件体积会急剧膨胀。我的做法是输出投影坐标时保留 3 位小数毫米级精度输出经纬度时保留 6 位小数约 0.1 m 精度对等值面分析完全够用文件却小很多。最后再分享一个我个人的经验习惯每个高斯羽烟项目上线前我都会用同一组参数分别跑一次模型代码和一次专业大气预测软件的结果做对比。很多商业软件的核心扩散公式其实跟高斯模型同源如果两条等值线形状差太远多半是我的参数输入或者坐标旋转出了问题而不是模型算错了。有了这一步校准基本可以把等值面格网点数据放心交给 GIS 端做后续分析和展示了。