如何平滑拟合400米跑道？分段拟合问题

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



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


希望大佬给出努力方向。

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

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

你这个问题很有意思，我花了一些时间来试着寻优了一下，我这版程序的结果和人为匹配的结果存在偏差，感觉是分段函数的部分不太理想，暂时也没有想到有好的改进方法

1. function draft
   
2. X=importdata('GPS_X.txt');
   
3. Y=importdata('GPS_Y.txt');
   
4. 
   
5. L0=max(X)-min(X)-max(Y)+min(Y);
   
6. xc0=min(X)+(max(Y)-min(Y))/2;
   
7. yc0=min(Y)+(max(Y)-min(Y))/2;
   
8. x0=[L0 xc0 yc0];
   
9. 
   
10. 
    
11. options = optimoptions(@lsqnonlin,'Display','iter','OptimalityTolerance',1e-24,'StepTolerance',1e-50);
    
12. [x,e0]=lsqnonlin(@myfun,x0,[0 min(X) min(Y)],[max(X)-min(X) max(X) max(Y)],options,X,Y);
    
13. 
    
14. x0=x+[0.00005 -0.000025 0];
    
15. e2=myfun(x0,X,Y);
    
16. disp([e0 norm(e2)^2]);
    
17. 
    
18. R=mean(abs(Y(X>x(2) & X<x(2)+x(1))-x(3)));
    
19. xf=linspace(x(2)-R,x(2)+x(1)+R,1000);
    
20. [yf1,yf2]=fun1(x(1),x(2),x(3),R,xf);
    
21. x=x0;
    
22. R=mean(abs(Y(X>x(2) & X<x(2)+x(1))-x(3)));
    
23. xa=linspace(x(2)-R,x(2)+x(1)+R,1000);
    
24. [ya1,ya2]=fun1(x(1),x(2),x(3),R,xa);
    
25. 
    
26. figure(1)
    
27. plot(X,Y,[xf xf],[yf1 yf2],[xa xa],[ya1 ya2]);
    
28. end
    
29. 
    
30. function fun=myfun(x,X,Y)
    
31. L=x(1);
    
32. xc=x(2);
    
33. yc=x(3);
    
34. X1=X(X<=xc);
    
35. Y1=Y(X<=xc);
    
36. Y2=Y(X>xc & X<xc+L);
    
37. X3=X(X>=xc+L);
    
38. Y3=Y(X>=xc+L);
    
39. 
    
40. D1=var((Y1-yc).^2+(X1-xc).^2);
    
41. D3=var((Y3-yc).^2+(X3-xc-L).^2);
    
42. D2=var(abs(Y2-yc).^2);
    
43. fun=[D1;D2;D3];
    
44. end
    
45. 
    
46. function [yf1,yf2]=fun1(L,xc,yc,R,xf)
    
47. 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);
    
48. 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);
    
49. end

复制代码

你这个问题很有意思，我花了一些时间来试着寻优了一下，我这版程序的结果和人为匹配的结果存在偏差，感觉是分段函数的部分不太理想，暂时也没有想到有好的改进方法

1. function draft
   
2. X=importdata('GPS_X.txt');
   
3. Y=importdata('GPS_Y.txt');
   
4. 
   
5. L0=max(X)-min(X)-max(Y)+min(Y);
   
6. xc0=min(X)+(max(Y)-min(Y))/2;
   
7. yc0=min(Y)+(max(Y)-min(Y))/2;
   
8. x0=[L0 xc0 yc0];
   
9. 
   
10. 
    
11. options = optimoptions(@lsqnonlin,'Display','iter','OptimalityTolerance',1e-24,'StepTolerance',1e-50);
    
12. [x,e0]=lsqnonlin(@myfun,x0,[0 min(X) min(Y)],[max(X)-min(X) max(X) max(Y)],options,X,Y);
    
13. 
    
14. x0=x+[0.00005 -0.000025 0];
    
15. e2=myfun(x0,X,Y);
    
16. disp([e0 norm(e2)^2]);
    
17. 
    
18. R=mean(abs(Y(X>x(2) & X<x(2)+x(1))-x(3)));
    
19. xf=linspace(x(2)-R,x(2)+x(1)+R,1000);
    
20. [yf1,yf2]=fun1(x(1),x(2),x(3),R,xf);
    
21. x=x0;
    
22. R=mean(abs(Y(X>x(2) & X<x(2)+x(1))-x(3)));
    
23. xa=linspace(x(2)-R,x(2)+x(1)+R,1000);
    
24. [ya1,ya2]=fun1(x(1),x(2),x(3),R,xa);
    
25. 
    
26. figure(1)
    
27. plot(X,Y,[xf xf],[yf1 yf2],[xa xa],[ya1 ya2]);
    
28. end
    
29. 
    
30. function fun=myfun(x,X,Y)
    
31. L=x(1);
    
32. xc=x(2);
    
33. yc=x(3);
    
34. X1=X(X<=xc);
    
35. Y1=Y(X<=xc);
    
36. Y2=Y(X>xc & X<xc+L);
    
37. X3=X(X>=xc+L);
    
38. Y3=Y(X>=xc+L);
    
39. 
    
40. D1=var((Y1-yc).^2+(X1-xc).^2);
    
41. D3=var((Y3-yc).^2+(X3-xc-L).^2);
    
42. D2=var(abs(Y2-yc).^2);
    
43. fun=[D1;D2;D3];
    
44. end
    
45. 
    
46. function [yf1,yf2]=fun1(L,xc,yc,R,xf)
    
47. 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);
    
48. 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);
    
49. end

复制代码

coolchen302 发表于 2022-3-4 21:41
你这个问题很有意思，我花了一些时间来试着寻优了一下，我这版程序的结果和人为匹配的结果存在偏差，感觉是 ...

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

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

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

coolchen302 发表于 2022-3-5 08:35
那在前面加个线性变换即可
我说的问题还不是这个，等你看完后自己也想想吧 ...

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

coolchen302 发表于 2022-3-5 08:35
那在前面加个线性变换即可
我说的问题还不是这个，等你看完后自己也想想吧 ...

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

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

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

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

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

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;

使用直角坐标系数据。首先将数据点坐标减去 ，其中 是各数据点 坐标的平均值， 同理，这样一来，数据点就基本以坐标原点为中心分布了。设跑道直线段半长为 ，半圆段半径为 ，跑道中心点为 ，则跑道方程可表示为

这个跑道方程保证在分段点处连续。损失函数设计多样，我使用了 作为损失函数进行拟合，得到如下一组结果





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

944371430 发表于 2022-3-5 23:24
使用直角坐标系数据。首先将数据点坐标减去 ，其中  是各数据点  坐标的平均值， 同理，这样一来，数据点就 ...

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

顾世梁 发表于 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( ...

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

像这种在直角坐标系中封闭的凸曲线，可以变换到极坐标空间进行数据拟合。具体做法如下：
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');        % 绘制拟合效果图

复制代码