模仿dbscan聚类算法过滤清洗轨迹数据

10 篇文章 2 订阅
10 篇文章 5 订阅

 一般情况下,获取到的轨迹数据展示效果不会很好,需要进行预处理,如何删除大量没用数据呢?或者如何提取停留点,研究了三天,刚开始打算使用dbscan聚类来实现,但是发现效率非常低,2000多个点要好几分钟,遂放弃,自己仿照dbscan的思想写了一个算法。后面再补算法的过程,先来看下结果,删除了1172个点(圆点就是删除的点,或者说是停留点,总点数2200个左右),效果非常好。

算法过程:

 猜想:如果窗口设置小一点,循环的时候判断窗口内少于N个点,是不是也可以进行异常值的剔除。

 下面补上代码:

其中

    <script type="text/javascript" src="Scripts/map/004a20190418.js"></script> 是轨迹数据
    <script type="text/javascript" src="Scripts/map/map.js"></script> openlayers加载底图代码
    <script type="text/javascript" src="Scripts/map/TrackPlaying.js"></script>轨迹播放代码,这个例子也可以忽略,删除相应代码即可
    <script type="text/javascript" src="Scripts/kalman.js"></script>卡尔曼算法,还在研究,这里可以忽略

<!doctype html>
<html>
<meta charset="utf-8">
<head>
    <title>轨迹回放调用示例</title> 
    <meta name="viewport" content="width=device-width,initial-scale=1.0"> 
    <link rel="stylesheet" type="text/css" href="./Css/ol.css"> 
	<style>
	html{
    width: 100%;
    height: 100%;

	}
	body{
		width: 100%;
		height: 100%;
	}
	</style>
