由勘探点生成三维地质层的算法实现-Jexonn's Blog

前言

这是我来公司后，部门领导给我安排的第一个任务，算法和业务相结合，还是很有挑战性的，值得记录一下。

这是产品经理给我画的核心算法的示意图：

大致意思是，勘探点数据中包含了每个勘探点的经纬度坐标和原始地表高程，而每个勘探点又包含了多个不同的岩土层数据，包括岩土层名称、岩土描述、岩土层深度、岩土层厚度，我需要根据前端传给我的勘探点 ID 列表，将其中隐含的地质层模型构建出来。由于地表是不平整的，构成每个地质层的上下表面点的高程不一定相同，所以构建出来的地质层是不规则的，需要将地质层的上下表面用三角网进行构建，侧面用四边形进行构建。

代码实现

实体类定义

三维坐标点 `Point3D`

package com.earthview.survey.domain.vo;

import lombok.Data;

/**
 * 3D点
 */
@Data
public class Point3D {
    public double x;
    public double y;
    public double z;

    public Point3D() {
    }

    public Point3D(double x, double y, double z) {
        this.x = x;
        this.y = y;
        this.z = z;
    }

}

地质点 `GeologicalPoint`

package com.earthview.survey.domain.vo;

import lombok.Data;

/**
 * 地质点
 */
@Data
public class GeologicalPoint extends Point3D {
    // 岩土名称（参与类别比较）
    String geotechnical_name;

    // 岩土描述（参与类别比较）
    String geotechnical_desc;

    // 地层某一层的上/下表面标记（用于后续算法）
    boolean isTop;

    // 勘探孔ID（用于后续算法，勘探孔ID和上/下表面标记均相同的两个地质点不可能为同一层地质层）
    Long holeId;

    public GeologicalPoint() {
        super();
    }

    public GeologicalPoint(double x, double y, double z, String geotechnical_name, String geotechnical_desc, String strata_number, boolean isTop, Long holeId) {
        super(x, y, z);
        this.geotechnical_name = geotechnical_name;
        this.geotechnical_desc = geotechnical_desc;
        this.isTop = isTop;
        this.holeId = holeId;
    }

    // 判断两个地质点是否属于同一类别
    public boolean isSameCategory(GeologicalPoint other) {
        return this.geotechnical_name.equals(other.geotechnical_name) &&
                this.geotechnical_desc.equals(other.geotechnical_desc);
    }

    public String toString() {
        return "GeologicalPoint{" +
                "x=" + x +
                ", y=" + y +
                ", z=" + z +
                ", geotechnical_name='" + geotechnical_name + '\'' +
                ", geotechnical_desc='" + geotechnical_desc + '\'' +
                ", isTop=" + isTop +
                ", holeId=" + holeId +
                '}';
    }
}

三角面片 `Triangle`

package com.earthview.survey.domain.vo;

import lombok.Data;

/**
 * 三角面片
 */
@Data
public  class Triangle {
    Point3D p1, p2, p3;

    public Triangle(Point3D p1, Point3D p2, Point3D p3) {
        this.p1 = p1;
        this.p2 = p2;
        this.p3 = p3;
    }
}

四边形面片 `Quad`

package com.earthview.survey.domain.vo;

import lombok.Data;

/**
 * 四边形面片
 */
@Data
public class Quad {
    Point3D p1, p2, p3, p4;

    public Quad(Point3D p1, Point3D p2, Point3D p3, Point3D p4) {
        this.p1 = p1;
        this.p2 = p2;
        this.p3 = p3;
        this.p4 = p4;
    }
}

地质层 `GeologicalLayer3D`

package com.earthview.survey.domain.vo;

import lombok.Data;

import java.util.List;

/**
 * 地质层
 */
@Data
public class GeologicalLayer3D {

    // 上表面地质点列表
    List<GeologicalPoint> topGeologicalPoints;

    // 下表面地质点列表
    List<GeologicalPoint> bottomGeologicalPoints;

    public String toString() {
        return "GeologicalLayer3D{" +
                "topGeologicalPoints=" + topGeologicalPoints +
                ", bottomGeologicalPoints=" + bottomGeologicalPoints +
                '}';
    }

    public String getGeotechnicalName() {
        return topGeologicalPoints.get(0).getGeotechnical_name();
    }

    public String getGeotechnicalDesc() {
        return topGeologicalPoints.get(0).getGeotechnical_desc();
    }
}

请求对象 `QueryByHoleIdListRequest`

package com.earthview.survey.request;  
  
import lombok.Data;  
  
import java.util.List;  
  
@Data  
public class QueryByHoleIdListRequest {  
    private List<Long> holeIds;  
}

响应对象 `GeologicalLayerVO`

