Java计算污染物扩散影响范围

高斯烟羽模型

高斯烟羽模型（Gaussian Plume Model）是一种广泛应用于大气污染扩散研究的数学模型。它主要用于估算大气中污染物从固定源（如烟囱）释放后，在不同风向和风速条件下污染物浓度的空间分布情况。该模型基于概率论和统计学原理，假设污染物在大气中的扩散遵循正态分布（高斯分布），因此得名“高斯烟羽模型”。

高斯烟羽模型的核心方程描述了污染物在下风向的浓度分布，主要受以下因素影响：

1. 风速：污染物扩散的主要动力来源。

2. 湍流扩散系数：描述垂直和水平方向上污染物随机扩散的程度。

3. 源强：单位时间内排放的污染物总量。

4. 源高：污染物排放口离地面的高度。

5. 大气稳定性：影响湍流扩散系数，与大气层结（稳定、不稳定或中性）有关。

6. 地形和障碍物：影响风场和湍流扩散。

高斯烟羽模型的特点是不需要显式地使用污染源的具体经纬度坐标来计算下风向某点的污染物浓度，因为模型的计算框架是基于相对坐标系的。在这个框架内，污染源通常被设定在坐标系的原点(0,0,0)，而所有的计算都是相对于这个原点进行的。这意味着无论污染源的真实地理位置在哪里，模型关注的是污染物扩散的方向和距离，即下风向的距离、横向和垂直方向的偏离。

高斯烟羽模型的基本假设之一是风向是已知的，并且在计算中将风向设定为x轴的方向。这样一来，模型只需要知道污染物释放点的高度（z方向的源高），以及污染物扩散点相对于释放点的下风向距离（x轴方向的距离）、横向距离（y轴方向的距离）和垂直距离（z轴方向的距离）。

当计算某个点的污染物浓度时，模型实际上是在计算相对于污染源原点的一个点的浓度，而不是关心这个点在地球表面上的绝对位置。因此，即使污染源和计算点的经纬度未知或未被直接使用，模型依然能够通过输入的相对坐标来预测污染物的浓度分布。

Maven依赖

    <dependency>
        <groupId>org.opengis</groupId>
        <artifactId>geoapi</artifactId>
        <version>2.3.1</version>
    </dependency>
    <dependency>
        <groupId>org.locationtech</groupId>
        <artifactId>jts-core</artifactId>
        <version>1.18.0</version>
    </dependency>
    <dependency>
        <groupId>org.apache.commons</groupId>
        <artifactId>commons-math3</artifactId>
        <version>3.6.1</version>
    </dependency>

环境条件类

@Data
@NoArgsConstructor
@AllArgsConstructor
public class EnvironmentalConditions {
    private double windSpeed;// 风速
    private double windDirection;// 风向
    private double humidity;// 湿度
    private double pressure;// 气压
    private double temperature;// 温度
    //...
}

污染源类

@Data
@NoArgsConstructor
@AllArgsConstructor
public class PollutantSource {
    private double emissionRate;// 排放速率
    private double height;// 排放口高度
    // ...
}

地形影响类

在计算污染物扩散影响范围时，地形数据对影响系数很重要，要考虑地面高度、地形遮挡、地形坡度，如果要更加精准的话还需要考虑山谷效应、城市热岛效应等因素。

public class TerrainEffect {
    
    /**
     * 计算地形对污染物扩散的影响，1.0表示无影响。
     */
    public static double getTerrainFactor(double x, double y, double z) {
        
        return 1.0;
    }
}

高斯扩散模型类

public class PollutionModel {
    
    public static double getElevation(double lat, double lon) {
        // 调用API ： https://api.opentopodata.org/v1/srtm90m?locations={lat},{lon}
        return 0.0;
    }
    
    /**
     * 将经纬度坐标转换为UTM坐标
     * @param targetCoord: 目标点的经纬度
     * @param sourceCoord: 污染源的经纬度
     */
    public static Coordinate convertToUTM(Coordinate targetCoord, Coordinate sourceCoord) {
        WKTReader reader = new WKTReader();
        
        // 定义WGS84地理坐标系
        CoordinateReferenceSystem crsWGS84 = CRS.decode("EPSG:4326");
        // 根据经纬度选择适当的UTM带
        CoordinateReferenceSystem crsUTM = CRS.decode(getUTMEPSGCode(sourceCoord.x, sourceCoord.y)); 
        
        // 创建几何对象
        GeometryFactory geometryFactory = new GeometryFactory(new PrecisionModel(), 0);
        Coordinate sourcePoint = reader.read("POINT (" + sourceCoord.y + " " + sourceCoord.x + ")").getCoordinate();
        Coordinate targetPoint = reader.read("POINT (" + targetCoord.y + " " + targetCoord.x + ")").getCoordinate();
        
        // 创建MathTransform对象
        MathTransform transform = CRS.findOperation(crsWGS84, crsUTM).getMathTransform();
        
        // 转换坐标
        Coordinate convertedSourcePoint = new Coordinate();
        Coordinate convertedTargetPoint = new Coordinate();
        transform.transform(sourcePoint, convertedSourcePoint);
        transform.transform(targetPoint, convertedTargetPoint);
        
        return convertedTargetPoint;
        
    }
    