</head>
<body style="margin:0;" scroll="no">
    <div id="map" style="height: 100%;width:100%;position: relative;">
        <div style="width:190px;position: absolute;z-index: 11;right:10px;top:10px;background-color: #0D9BF2;opacity: 0.7;border-radius: 3px;padding: 10px;font-size: 12px">
            <table>
                <tr>
                    <td colspan="3">回放速度:<input type="range" id="speed" min="1" max="20" /></td>
                </tr>
                <tr>
                    <td><input type="button" value="开始回放" id="startButton" onclick="startAnimation();"style="cursor: pointer;"></td>
                    <td><input type="button" value="暂停" id="pauseButton" onclick="pauseAnimation();"style="cursor: pointer;"></td>
                </tr>
            </table>
        </div>
    </div>
     <script type="text/javascript" src="Scripts/map/004a20190418.js"></script>
    <script language="javascript" src="Scripts/ol-debug.js"></script>
    <script language="javascript" src="Scripts/jquery-1.4.4.js"></script>
    <script type="text/javascript" src="Scripts/map/map.js"></script>
    <script type="text/javascript" src="Scripts/map/TrackPlaying.js"></script>
	<script type="text/javascript" src="Scripts/sylvester.src.js"></script>
	<script type="text/javascript" src="Scripts/kalman.js"></script>
    <script type="text/javascript">
	var Collection_Frequency=10;//轨迹上报频率10s
	var EPS=24;//单位范围内的点数
	var RADIUS=0.0004;//正方形半径(半边长)
	
        var styles = [//建立样式数组
            new ol.style.Style({//0、轨迹路线样式
                stroke: new ol.style.Stroke({
                    color: 'azureblue',
                    width: 2
                })  
            }),
            new ol.style.Style({//轨迹点图标
                image: new ol.style.Icon({  
                    src: 'Images/taxi.png',                    
                    rotateWithView: false,  
                })
            })];
			
		var routeStyle =new ol.style.Style({
			stroke: new ol.style.Stroke({
				color: 'blue',
				width: 2
			})  
		});
		
		var kalmanStyle =new ol.style.Style({
			stroke: new ol.style.Stroke({
				color: 'red',
				width: 2
			})  
		});
		var pointStyle =new ol.style.Style({
			image: new ol.style.Circle({
				radius: 2,
				stroke: new ol.style.Stroke({
					color: '#60A090',
					width: 1
				}) 
			})  
		});
		var lineSource = new ol.source.Vector();//线条数据源
		var lineLayer = new ol.layer.Vector();//定义线条图层
		map.addLayer(lineLayer);
		
		var kalmanSource = new ol.source.Vector();//kalman线条数据源
		var kalmanLayer = new ol.layer.Vector();//定义kalman线条图层
		map.addLayer(kalmanLayer);
		
		var lineSourcePolygon = new ol.source.Vector();//采用的正方形图形
		var lineLayerPolygon = new ol.layer.Vector({
			type:'fanwei'
		});//
		lineLayerPolygon.setSource(lineSourcePolygon);
				map.addLayer(lineLayerPolygon);
		
		var routeStyle1 =new ol.style.Style({
			stroke: new ol.style.Stroke({
				color: 'red',
				width: 2
			})  
		});
		var deletePointSource = new ol.source.Vector();//剔除的点数据源
		var deletePointLayer = new ol.layer.Vector({
			type:'delePoint',
			//style:pointStyle
		});//剔除的点图层
		deletePointLayer.setSource(deletePointSource);
		map.addLayer(deletePointLayer);
		
		var vectorSourceHeatmap = new ol.source.Vector();
		// Heatmap            
		var vectorHeatmap = new ol.layer.Heatmap({
		  source: vectorSourceHeatmap,
		  blur: 15,
		  radius: 3,
		});
		map.addLayer(vectorHeatmap);
		var features=[];
		
		
        var routeCoords=[]
        var rowTrack=trackData.ROWDATA.ROW;
		var geometry = new ol.geom.LineString();
		var kalmangeometry = new ol.geom.LineString();
		
			
			
		var countThis=0;	//聚类数量
		//var countLast=0;
		var countDelete=0;//删除数量
		var countOutliers=0;//异常值数量
		var tempPre=[];
		var tempAft=[];
		var classArr=[];//每个聚类包含的点
		var classArrAll=[];//所有聚类
		//var indexArr=[];
		var countAll=rowTrack.length;
		//var tempCountY=0;

		//这里为了加载热力图
		for (var k = rowTrack.length-1; k >= 0 ; k--) {

			var feature=new ol.Feature(new ol.geom.Point([parseFloat(rowTrack[k].C_COORDINATEX),parseFloat(rowTrack[k].C_COORDINATEY)]));
			features.push(feature);	
		}
		//vectorSourceHeatmap.addFeatures(features); 

		for (var k = rowTrack.length-1; k >= 0 ; k--) {
			//想法1:这里可以先判断是否等于上一个点,等于的话直接删除(好像不合理)
			//想法2:允许出现某几次阈值距离外的畸变点,这样的畸变点也算入聚类中
			countLast=countThis;
			countThis=0;
			var geometryPolygon = new ol.geom.Polygon()
			var x1=parseFloat(rowTrack[k].C_COORDINATEX)-RADIUS;
			var y1=parseFloat(rowTrack[k].C_COORDINATEY)-RADIUS;
			var x2=parseFloat(rowTrack[k].C_COORDINATEX)+RADIUS;
			var y2=parseFloat(rowTrack[k].C_COORDINATEY)+RADIUS;
			geometryPolygon.setCoordinates([[[x1,y1],[x2,y1],[x2,y2],[x1,y2]]]);  
			for (var i =Math.min(k+EPS,rowTrack.length)-1 ; i >=Math.max(k-EPS,0) ; i--) {//只循环当前点前后EPS个点
				/*if(geometryPolygon.intersectsCoordinate([parseFloat(rowTrack[i].C_COORDINATEX),parseFloat(rowTrack[i].C_COORDINATEY)])){
					countThis++;
					temp.push(i);
				}*/
				var xi=parseFloat(rowTrack[i].C_COORDINATEX);
				var yi=parseFloat(rowTrack[i].C_COORDINATEY);

				if(xi<x2&&yi<y2&&xi>x1&&yi>y1){
					countThis++;
					tempPre.push(i);
					//tempCountX+=parseFloat(rowTrack[i].C_COORDINATEX);
					//tempCountY+=parseFloat(rowTrack[i].C_COORDINATEY);
				}
			}
			if(countThis>=EPS){
				countThis=kuoZhan(k,countThis,EPS,tempPre,[x1,y1,x2,y2]);
			}
			if(countThis>=EPS){
				var routeFeaturePolygon = new ol.Feature({
					geometry:geometryPolygon
				});
				routeFeaturePolygon.setStyle(routeStyle);
				//classArr=tempAft;
				//lineSourcePolygon.addFeature(routeFeaturePolygon);  
				if(tempPre[0]<=tempAft[tempAft.length-1]-1){//如果相邻两次有重合或者挨着
					//将两个数组合并
					classArr=hebing(classArr,tempPre);
				}else{
					classArrAll.push(classArr);
					classArr=[];
				}
				//console.log(tempPre+'-------'+k);
				for(var j=0;j<tempPre.length;j++){
					var featureDelete=new ol.Feature(new ol.geom.Point([parseFloat(rowTrack[tempPre[j]].C_COORDINATEX),parseFloat(rowTrack[tempPre[j]].C_COORDINATEY)]));
					featureDelete.set('geo',[x1,y1,x2,y2]);
					featureDelete.setStyle(pointStyle)
					deletePointSource.addFeature(featureDelete);
					console.log('====='+tempPre[j]);
					rowTrack.splice(tempPre[j],1);
					if(tempPre[j]<=k){
						k--;
					}
					countDelete++;
					
					//rowTrack[temp[j]].C_COORDINATEX=tempCountX/temp.length;
					//rowTrack[temp[j]].C_COORDINATEY=tempCountY/temp.length;
				}
			}
			
			//console.log(k+'-----------------------'+countThis);
			//countLast=0;
			//tempCountX=0;
			//tempCountY=0;
			tempAft=tempPre;
			tempPre=[];
			//routeCoords.push([(parseFloat(rowTrack[i-4].C_COORDINATEX)+parseFloat(rowTrack[i-3].C_COORDINATEX)+parseFloat(rowTrack[i-2].C_COORDINATEX)+parseFloat(rowTrack[i-1].C_COORDINATEX)+parseFloat(rowTrack[i].C_COORDINATEX))/5,(parseFloat(rowTrack[i-4].C_COORDINATEY)+parseFloat(rowTrack[i-3].C_COORDINATEY)+parseFloat(rowTrack[i-2].C_COORDINATEY)+parseFloat(rowTrack[i-1].C_COORDINATEY)+parseFloat(rowTrack[i].C_COORDINATEY))/5])
			
			
			 
			//geometry1.appendCoordinate([(parseFloat(rowTrack[i-4].C_COORDINATEX)+parseFloat(rowTrack[i-3].C_COORDINATEX)+parseFloat(rowTrack[i-2].C_COORDINATEX)+parseFloat(rowTrack[i-1].C_COORDINATEX)+parseFloat(rowTrack[i].C_COORDINATEX))/5,(parseFloat(rowTrack[i-4].C_COORDINATEY)+parseFloat(rowTrack[i-3].C_COORDINATEY)+parseFloat(rowTrack[i-2].C_COORDINATEY)+parseFloat(rowTrack[i-1].C_COORDINATEY)+parseFloat(rowTrack[i].C_COORDINATEY))/5]); 	
			
			
		}
		console.log(classArrAll);
		for(var j=0;j<rowTrack.length;j++){
				geometry.appendCoordinate([parseFloat(rowTrack[j].C_COORDINATEX),parseFloat(rowTrack[j].C_COORDINATEY)]);  
				routeCoords.push([parseFloat(rowTrack[j].C_COORDINATEX),parseFloat(rowTrack[j].C_COORDINATEY)])
			}
		console.log(countAll+'=====共删除====='+countDelete);
		//console.log(indexArr.join());
		var routeFeature = new ol.Feature({
			geometry:geometry
		});
		routeFeature.setStyle(styles[0]);
		lineSource.addFeature(routeFeature);  
		lineLayer.setSource(lineSource);
		
		
		
		
		
		//下面是卡尔曼滤波
		/* var x_0 = $V([111.47975954297902,33.30957616113584,0.6561951944928555,0.6561951944928555]);
		var P_0 = $M([[1,0,0,0],
						[0,1,0,0],
						[0,0,1,0],
						[0,0,0,1],
						]);
		var F_k=$M([[1.0,0.0,10.0,0.0],
					[0.0,1.0,0.0,10.0],
					[0.0,0.0,1.0,0.0],
					[0.0,0.0,0.0,1.0],
					]);
					
		var Q_k=$M([[0.5*(0.5*10*10)*(0.5*10*10),0.5*(0.5*10*10)*(0.5*10*10),0.5*(0.5*10*10)*(10*10),0.5*(0.5*10*10)*(10*10)],
					[0.5*(0.5*10*10)*(0.5*10*10),0.5*(0.5*10*10)*(0.5*10*10),0.5*(0.5*10*10)*(10*10),0.5*(0.5*10*10)*(10*10)],
					[0.5*(0.5*10*10)*(10*10),0.5*(0.5*10*10)*(10*10),0.5*10.0*10.0,0.5*10.0*10.0],
					[0.5*(0.5*10*10)*(10*10),0.5*(0.5*10*10)*(10*10),0.5*10.0*10.0,0.5*10.0*10.0],
					]);
		var KM = new KalmanModel(x_0,P_0,F_k,Q_k);

		var z_k = $V([111.4797595,33.309576,0.6561951944928555,0.6561951944928555]);
		var H_k = $M([[1.0,0.0,0.0,0.0],
					  [0.0,1.0,0.0,0.0]]);
		var R_k = $M([[0.01,0.0],[0.0,0.01]]);
		var KO = new KalmanObservation(z_k,H_k,R_k);

		var rowTrack=trackData.ROWDATA.ROW;

		var trackNew=[];
		for (var i=0;i<rowTrack.length;i++){
		  z_k = $V([parseFloat(rowTrack[i].C_COORDINATEX),parseFloat(rowTrack[i].C_COORDINATEY)]);
		  KO.z_k=z_k;
		  KM.update(KO);
		  console.log(KM.x_k.elements);
		  //trackNew.push([KM.x_k.elements[0],KM.x_k.elements[1]]);
		  kalmangeometry.appendCoordinate([KM.x_k.elements[0],KM.x_k.elements[1]]); 
		}
		var kalmanFeature = new ol.Feature({
			geometry:kalmangeometry
		});
		kalmanFeature.setStyle(kalmanStyle);
		kalmanSource.addFeature(kalmanFeature);  
		kalmanLayer.setSource(kalmanSource); */


		//下面是道格拉斯——普客抽希算法

		var dpResult=douglasPeucker(trackData.ROWDATA.ROW,2);
		console.log('-=-=-=-=-='+dpResult.length);
		var featuresDP=[];
		for (var k = dpResult.length-1; k >= 0 ; k--) {

			kalmangeometry.appendCoordinate([parseFloat(dpResult[k].C_COORDINATEX),parseFloat(dpResult[k].C_COORDINATEY)]);
				
		}
		var kalmanFeature = new ol.Feature({
			geometry:kalmangeometry
		});
		kalmanFeature.setStyle(kalmanStyle);
		kalmanSource.addFeature(kalmanFeature);  
		kalmanLayer.setSource(kalmanSource);
		//计算距离
		function calculationDistance(point1, point2){
			let lat1 = point1.C_COORDINATEY;
			let lat2 = point2.C_COORDINATEY;
			let lng1 = point1.C_COORDINATEX;
			let lng2 = point2.C_COORDINATEX;
			let radLat1 = lat1 * Math.PI / 180.0;
			let radLat2 = lat2 * Math.PI / 180.0;
			let a = radLat1 - radLat2;
			let b = (lng1 * Math.PI / 180.0) - (lng2 * Math.PI / 180.0);
			let s = 2 * Math.asin(Math.sqrt(Math.pow(Math.sin(a / 2), 2)
				+ Math.cos(radLat1) * Math.cos(radLat2) * Math.pow(Math.sin(b / 2), 2)));
			return s * 6370996.81;
		};
		//计算垂距
		function distToSegment(start, end, center){
			//下面用海伦公式计算面积
			let a = Math.abs(calculationDistance(start, end));
			let b = Math.abs(calculationDistance(start, center));
			let c = Math.abs(calculationDistance(end, center));
			let p = (a + b + c) / 2.0;
			let s = Math.sqrt(Math.abs(p * (p - a) * (p - b) * (p - c)));
			return s * 2.0 / a;
		};
		//递归方式压缩轨迹
		function compressLine (coordinate, result, start, end, dMax){
			if (start < end) {
				let maxDist = 0;
				let currentIndex = 0;
				let startPoint = coordinate[start];
				let endPoint = coordinate[end];
				for (let i = start + 1; i < end; i++) {
					let currentDist = distToSegment(startPoint, endPoint, coordinate[i]);
					if (currentDist > maxDist) {
						maxDist = currentDist;
						currentIndex = i;
					}
				}
				if (maxDist >= dMax) {
					//将当前点加入到过滤数组中
					result.push(coordinate[currentIndex]);
					//将原来的线段以当前点为中心拆成两段,分别进行递归处理
					compressLine(coordinate, result, start, currentIndex, dMax);
					compressLine(coordinate, result, currentIndex, end, dMax);
				}
			}
			return result;
		};
		/**
		 *
		 *@param coordinate 原始轨迹Array<{latitude,longitude}>
		 *@param dMax 允许最大距离误差
		 *@return douglasResult 抽稀后的轨迹
		 *
		 */
		function douglasPeucker(coordinate, dMax ){
			if (!coordinate || !(coordinate.length > 2)) {
				return null;
			}
			coordinate.forEach((item, index) => {
				item.key = index;
			});
			let result = compressLine(coordinate, [], 0, coordinate.length - 1, dMax);
			result.push(coordinate[0]);
			result.push(coordinate[coordinate.length - 1]);
			let resultLatLng = result.sort((a, b) => {
				if (a.key < b.key) {
					return -1;
				} else if (a.key > b.key)
					return 1;
				return 0;
			});
			resultLatLng.forEach((item) => {
				item.key = undefined;
			});
			return resultLatLng;
		};


		map.on('click',function( {pixel} ){
			var feature = map.forEachFeatureAtPixel(pixel,
			function(feature) {
			  return feature;
			});
			if(feature.get('geo')){
				console.log('------'+feature.get('geo'));
			}
		})
		
		
        var trackPlaying = new TrackPlaying(map, routeCoords, styles[0], styles[1]);//new一个轨迹播放对象,传入所需参数,map, routeCoords必须传入数据进去

        function startAnimation(){
            var speed = document.getElementById("speed").value;
            trackPlaying.ClearLayer();//每次点击先清空一下图层
            trackPlaying.speed =  speed;//设置速度
            trackPlaying.dealCoordsData();//过滤不合格经纬度数据
            trackPlaying.DrawLine();//画线
            trackPlaying.TrackMoving();//开始移动
            map.render();
            document.getElementById("startButton").value="停止回放";
            document.getElementById("startButton").onclick=stopAnimation
            document.getElementById("pauseButton").disabled=false;
        };

        function stopAnimation(){
            trackPlaying.StopMoving();//停止移动          
            document.getElementById("startButton").value="开始回放";
            document.getElementById("startButton").onclick=startAnimation
            document.getElementById("pauseButton").value="暂停";
            document.getElementById("pauseButton").onclick=pauseAnimation
            document.getElementById("pauseButton").disabled=true;
        }

        function pauseAnimation(){
            trackPlaying.PauseMoving();//暂停移动
            document.getElementById("pauseButton").value="继续";
            document.getElementById("pauseButton").onclick=continueAnimation
        }

        function continueAnimation(){
            trackPlaying.ContinueMoving();//继续移动
            document.getElementById("pauseButton").value="暂停";
            document.getElementById("pauseButton").onclick=pauseAnimation
        }
		
		
		function kuoZhan(k,count,eps,temp,polygonPoint){
			var eps2=eps+12;
			//console.log(polygonPoint);
			for (var i =Math.min(k+eps2,rowTrack.length)-1 ; i >=Math.min(k+eps2-12,rowTrack.length); i--) {//只循环向后扩展的30个点
				var xi=parseFloat(rowTrack[i].C_COORDINATEX);
				var yi=parseFloat(rowTrack[i].C_COORDINATEY);
				if(xi<polygonPoint[2]&&yi<polygonPoint[3]&&xi>polygonPoint[0]&&yi>polygonPoint[1]){
					count++;
					temp.push(i);
				}
			}
			for (var j =Math.max(k-eps2-12,0)-1 ; j >=Math.max(k-eps2,0); j--) {//只循环向前扩展的30个点
				var xj=parseFloat(rowTrack[j].C_COORDINATEX);
				var yj=parseFloat(rowTrack[j].C_COORDINATEY);
				if(xj<polygonPoint[2]&&yj<polygonPoint[3]&&xj>polygonPoint[0]&&yj>polygonPoint[1]){
					count++;
					temp.push(j);
				}
			}
			if(count>=eps2){
				kuoZhan(k,count,eps2,temp,polygonPoint);
			}
			return count;
		}
		
		//合并去重两个数组
		function hebing(list1,list2){
			var result=[];
			var i=0,j=0,k=0;//定义三个变量。i  j  k分别控制list1 list2  result三个数组的下标
			while(i<list1.length && j< list2.length){//两个数组都不为空的时候
				if(list1[i]<list2[j]){//若list1的元素小,加入result
					result.push(list1[i]);
					i++;//list1的下标后移
				}else if(list1[i]==list2[j]){//若二者相等 这一步就是去重复。
					result.push(list1[i]);//这里两个数组的元素把哪个加入result都行
					i++;//这里要注意的是  两个重复了,去重之后,两个数组的下标都要后移
					j++;
				}else {
					result.push(list2[j]);//若list2的元素小,加入result
					j++;
				}	
			}
			
			//下面是当其中一个数组元素全部添加到result了,另一个还没添加完  继续添加。
			while(i<list1.length){//
				result.push(list1[i]);
				i++;
			}
			while(j<list2.length){
				result.push(list2[j]);
				j++;
			}
			return result;
		}
		
    </script>
</body>
</html>

 

评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值