package com.earthview.survey.response;

import com.earthview.survey.domain.vo.Point3D;
import com.earthview.survey.domain.vo.Quad;
import com.earthview.survey.domain.vo.Triangle;
import lombok.Data;

import java.util.List;

@Data
public class GeologicalLayerVO {
    private String geotechnicalName;
    private String geotechnicalDesc;
    private List<Triangle> topSurfaceTriangulation;
    private List<Triangle> bottomSurfaceTriangulation;
    private List<Quad> sideFaces;
    // 上表面中心点
    private Point3D topCenter;

    public String toString() {
        return "GeologicalLayerVO(geotechnicalName=" + this.getGeotechnicalName() + ", geotechnicalDesc=" + this.getGeotechnicalDesc() + ", topSurfaceTriangulation=" + this.getTopSurfaceTriangulation() + ", bottomSurfaceTriangulation=" + this.getBottomSurfaceTriangulation() + ", sideFaces=" + this.getSideFaces() + ")";
    }
}

算法接口 `buildGeologicalLayer`

@ApiOperation("生成三维地质层")  
@PostMapping("/buildGeologicalLayer")  
public AjaxResult buildGeologicalLayer(@RequestBody QueryByHoleIdListRequest request) {  
    List<Long> holeIds = request.getHoleIds();  
  
    // 如果 holeIds 长度小于3，则返回空列表  
    if (holeIds.size() < 3) {  
        return null;  
    }  
  
    // 地质点列表  
    List<GeologicalPoint> geologicalPoints = generateGeologicalLayerService.convertHoleIdsToGeologicalPoints(holeIds);  
  
    // 地质层列表  
    List<GeologicalLayer3D> geologicalLayer3DS = GeologicalLayerUtil.points2Layers(geologicalPoints);  
  
    // 构建VO对象  
    List<GeologicalLayerVO> geologicalLayerVOS = convertGeologicalLayer3DListToVOList(geologicalLayer3DS);  
    return AjaxResult.success(geologicalLayerVOS);  
  
}

服务类

GenerateGeologicalLayerServiceImpl：

package com.earthview.survey.service.impl;  
  
import com.baomidou.mybatisplus.core.conditions.query.QueryWrapper;  
import com.earthview.survey.domain.EvSurveyExplorationHoleGeology;  
import com.earthview.survey.domain.EvSurveyExplorationPointData;  
import com.earthview.survey.domain.vo.GeologicalPoint;  
import com.earthview.survey.service.GenerateGeologicalLayerService;  
import com.earthview.survey.service.IEvSurveyExplorationHoleGeologyService;  
import com.earthview.survey.service.IEvSurveyExplorationPointDataService;  
import org.springframework.beans.factory.annotation.Autowired;  
import org.springframework.stereotype.Service;  
  
import java.math.BigDecimal;  
import java.util.ArrayList;  
import java.util.List;  
  
@Service  
public class GenerateGeologicalLayerServiceImpl implements GenerateGeologicalLayerService {  
  
    @Autowired  
    private IEvSurveyExplorationPointDataService pointDataService;  
  
    @Autowired  
    private IEvSurveyExplorationHoleGeologyService holeGeologyService;  
  
    /***  
     * 根据给定的 holeId 进行生成该勘探点的所有地质点  
     * @param holeId 勘探点ID  
     * @return 地质点列表  
     */  
    @Override  
    public List<GeologicalPoint> convertHoleIdToGeologicalPoints(Long holeId) {  
        // 返回结果  
        List<GeologicalPoint> geologicalPoints = new ArrayList<>();  
  
        // 查询勘探点数据  
        EvSurveyExplorationPointData pointData = pointDataService.getById(holeId);  
        // 勘探点编号  
        String point_number = pointData.getPointNumber();  
        // 经度坐标  
        BigDecimal longitude = pointData.getLongitude();  
        // 纬度坐标  
        BigDecimal latitude = pointData.getLatitude();  
        // 勘探点原始地表高程  
        BigDecimal height = pointData.getHeight();  
  
  
        // 查询地质数据  
        QueryWrapper<EvSurveyExplorationHoleGeology> objectQueryWrapper = new QueryWrapper<>();  
        objectQueryWrapper.eq("hole_id", holeId);  
        List<EvSurveyExplorationHoleGeology> HoleGeologyList = holeGeologyService.list(objectQueryWrapper);  
        // 遍历 HoleGeologyList，并计算其高程  
        for (EvSurveyExplorationHoleGeology holeGeology : HoleGeologyList) {  
            // 地层编号 T-1            
            String strata_number = holeGeology.getStrataNumber();  
            // 岩土名称 填土  
            String geotechnical_name = holeGeology.getGeotechnicalName();  
            // 岩土描述  
            String geotechnical_desc = holeGeology.getGeotechnicalDesc();  
            // 岩土层深度(m)  
            BigDecimal rock_depth = holeGeology.getRockDepth();  
            // 岩土层厚度(m)  
            BigDecimal rock_thickness = holeGeology.getRockThickness();  
  
            // 根据勘探点原始地表高程和岩土层深度计算地质点高程  
            // 上表面高程(地表高程-（深度-厚度）)  
            BigDecimal z_top = height.subtract(rock_depth.subtract(rock_thickness));  
            // 下表面高程(地表高程-深度)  
            BigDecimal z_bottom = height.subtract(rock_depth);  
            // 创建地质点  
            // 上表面地质点  
            GeologicalPoint topGeologicalPoint = new GeologicalPoint(longitude.doubleValue(), latitude.doubleValue(), z_top.doubleValue(), geotechnical_name, geotechnical_desc, strata_number, true, holeId);  
            // 下表面地质点  
            GeologicalPoint bottomGeologicalPoint = new GeologicalPoint(longitude.doubleValue(), latitude.doubleValue(), z_bottom.doubleValue(), geotechnical_name, geotechnical_desc, strata_number, false, holeId);  
  
            geologicalPoints.add(topGeologicalPoint);  
            geologicalPoints.add(bottomGeologicalPoint);  
        }  
        return geologicalPoints;  
    }  
  
