豪斯多夫距离（Hausdorff Distance）计算轨迹相似度

最近项目需求，要对运营车辆的运营报表，有个需求是按照运营线路的维度进行统计，需要将所有车辆的运营轨迹进行线路的匹配。平台中已经维护了所有的运营线路，所以只需要将车辆轨迹与运营线路进行相似度计算就可以了。

判断两条轨迹的相似性方法有很多:

基于点方法： EDR，LCSS，DTW等
基于形状的方法： Frechet， Hausdorff
基于分段的方法：One Way Distance, LIP distance
基于特定任务的方法：TRACLUS， Road Network，grid等
摘自如何判断两条轨迹（或曲线）的相似度？

这里用Hausdorff Distance来实现。

什么是豪斯多夫距离

Hausdorff距离是描述两组点集之间相似程度的一种量度，它是两个点集之间距离的一种定义形式。https://www.cnblogs.com/icmzn/p/8531719.html

豪斯多夫距离使用场景

根据资料，Hausdorff Distance较多用于物体轮廓的相似度判断，在计算机图像识别领域应用较多，如医疗领域肿瘤识别、军事领域打击物判断等。

具体代码(java + geotools)：
GeoTools中提供的HausdorffSimilarityMeasure.measure()方法，返回的是两条线路的相似度，返回值0-1，1表示完全一致。

POM仓库：

<repositories>
    <repository>
      <id>osgeo</id>
      <name>OSGeo Release Repository</name>
      <url>https://repo.osgeo.org/repository/release/</url>
      <snapshots><enabled>true</enabled></snapshots>
      <releases><enabled>true</enabled></releases>
    </repository>
  </repositories>

GeoTools依赖：

   <dependency>
      <groupId>org.geotools</groupId>
      <artifactId>gt-main</artifactId>
      <version>25.0</version>
    </dependency>

测试代码：

package com.yeetrack;

import com.alibaba.fastjson.JSONArray;
import com.alibaba.fastjson.JSONObject;
import lombok.extern.slf4j.Slf4j;
import org.geotools.geometry.jts.JTSFactoryFinder;
import org.geotools.referencing.GeodeticCalculator;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.junit.Test;
import org.locationtech.jts.algorithm.match.HausdorffSimilarityMeasure;
import org.locationtech.jts.geom.*;

/**
 * 文档：https://shengshifeiyang.gitee.io/geotools-learning/
 */
@Slf4j
public class GisUtil {

    /**
     * pi
     */
    public final static double pi = 3.1415926535897932384626;
    /**
     * a
     */
    public final static double a = 6378245.0;
    /**
     * ee
     */
    public final static double ee = 0.00669342162296594323;