    /**
     * 根据给定的污染源、环境条件和空间坐标计算污染物浓度。
     * 
     * @param source 污染源对象，包含排放率和排放口高度信息。
     * @param conditions 环境条件对象，包括风速、风向、温度、压力和湿度。
     * @param x 下风向的距离 米。
     * @param y 横向距离 米。
     * @param z 垂直距离 米。
     * @return 返回指定坐标处的污染物浓度。
     */
    public static double calculateConcentration(PollutantSource source, EnvironmentalConditions conditions, double x, double y, double z) {
        
        // 将坐标系转换成和风向一致
        double[] rotatedCoordinates = rotateCoordinates(x, y, conditions.getWindDirection());
        double rotatedX = rotatedCoordinates[0]; // 旋转后的下风向距离
        double rotatedY = rotatedCoordinates[1]; // 旋转后的横向距离
        
        // 计算横向和垂直方向上的标准差
        double sigmaY = calculateSigmaY(x, conditions);
        double sigmaZ = calculateSigmaZ(x, conditions);
        
        // 获取风速
        double windSpeed = conditions.getWindSpeed();
        
        // 高斯烟羽模型公式计算浓度
        double C = (source.getEmissionRate() / (2 * Math.PI * sigmaY * sigmaZ * windSpeed)) *
                   FastMath.exp(-rotatedY * rotatedY / (2 * sigmaY * sigmaY)) *
                   FastMath.exp(-(z - source.getHeight()) * (z - source.getHeight()) / (2 * sigmaZ * sigmaZ)) *
                   TerrainEffect.getTerrainFactor(rotatedX, rotatedY, z);
        
        return C;
    }
    
    /**
     * 根据经纬度获取UTM带编号
     */
    private static String getUTMEPSGCode(double longitude, double latitude) {
        int bandNumber = ((int)(longitude + 180) / 6) + 1;
        String hemispherePrefix = latitude >= 0 ? "326" : "327"; // 区分南半球或者北半球
        return "EPSG:" + hemispherePrefix + String.format("%02d", bandNumber); // 编码
    }
    
    /**
     * 旋转坐标系，使得污染物在新的坐标系中进行计算，考虑风向的影响。
     * 风向以度为单位，通过三角函数将原始坐标转换为旋转后的坐标。
     * 
     * @param x 原始坐标系中的x轴距离。
     * @param y 原始坐标系中的y轴距离。
     * @param windDirection 风向角度，单位为度。
     * @return 返回一个新的坐标数组，包含旋转后的x和y坐标。
     */
    private static double[] rotateCoordinates(double x, double y, double windDirection) {
        double windDirectionRadians = Math.toRadians(windDirection); // 旋转角度为弧度
        double rotatedX = x * Math.cos(windDirectionRadians) + y * Math.sin(windDirectionRadians);
        double rotatedY = -x * Math.sin(windDirectionRadians) + y * Math.cos(windDirectionRadians);
        return new double[]{rotatedX, rotatedY};
    }
    
     /**
     * 计算横向扩散系数sigmaY，基于下风向距离和环境条件。
     * 
     * @param x 下风向距离。
     * @param conditions 当前环境条件。
     * @return 返回横向扩散系数。
     */
    private static double calculateSigmaY(double x, EnvironmentalConditions conditions) {
        double baseSigmaY = 0.1 * x; // 初始横向扩散系数
        return adjustForConditions(baseSigmaY, conditions);
    }
    
    /**
     * 计算垂直扩散系数sigmaZ，基于下风向距离和环境条件。
     * 
     * @param x 下风向距离。
     * @param conditions 当前环境条件。
     * @return 返回垂直扩散系数。
     */
    private static double calculateSigmaZ(double x, EnvironmentalConditions conditions) {
        double baseSigmaZ = 0.05 * x; // 初始垂直扩散系数
        return adjustForConditions(baseSigmaZ, conditions);
    }
    
    /**
     * 根据环境条件调整扩散系数，考虑到温度、湿度和压力等因素。
     * 
     * @param sigma 初始扩散系数。
     * @param conditions 当前环境条件。
     * @return 返回调整后的扩散系数。
     */
    private static double adjustForConditions(double sigma, EnvironmentalConditions conditions) {
        double temperatureFactor = 1.0 + 0.01 * (conditions.getTemperature() - 20);
        double humidityFactor = 1.0 + 0.01 * (conditions.getHumidity() - 50);
        double pressureFactor = 1.0 + 0.001 * (conditions.getPressure() - 1013);
        
        return sigma * temperatureFactor * humidityFactor * pressureFactor;
    }
}

计算示例

public static void main(String[] args) {
        // 假设单位时间内排放量为100，排放口高度为10米
        PollutantSource source = new PollutantSource(100, 10);
        // 环境条件
        EnvironmentalConditions conditions = new EnvironmentalConditions(5, 0, 50, 1013, 25);
        
        // 污染源的经纬度
        Coordinate sourceCoord = new Coordinate(48.8566, 2.3522);
        // 目标点的经纬度
        Coordinate targetCoord = new Coordinate(48.8580, 2.3540);
        // 将经纬度坐标转换为UTM坐标
        Coordinate utmTargetCoord = PollutionModel.convertToUTM(sourceCoord, targetCoord);
       
        // UTM坐标转换为(x, y)
        double x = utmTargetCoord.x;
        double y = utmTargetCoord.y;
        // 获取目标点的海拔高度
        double z = PollutionModel.getElevation(targetCoord.y, targetCoord.x);
        
        double concentration = PollutionModel.calculateConcentration(source, conditions, x, y, z);
        
        System.out.println("污染物浓度为： (" + x + ", " + y + ", " + z + ") : " + concentration);
}