    /***  
     * 根据给定的多个 holeId 进行生成这些勘探点的所有地质点  
     * @param holeIds 勘探点ID列表  
     * @return 地质点列表  
     */  
    @Override  
    public List<GeologicalPoint> convertHoleIdsToGeologicalPoints(List<Long> holeIds) {  
        List<GeologicalPoint> geologicalPoints = new ArrayList<>();  
        for (Long holeId : holeIds) {  
            geologicalPoints.addAll(convertHoleIdToGeologicalPoints(holeId));  
        }  
        return geologicalPoints;  
    }  
  
}

EvSurveyExplorationHoleGeologyServiceImpl：

package com.earthview.survey.service.impl;  
  
import com.baomidou.mybatisplus.core.conditions.query.LambdaQueryWrapper;  
import com.baomidou.mybatisplus.core.conditions.query.QueryWrapper;  
import com.baomidou.mybatisplus.extension.plugins.pagination.Page;  
import com.baomidou.mybatisplus.extension.service.impl.ServiceImpl;  
import com.earthview.common.core.enums.ResponseEnum;  
import com.earthview.common.core.exception.ServiceException;  
import com.earthview.common.core.utils.survey.SnowFlakeUtil;  
import com.earthview.common.core.web.domain.AjaxResult;  
import com.earthview.common.core.web.page.PageDataResponse;  
import com.earthview.survey.domain.EvSurveyExplorationHoleGeology;  
import com.earthview.survey.domain.EvSurveyExplorationPointData;  
import com.earthview.survey.mapper.EvSurveyExplorationHoleGeologyMapper;  
import com.earthview.survey.mapper.EvSurveyExplorationPointDataMapper;  
import com.earthview.survey.request.IdListRequest;  
import com.earthview.survey.request.QueryByHoleIdListRequest;  
import com.earthview.survey.request.QueryHoleGeologyRequest;  
import com.earthview.survey.service.IEvSurveyExplorationHoleGeologyService;  
import org.springframework.beans.factory.annotation.Autowired;  
import org.springframework.stereotype.Service;  
  
import java.util.ArrayList;  
import java.util.HashMap;  
import java.util.List;  
import java.util.Map;  
import java.util.stream.Collectors;  
  
/**  
 * 岩土工程勘探孔地质成果表 服务实现类  
 */
@Service  
public class EvSurveyExplorationHoleGeologyServiceImpl extends ServiceImpl<EvSurveyExplorationHoleGeologyMapper, EvSurveyExplorationHoleGeology> implements IEvSurveyExplorationHoleGeologyService {  
  
    @Override  
    public List<EvSurveyExplorationHoleGeology> getHoleDataByIds(QueryByHoleIdListRequest request) {  
        LambdaQueryWrapper<EvSurveyExplorationHoleGeology> queryWrapper = new LambdaQueryWrapper<>();  
        queryWrapper.in(EvSurveyExplorationHoleGeology::getHoleId, request.getHoleIds().stream().map(Long::valueOf).collect(Collectors.toList()));  
        return this.baseMapper.selectList(queryWrapper);  
    }  
  
}

算法工具类

核心算法实现不公开，输入口令可见。

算法实现逻辑

数据准备