    /**
     * @Description WGS84 to 火星坐标系 (GCJ-02)
     * @param lon 经度
     * @param lat 纬度
     * @return
     */
    public static double[] wgs84_To_Gcj02(double lon, double lat) {
        if (outOfChina(lat, lon)) {
            return null;
        }
        double dLat = transformLat(lon - 105.0, lat - 35.0);
        double dLon = transformLon(lon - 105.0, lat - 35.0);
        double radLat = lat / 180.0 * pi;
        double magic = Math.sin(radLat);
        magic = 1 - ee * magic * magic;
        double sqrtMagic = Math.sqrt(magic);
        dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi);
        dLon = (dLon * 180.0) / (a / sqrtMagic * Math.cos(radLat) * pi);
        double mgLat = lat + dLat;
        double mgLon = lon + dLon;
        return new double[] { mgLon, mgLat };
    }

    /**
     * @Description 火星坐标系 (GCJ-02) to WGS84
     * @param lon
     * @param lat
     * @return
     */
    public static double[] gcj02_To_Wgs84(double lon, double lat) {
        double[] gps = transform(lat, lon);
        double lontitude = lon * 2 - gps[1];
        double latitude = lat * 2 - gps[0];
        return new double[] { lontitude, latitude };
    }

    /**
     * @Description 火星坐标系 (GCJ-02) to 百度坐标系 (BD-09)
     * @param gg_lon
     * @param gg_lat
     * @return
     */
    public static double[] gcj02_To_Bd09(double gg_lon, double gg_lat) {
        double x = gg_lon, y = gg_lat;
        double z = Math.sqrt(x * x + y * y) + 0.00002 * Math.sin(y * pi);
        double theta = Math.atan2(y, x) + 0.000003 * Math.cos(x * pi);
        double bd_lon = z * Math.cos(theta) + 0.0065;
        double bd_lat = z * Math.sin(theta) + 0.006;
        return new double[] { bd_lon, bd_lat };
    }

    /**
     * @Description 百度坐标系 (BD-09) to 火星坐标系 (GCJ-02)
     * @param bd_lon
     * @param bd_lat
     * @return
     */
    public static double[] bd09_To_Gcj02(double bd_lon, double bd_lat) {
        double x = bd_lon - 0.0065, y = bd_lat - 0.006;
        double z = Math.sqrt(x * x + y * y) - 0.00002 * Math.sin(y * pi);
        double theta = Math.atan2(y, x) - 0.000003 * Math.cos(x * pi);
        double gg_lon = z * Math.cos(theta);
        double gg_lat = z * Math.sin(theta);
        return new double[] { gg_lon, gg_lat };
    }
    /**
     * @Description 百度坐标系 (BD-09) to WGS84
     * @param bd_lat
     * @param bd_lon
     * @return
     */
    public static double[] bd09_To_Wgs84(double bd_lon,double bd_lat) {
        double[] gcj02 = GisUtil.bd09_To_Gcj02(bd_lon, bd_lat);
        double[] map84 = GisUtil.gcj02_To_Wgs84(gcj02[0], gcj02[1]);
        return map84;

    }

    /**
     * @Description 判断是否在中国范围内
     * @param lat
     * @param lon
     * @return
     */
    public static boolean outOfChina(double lat, double lon) {
        if (lon < 72.004 || lon > 137.8347)
            return true;
        if (lat < 0.8293 || lat > 55.8271)
            return true;
        return false;
    }

    /**
     * @Description transform
     * @param lat
     * @param lon
     * @return
     */
    private static double[] transform(double lat, double lon) {
        if (outOfChina(lat, lon)) {
            return new double[] { lat, lon };
        }
        double dLat = transformLat(lon - 105.0, lat - 35.0);
        double dLon = transformLon(lon - 105.0, lat - 35.0);
        double radLat = lat / 180.0 * pi;
        double magic = Math.sin(radLat);
        magic = 1 - ee * magic * magic;
        double sqrtMagic = Math.sqrt(magic);
        dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi);
        dLon = (dLon * 180.0) / (a / sqrtMagic * Math.cos(radLat) * pi);
        double mgLat = lat + dLat;
        double mgLon = lon + dLon;
        return new double[] { mgLat, mgLon };
    }

    /**
     * @Description transformLat
     * @param x
     * @param y
     * @return
     */
    private static double transformLat(double x, double y) {
        double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * Math.sqrt(Math.abs(x));
        ret += (20.0 * Math.sin(6.0 * x * pi) + 20.0 * Math.sin(2.0 * x * pi)) * 2.0 / 3.0;
        ret += (20.0 * Math.sin(y * pi) + 40.0 * Math.sin(y / 3.0 * pi)) * 2.0 / 3.0;
        ret += (160.0 * Math.sin(y / 12.0 * pi) + 320 * Math.sin(y * pi / 30.0)) * 2.0 / 3.0;
        return ret;
    }

    /**
     * @Description transformLon
     * @param x
     * @param y
     * @return
     */
    public static double transformLon(double x, double y) {
        double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1* Math.sqrt(Math.abs(x));
        ret += (20.0 * Math.sin(6.0 * x * pi) + 20.0 * Math.sin(2.0 * x * pi)) * 2.0 / 3.0;
        ret += (20.0 * Math.sin(x * pi) + 40.0 * Math.sin(x / 3.0 * pi)) * 2.0 / 3.0;
        ret += (150.0 * Math.sin(x / 12.0 * pi) + 300.0 * Math.sin(x / 30.0* pi)) * 2.0 / 3.0;
        return ret;
    }

    /**
     * @Description WGS84 to 高斯投影(6度分带)
     * @param longitude 经度
     * @param latitude 纬度
     * @return double[] x y
     */
    public static double[] wgs84_To_Gauss6(double longitude, double latitude) {
        int ProjNo = 0;
        int ZoneWide; // //带宽
        double[] output = new double[2];
        double longitude1, latitude1, longitude0, X0, Y0, xval, yval;
        double a, f, e2, ee, NN, T, C, A, M, iPI;
        iPI = 0.0174532925199433; // //3.1415926535898/180.0;
        ZoneWide = 6; //6度带宽
        a = 6378137.0;
        f = 1.0 / 298.257223563; //WGS84坐标系参数
        //a = 6378245.0;f = 1.0 / 298.3; // 54年北京坐标系参数
        // //a=6378140.0; f=1/298.257; //80年西安坐标系参数
        ProjNo = (int) (longitude / ZoneWide);
        longitude0 = (double)(ProjNo * ZoneWide + ZoneWide / 2);
        longitude0 = longitude0 * iPI;
        longitude1 = longitude * iPI; // 经度转换为弧度
        latitude1 = latitude * iPI; // 纬度转换为弧度
        e2 = 2 * f - f * f;
        ee = e2 / (1.0 - e2);
        NN = a
                / Math.sqrt(1.0 - e2 * Math.sin(latitude1)
                * Math.sin(latitude1));
        T = Math.tan(latitude1) * Math.tan(latitude1);
        C = ee * Math.cos(latitude1) * Math.cos(latitude1);
        A = (longitude1 - longitude0) * Math.cos(latitude1);
        M = a
                * ((1 - e2 / 4 - 3 * e2 * e2 / 64 - 5 * e2 * e2 * e2 / 256)
                * latitude1
                - (3 * e2 / 8 + 3 * e2 * e2 / 32 + 45 * e2 * e2 * e2
                / 1024) * Math.sin(2 * latitude1)
                + (15 * e2 * e2 / 256 + 45 * e2 * e2 * e2 / 1024)
                * Math.sin(4 * latitude1) - (35 * e2 * e2 * e2 / 3072)
                * Math.sin(6 * latitude1));
        // 因为是以赤道为Y轴的，与我们南北为Y轴是相反的，所以xy与高斯投影的标准xy正好相反;
        xval = NN
                * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 14
                * C - 58 * ee)
                * A * A * A * A * A / 120);
        yval = M
                + NN
                * Math.tan(latitude1)
                * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24 + (61
                - 58 * T + T * T + 270 * C - 330 * ee)
                * A * A * A * A * A * A / 720);
        X0 = 1000000L * (ProjNo + 1) + 500000L;
        Y0 = 0;
        xval = xval + X0;
        yval = yval + Y0;
        output[0] = xval;
        output[1] = yval;
        return output;
    }

    /**
     * 计算两个点之间的距离（国测局02即火星坐标系）
     * @param x1
     * @param y1
     * @param x2
     * @param y2
     * @return
     */
    public static double distanceOfTwoGCJ02Points(double x1, double y1, double x2, double y2) {
        //火星坐标转wgs84坐标
        double[] wgsPointA = gcj02_To_Wgs84(x1, y1);
        double[] wgsPointB = gcj02_To_Wgs84(x2, y2);

        //gps84转高斯6度分带投影
        double[] gaussPointA = wgs84_To_Gauss6(wgsPointA[0], wgsPointA[1]);
        double[] gaussPointB = wgs84_To_Gauss6(wgsPointB[0], wgsPointB[1]);

        GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory();
        Point a = geometryFactory.createPoint(new Coordinate(gaussPointA[0], gaussPointA[1]));
        Point b = geometryFactory.createPoint(new Coordinate(gaussPointB[0], gaussPointB[1]));
        return a.distance(b);
    }

    /**
     * 计算两个点之间的距离（百度坐标系，即bd09）
     * @param x1
     * @param y1
     * @param x2
     * @param y2
     * @return
     */
    public static double distanceOfTwoBaiduPoints(double x1, double y1, double x2, double y2) {
        //bd09坐标转wgs84坐标
        double[] wgsPointA = bd09_To_Wgs84(x1, y1);
        double[] wgsPointB = bd09_To_Wgs84(x2, y2);

        //gps84转高斯6度分带投影
        double[] gaussPointA = wgs84_To_Gauss6(wgsPointA[0], wgsPointA[1]);
        double[] gaussPointB = wgs84_To_Gauss6(wgsPointB[0], wgsPointB[1]);

        GeodeticCalculator geodeticCalculator = new GeodeticCalculator(DefaultGeographicCRS.WGS84);
        geodeticCalculator.setStartingGeographicPoint(wgsPointA[0], wgsPointA[1]);
        geodeticCalculator.setDestinationGeographicPoint(wgsPointB[0], wgsPointB[1]);
        double distance = geodeticCalculator.getOrthodromicDistance();
        return distance;
    }

    /**
     * 两条线之间最短的距离
     * @param line1
     * @param line2
     * @return
     */
    public static double distanceOfTwoBaiduLines(LineString line1, LineString line2) {
        double minDis = Long.MAX_VALUE;
        long pointCount1 = line1.getCoordinates().length;
        long pointCount2 = line2.getCoordinates().length;
        for (int i=0; i<pointCount1; i++) {
            Coordinate point1 = line1.getCoordinateN(i);
            for (int j = 0; j<pointCount2; j++) {
                Coordinate point2 = line2.getCoordinateN(j);
                double temp = distanceOfTwoBaiduPoints(point1.getX(), point1.getY(), point2.getX(), point2.getY());
                if (temp < minDis) {
                    minDis = temp;
                }
            }
        }
        return minDis;
    }

    /**
     * 计算两点间的距离，百度坐标系，bd09
     * @param p1
     * @param p2
     * @return
     */
    public static double distanceOfTwoBaiduPoints(Point p1, Point p2) {
        return distanceOfTwoBaiduPoints(p1.getX(), p1.getY(), p2.getX(), p2.getY());
    }

    @Test
    public void testHausdorff() {
        GeometryFactory geometryFactory = JTSFactoryFinder.getGeometryFactory();
        String gps1 = "[{\"uploadTime\":\"2023-11-30 02:15:41\",\"longitude\":117.34200831126792,\"latitude\":39.12830255917371,\"speed\":9,\"accumulativeMileage\":43479},{\"uploadTime\":\"2023-11-30 02:15:51\",\"longitude\":117.34164337048267,\"latitude\":39.1282202856569,\"speed\":14,\"accumulativeMileage\":43479},。。。";
        JSONArray array1 = JSONArray.parseArray(gps1);
        Coordinate[] coordinates1 = new Coordinate[array1.size()];
        for (int i=0; i<array1.size(); i++) {
            JSONObject gps = array1.getJSONObject(i);
            double lon = gps.getDoubleValue("longitude");
            double lat = gps.getDoubleValue("latitude");
            coordinates1[i] = new Coordinate(lon, lat);
        }
        LineString lineString1 = geometryFactory.createLineString(coordinates1); //路线1

        String gps2 = "[{\"uploadTime\":\"2023-12-01 01:49:30\",\"longitude\":117.34165433596777,\"latitude\":39.12822050266268,\"speed\":18,\"accumulativeMileage\":43921},。。。。";
        JSONArray array3 = JSONArray.parseArray(gps3);
        Coordinate[] coordinates3 = new Coordinate[array3.size()];
        for (int i=0; i<array3.size(); i++) {
            JSONObject gps = array3.getJSONObject(i);
            double lon = gps.getDoubleValue("longitude");
            double lat = gps.getDoubleValue("latitude");
            coordinates3[i] = new Coordinate(lon, lat);
        }
        LineString lineString3 = geometryFactory.createLineString(coordinates3); //路线3，跟线路1完全没关系

        Coordinate[] coordinates4 = new Coordinate[array2.size()/2];
        for (int i=0; i<array2.size()/2; i++) {
            JSONObject gps = array2.getJSONObject(i);
            double lon = gps.getDoubleValue("longitude");
            double lat = gps.getDoubleValue("latitude");
            coordinates4[i] = new Coordinate(lon, lat);
        }
        LineString lineString4 = geometryFactory.createLineString(coordinates4); //线路4，取路线2的一段

        HausdorffSimilarityMeasure hausdorffSimilarityMeasure = new HausdorffSimilarityMeasure();
        double hausdorffDis = hausdorffSimilarityMeasure.measure(lineString1, lineString2);
        log.info("hausdorff message: {}", hausdorffDis);

        double hausdorffDis2 = hausdorffSimilarityMeasure.measure(lineString1, lineString3);
        log.info("hausdorff message: {}", hausdorffDis2);

        double hausdorffDis3 = hausdorffSimilarityMeasure.measure(lineString2, lineString4);

        log.info("hausdorff messure: {}", hausdorffDis3);
    }
}

版权声明

本站文章、图片、视频等（除转载外），均采用知识共享署名 4.0 国际许可协议（CC BY-NC-SA 4.0），转载请注明出处、非商业性使用、并且以相同协议共享。

© 空空博客，本文链接：https://www.yeetrack.com/?p=1561

Tags: java

豪斯多夫距离（Hausdorff Distance）计算轨迹相似度

什么是豪斯多夫距离

豪斯多夫距离使用场景

版权声明

You may also like...

发表回复取消回复

热门文章

近期文章

近期评论

什么是豪斯多夫距离

豪斯多夫距离使用场景

版权声明

You may also like...

php中文件的属性及方法

Java 使用Jedis操作Redis

springboot使用logback

发表回复 取消回复

热门文章

近期文章

近期评论

标签云

发表回复取消回复