[已解决] 如何平滑拟合400米跑道?分段拟合问题

[复制链接]
Lantor_ 发表于 2022-3-3 20:28:25
本帖最后由 Lantor_ 于 2022-3-21 08:09 编辑

大家好,目前我遇到一个实际问题。
目的是获得一个400米操场的较为准确的坐标位置数据,已经使用gps采点实际采点,相当于有一组散点。然后,想用这些数据去拟合直线和半圆方程。

untitled.png

遇到的问题是,如果曲线拟合,这些数据不是函数;如果分段拟合的话,分段地方出现线段错位,图形不平滑。请问该如何解决?


希望大佬给出努力方向。

//2022.3.4 问题更新
今天在解决过程中意识到,因为散点是经纬度坐标(地球经度有180度,纬度之和有360度),所以两边部分可能在方程上不是一个圆形。我观察到直线也可能非水平,有一点倾斜度。

PS.文件里数据经过平移,保护隐私



GPS_Y.txt

23.02 KB, 下载次数: 23

GPS_X.txt

23.18 KB, 下载次数: 18

最佳答案


coolchen302 发表于 2022-3-4 21:41:56
你这个问题很有意思,我花了一些时间来试着寻优了一下,我这版程序的结果和人为匹配的结果存在偏差,感觉是分段函数的部分不太理想,暂时也没有想到有好的改进方法
  1. function draft
  2. X=importdata('GPS_X.txt');
  3. Y=importdata('GPS_Y.txt');

  4. L0=max(X)-min(X)-max(Y)+min(Y);
  5. xc0=min(X)+(max(Y)-min(Y))/2;
  6. yc0=min(Y)+(max(Y)-min(Y))/2;
  7. x0=[L0 xc0 yc0];


  8. options = optimoptions(@lsqnonlin,'Display','iter','OptimalityTolerance',1e-24,'StepTolerance',1e-50);
  9. [x,e0]=lsqnonlin(@myfun,x0,[0 min(X) min(Y)],[max(X)-min(X) max(X) max(Y)],options,X,Y);

  10. x0=x+[0.00005 -0.000025 0];
  11. e2=myfun(x0,X,Y);
  12. disp([e0 norm(e2)^2]);

  13. R=mean(abs(Y(X>x(2) & X<x(2)+x(1))-x(3)));
  14. xf=linspace(x(2)-R,x(2)+x(1)+R,1000);
  15. [yf1,yf2]=fun1(x(1),x(2),x(3),R,xf);
  16. x=x0;
  17. R=mean(abs(Y(X>x(2) & X<x(2)+x(1))-x(3)));
  18. xa=linspace(x(2)-R,x(2)+x(1)+R,1000);
  19. [ya1,ya2]=fun1(x(1),x(2),x(3),R,xa);

  20. figure(1)
  21. plot(X,Y,[xf xf],[yf1 yf2],[xa xa],[ya1 ya2]);
  22. end

  23. function fun=myfun(x,X,Y)
  24. L=x(1);
  25. xc=x(2);
  26. yc=x(3);
  27. X1=X(X<=xc);
  28. Y1=Y(X<=xc);
  29. Y2=Y(X>xc & X<xc+L);
  30. X3=X(X>=xc+L);
  31. Y3=Y(X>=xc+L);

  32. D1=var((Y1-yc).^2+(X1-xc).^2);
  33. D3=var((Y3-yc).^2+(X3-xc-L).^2);
  34. D2=var(abs(Y2-yc).^2);
  35. fun=[D1;D2;D3];
  36. end

  37. function [yf1,yf2]=fun1(L,xc,yc,R,xf)
  38. yf1=(sqrt(R^2-(xf-xc).^2)+yc).*(xf<=xc)+(yc+R).*(xf>xc & xf<xc+L)+(sqrt(R^2-(xf-xc-L).^2)+yc).*(xf>=xc+L);
  39. yf2=(-sqrt(R^2-(xf-xc).^2)+yc).*(xf<=xc)+(yc-R).*(xf>xc & xf<xc+L)+(-sqrt(R^2-(xf-xc-L).^2)+yc).*(xf>=xc+L);
  40. end
复制代码
回复此楼

12 条回复