通过 GenerateGeologicalLayerServiceImpl 类的 convertHoleIdsToGeologicalPoints 方法，根据传入的勘探孔 ID 列表，提取出其中包含的所有地质点列表。
每个勘探孔会生成多个地质点，包括上表面和下表面的点，每个点包含三维坐标、岩土名称、描述等信息。

地质层划分

地质层划分是核心算法，主要在 GeologicalLayerUtil 类的 points2Layers 方法中实现。整个过程大致如下：

初始分组：
- 首先，将所有地质点按照地质类别（岩土名称和描述都相同的为一组）进行分组，确保相同类型的地质点被归为一组。
- 这一步为后续地质层的划分做了铺垫，因为如果每一类地质层只有一层，那么划分后的每一组地质点就可以直接构成一层地质层。但是实际上这一组相同类型的地质点可能在空间维度上构成了多层地质层，还需要进一步进行划分。
- 后续的算法都是针对这一步分好的每一组进行的。
勘探孔处理：
- 对于每一组地质点，再按勘探孔 ID 进行分组。
- 在每个勘探孔内，将地质点按 z 坐标降序排序。
- 然后，将相邻的两个点配对为上下表面点对（PointPair）。
由上下表面点对生成地质层列表：
- 使用 splitIntoSeparateLayers 方法将点对列表分割成不同的地质层。

splitIntoSeparateLayers 方法的详细实现逻辑如下：

初始排序：
- 将所有点对按照上表面点的 z 坐标进行降序排序。也就是说，后续的处理是从离地面最近的点对开始，向下遍历点对进行处理。
初始化：
- 创建一个用于存储最终地质层的列表 layers。
- 创建一个用于存储当前正在处理的地质层点对的列表 currentLayerPairs。
- 将排序后的第一个点对添加到 currentLayerPairs。
遍历点对：
- 从第二个点对开始，遍历所有剩余的点对。
- 对于每个点对：
  a. 获取当前点对 currentPair。
  b. 获取 currentLayerPairs 中的最后一个点对 lastPair。
收集构成同一地质层的所有点对：
- 使用 isInSameLayer 方法判断当前点对是否与上一个点对属于同一地质层。如果两个点对的上下表面范围内有重叠，则认为是同一层地质层。
- 如果属于同一层：
  - 将当前点对添加到 currentLayerPairs。
- 如果不属于同一层：
  - 即当前点对和currentLayerPairs中的点对不属于同一层地质层，说明已经收集好了一个地质层中包含的所有点对，使用 createLayer 方法创建一个新的地质层，并添加到 layers 列表中。
  - 清空 currentLayerPairs，并将当前点对作为新层的第一个点对。
处理最后一层：
- 遍历结束后，如果 currentLayerPairs 不为空，将这最后一层也添加到 layers 列表中。
返回结果：
- 返回包含所有分割好的地质层的 layers 列表。

这个方法的核心思想是通过比较相邻点对的空间关系来确定地质层的边界。它假设在同一地质层内，点对应该在空间上是连续的。当发现不连续时，就认为遇到了新的地质层。

三维模型构建

对每个地质层，使用 GeometryUtil 类的方法构建三维模型：

a. 使用 createDelaunayTriangulation 方法为上下表面创建 Delaunay 三角网格。

这里使用 JTS 库进行 Delaunay 三角剖分。Delaunay 三角剖分只处理二维坐标，而原始数据是三维点。通过 findClosestPoint 方法，找到与二维结果最接近的原始三维点，从而恢复Z坐标值。

为什么要用最 “接近” 的方法查找？因为在计算 Delaunay 三角剖分的过程中，可能会出现浮点数精度误差，导致三角剖分结果中的坐标与原始输入坐标有细微差异。使用 findClosestPoint 可以确保最终的三角形顶点使用的是原始输入的点，而不是经过计算后可能存在精度误差的点。

b. 使用 createSideFaces 方法创建地质层的侧面四边形。

这里使用 JTS 库创建上/下表面点集的凸包（包含所有点的最小凸多边形），从而获取地质层侧面一圈的边界点，用于创建地质层的侧面四边形。

数据转换

使用 GeologicalLayerUtil 类的 convertGeologicalLayer3DToVO 方法将 GeologicalLayer3D 对象转换为 GeologicalLayerVO 对象。
这个 VO 对象包含了地质层的名称、描述、上下表面的三角网格、侧面的四边形以及上表面的中心点（中心点方便前端渲染 “飞入”效果）。

结果返回

将处理后的 GeologicalLayerVO 列表作为结果返回。

前言

代码实现

实体类定义

三维坐标点 Point3D

地质点 GeologicalPoint

三角面片 Triangle

四边形面片 Quad

地质层 GeologicalLayer3D

请求对象 QueryByHoleIdListRequest

响应对象 GeologicalLayerVO

算法接口 buildGeologicalLayer

服务类