coolchen302 发表于 2022-3-4 21:41:56
你这个问题很有意思,我花了一些时间来试着寻优了一下,我这版程序的结果和人为匹配的结果存在偏差,感觉是分段函数的部分不太理想,暂时也没有想到有好的改进方法
  1. function draft
  2. X=importdata('GPS_X.txt');
  3. Y=importdata('GPS_Y.txt');

  4. L0=max(X)-min(X)-max(Y)+min(Y);
  5. xc0=min(X)+(max(Y)-min(Y))/2;
  6. yc0=min(Y)+(max(Y)-min(Y))/2;
  7. x0=[L0 xc0 yc0];


  8. options = optimoptions(@lsqnonlin,'Display','iter','OptimalityTolerance',1e-24,'StepTolerance',1e-50);
  9. [x,e0]=lsqnonlin(@myfun,x0,[0 min(X) min(Y)],[max(X)-min(X) max(X) max(Y)],options,X,Y);

  10. x0=x+[0.00005 -0.000025 0];
  11. e2=myfun(x0,X,Y);
  12. disp([e0 norm(e2)^2]);

  13. R=mean(abs(Y(X>x(2) & X<x(2)+x(1))-x(3)));
  14. xf=linspace(x(2)-R,x(2)+x(1)+R,1000);
  15. [yf1,yf2]=fun1(x(1),x(2),x(3),R,xf);
  16. x=x0;
  17. R=mean(abs(Y(X>x(2) & X<x(2)+x(1))-x(3)));
  18. xa=linspace(x(2)-R,x(2)+x(1)+R,1000);
  19. [ya1,ya2]=fun1(x(1),x(2),x(3),R,xa);

  20. figure(1)
  21. plot(X,Y,[xf xf],[yf1 yf2],[xa xa],[ya1 ya2]);
  22. end

  23. function fun=myfun(x,X,Y)
  24. L=x(1);
  25. xc=x(2);
  26. yc=x(3);
  27. X1=X(X<=xc);
  28. Y1=Y(X<=xc);
  29. Y2=Y(X>xc & X<xc+L);
  30. X3=X(X>=xc+L);
  31. Y3=Y(X>=xc+L);

  32. D1=var((Y1-yc).^2+(X1-xc).^2);
  33. D3=var((Y3-yc).^2+(X3-xc-L).^2);
  34. D2=var(abs(Y2-yc).^2);
  35. fun=[D1;D2;D3];
  36. end

  37. function [yf1,yf2]=fun1(L,xc,yc,R,xf)
  38. yf1=(sqrt(R^2-(xf-xc).^2)+yc).*(xf<=xc)+(yc+R).*(xf>xc & xf<xc+L)+(sqrt(R^2-(xf-xc-L).^2)+yc).*(xf>=xc+L);
  39. yf2=(-sqrt(R^2-(xf-xc).^2)+yc).*(xf<=xc)+(yc-R).*(xf>xc & xf<xc+L)+(-sqrt(R^2-(xf-xc-L).^2)+yc).*(xf>=xc+L);
  40. end
复制代码
回复此楼

Lantor_ 发表于 2022-3-4 23:35:46
本帖最后由 Lantor_ 于 2022-3-4 23:47 编辑
coolchen302 发表于 2022-3-4 21:41
你这个问题很有意思,我花了一些时间来试着寻优了一下,我这版程序的结果和人为匹配的结果存在偏差,感觉是 ...

小弟先行谢过,万分感谢您的回答。  因为初学matlab,我还得认真研究您的代码。
针对您提到的衔接问题:今天我在尝试过程中意识到,因为GPS坐标散点是经纬度坐标(规定地球纬度之和有180度,经度之和有360度,1m距离在经纬度上对应不同的度数,猜想可能是2倍)。所以,当投射到直角坐标上石,两边部分可能在方程上不是一个圆形。  

coolchen302 发表于 2022-3-5 08:35:24
Lantor_ 发表于 2022-3-4 23:35
小弟先行谢过,万分感谢您的回答。  因为初学matlab,我还得认真研究您的代码。
针对您提到的衔接问题:今 ...

那在前面加个线性变换即可
我说的问题还不是这个,等你看完后自己也想想吧

TECHNOLOGYBOY 发表于 2022-3-5 12:37:58
本帖最后由 TECHNOLOGYBOY 于 2022-3-5 12:41 编辑
coolchen302 发表于 2022-3-5 08:35
那在前面加个线性变换即可
我说的问题还不是这个,等你看完后自己也想想吧 ...

刚刚我也试了您的程序,在直道和弯道衔接部分会早于实际gps踩点进行圆形拟合,我在想数据分成四个两组,弯道和直道,再进行拟合。弯道用cftool的高斯钟形二阶函数能达到0.95,直线的话用一次函数拟合 。

Lantor_ 发表于 2022-3-5 13:47:09
本帖最后由 Lantor_ 于 2022-3-5 13:49 编辑
coolchen302 发表于 2022-3-5 08:35
那在前面加个线性变换即可
我说的问题还不是这个,等你看完后自己也想想吧 ...

大佬好,感谢您的回复。我印证了一下坐标系的这个想法:就是把经纬度转化为直角坐标,再用你的源代码拟合如图,不知您怎么看,还有优化的空间么。




L}]EJ5ABSM8WLDG@H58AWNX.png
TS5OLU3Z{$~UUY`40RA{TCU.png

GPS_Y.txt

26.92 KB, 下载次数: 3

直角坐标

GPS_X.txt

26.92 KB, 下载次数: 3

直角坐标

untitled.fig

105.32 KB, 下载次数: 6


coolchen302 发表于 2022-3-5 16:25:37
Lantor_ 发表于 2022-3-5 13:47
大佬好,感谢您的回复。我印证了一下坐标系的这个想法:就是把经纬度转化为直角坐标,再用你的源代码拟合 ...

这个差不多了吧,肉眼看来已经比较接近了~~

coolchen302 发表于 2022-3-5 16:27:13
TECHNOLOGYBOY 发表于 2022-3-5 12:37
刚刚我也试了您的程序,在直道和弯道衔接部分会早于实际gps踩点进行圆形拟合,我在想数据分成四个两组,弯 ...

就是看怎么分组了,我的程序也是分成圆孤段和直线段来进行寻优的~

顾世梁 发表于 2022-3-5 21:14:10
本帖最后由 顾世梁 于 2022-3-5 21:30 编辑

SSy =       319.93
fx1=@(b,x)(x<=b(6)).*(b(1)+sqrt(b(2)+b(3)*x+b(4)*x.^2))+(x>b(6)).*(b(1)+b(5)).*(x<=b(7))+(x>b(7)).*(b(1)+sqrt(b(2)+b(3)*(x-b(8))+b(4)*(x-b(8)).^2));
fx2=@(b,x)(x<=b(6)).*(b(1)-sqrt(b(2)+b(3)*x+b(4)*x.^2))+(x>b(6)).*(b(1)-b(5)).*(x<=b(7))+(x>b(7)).*(b(1)-sqrt(b(2)+b(3)*(x-b(8))+b(4)*(x-b(8)).^2));
b =[5.070578128  -0.1793462142  1.555501295  -1.577245234  0.4559707586  0.4136961087  1.270876033  0.7356334918];
RSS = 0.167004092808
MSe = 7.7067e-05
R^2 =  0.99948note: 为便于显示和拟合,将原数据的尺度作了如下变换:
x=(x-84.08)*1000;
y=(y-91.46)*1000;
nh1737.png

944371430 发表于 2022-3-5 23:24:41
本帖最后由 944371430 于 2022-3-6 11:55 编辑

使用直角坐标系数据。首先将数据点坐标减去 ,其中 是各数据点 坐标的平均值, 同理,这样一来,数据点就基本以坐标原点为中心分布了。设跑道直线段半长为 ,半圆段半径为 ,跑道中心点为 ,则跑道方程可表示为
1.jpg
这个跑道方程保证在分段点处连续。损失函数设计多样,我使用了 作为损失函数进行拟合,得到如下一组结果

2.jpg
跑道.png


上述设计假设跑道相对于直角坐标系没有倾斜角。若跑道有倾斜角,则设计相对更复杂。







Lantor_ 发表于 2022-3-6 19:45:23
本帖最后由 Lantor_ 于 2022-3-6 19:50 编辑
944371430 发表于 2022-3-5 23:24
使用直角坐标系数据。首先将数据点坐标减去 ,其中  是各数据点  坐标的平均值, 同理,这样一来,数据点就 ...

大佬讲解的很清晰,图像拟合也很好,十分感谢!您和一楼的处理在图像上让我也意识到,转换坐标后,实际直线可能不是水平的。我拿去再研究学习,感谢感谢!!!

Lantor_ 发表于 2022-3-6 20:06:26
顾世梁 发表于 2022-3-5 21:14
SSy =       319.93
fx1=@(b,x)(xb(6)).*(b(1)+b(5)).*(xb(7)).*(b(1)+sqrt(b(2)+b(3)*(x-b(8))+b(4)*(x-b( ...

感谢大佬花时间解答,看图像的拟合程度相当接近,赞!小弟先拿去研究研究~

谢中华 发表于 2022-5-13 08:26:02
像这种在直角坐标系中封闭的凸曲线,可以变换到极坐标空间进行数据拟合。具体做法如下:
1. 将直角坐标数据 (xi, yi) 减去其均值,进行中心化变换;
2. 将中心化变换后的直角坐标数据转为极坐标数据 (ti, ri);
3. 对极坐标数据 (ti, ri) 进行数据拟合,可以选择合适的函数进行拟合,也可以做光滑样条插值拟合;
4. 将拟合后的极坐标数据转为直角坐标数据;
5. 最后对直角坐标数据进行逆中心化变换并绘图,得到光滑曲线。
相应的代码如下:
  1. x = readmatrix('GPS_X.txt');    % 读取数据
  2. y = readmatrix('GPS_Y.txt');    % 读取数据
  3. xm = mean(x);                   % 求均值
  4. ym = mean(y);                   % 求均值
  5. x1 = x-xm                       % 中心化
  6. y1 = y-ym;                      % 中心化
  7. [t,r] = cart2pol(x1,y1);        % 直角坐标转为极坐标
  8. ti = linspace(-pi,pi,2000);     % 定义新的角度向量
  9. ri = csaps(t,r,0.99,ti);         % 光滑样条插值拟合
  10. [xi,yi] = pol2cart(ti,ri);      % 极坐标转为直角坐标
  11. xi = xi+xm;                     % 逆中心化变换
  12. yi = yi+ym;                     % 逆中心化变换
  13. plot(x,y,'b',xi,yi,'r');        % 绘制拟合效果图
复制代码
拟合跑道.png

您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表