【Matlab&Mathematica】对三维空间上的点进行椭圆拟合

OStack程序员社区-中国程序员成长平台 › 门户 › 编程›综合其他

原作者: [db:作者] 来自: [db:来源] 收藏邀请

问题是这样：比如有一个地心惯性系的轨道，然后从轨道上取了几个点，问能不能根据这几个点把轨道还原了？

当然，如果知道轨道这几个点的速度的情况下，根据轨道六根数也是能计算轨道的，不过真近点角是随时间变动的。

下面我会用数学的方法来解这个问题，基本思想是通过拟合空间上点的平面与椭球平面的交线将该轨道计算出来，算是一种思路吧。

首先需要有轨道数据，我们就从STK上获得，我使用默认参数生成了一个轨道，如下图：

输出j2000下的位置速度：

取其中5个点进行拟合：

可以先计算椭球，设椭球方程为x^2/a+y^2/b+z^2/c+d=0，然后求其最小二乘函数f(a,b,c,d) = sum((x^2/a+y^2/b+z^2/c+d)^2)。

通过单纯性法求该函数符合上面5个点的最小值时的a,b,c,d四个参数。

matlab里就是用fminsearch函数就行了。

代码如下：

clear all;
close all;
clc;

x=[4272.199656,8091.936548,9250.537919,7326.513290,4775.406342];
y=[-9122.847094,-8215.105235,-4377.145579,3801.632308,7889.177050];
z=[332.461913,8019.448934,13254.361230,16413.922920,15098.049929];

fun=@(t)sum((x.^2/t(1)+y.^2/t(2)+z.^2/t(3)+t(4)).^2);

[t fval]= fminsearch( fun , rand(4,1)) ; %单纯形法多元优化

t

这里的t就是要求的a,b,c,d四个参数。

我求的参数是

11831.9652077381
7748.53668377191
37674.4028071175
-14504.0838741735

得到椭球方程为：

x^2/11831.96520773810 + y^2/7748.536683771910 + z^2/37674.40280711745 - 14504.08387417353 = 0

Mathematica中画出来就是：

ContourPlot3D[
 x^2/11831.96520773810 + y^2/7748.536683771910 + 
   z^2/37674.40280711745 - 14504.08387417353 == 0, {x, -20000, 
  20000}, {y, -20000, 20000}, {z, -30000, 30000}]

椭球有了，下面我们求平面，使用ransac进行平面拟合。

我参考了这篇博客：https://blog.csdn.net/u010128736/article/details/53422070

不过他好像也参考了我过去写的ransac直线拟合：https://www.cnblogs.com/tiandsp/archive/2013/06/03/3115743.html，有趣 : )

代码如下：

clc;clear all;close all;

iter = 1000; 

x=[4272.199656,8091.936548,9250.537919,7326.513290,4775.406342];
y=[-9122.847094,-8215.105235,-4377.145579,3801.632308,7889.177050];
z=[332.461913,8019.448934,13254.361230,16413.922920,15098.049929];

data=[x;y;z];

%%% 绘制数据点
 figure;plot3(data(1,:),data(2,:),data(3,:),'o');hold on; % 显示数据点
 number = size(data,2); % 总点数
 bestParameter1=0; bestParameter2=0; bestParameter3=0; % 最佳匹配的参数
 sigma = 1;
 pretotal=0;     %符合拟合模型的数据的个数

for i=1:iter
 %%% 随机选择三个点
     idx = randperm(number,3); 
     sample = data(:,idx); 

     %%%拟合直线方程 z=ax+by+c
     plane = zeros(1,3);
     x = sample(1, :);
     y = sample(2, :);
     z = sample(3, :);

     a = ((z(1)-z(2))*(y(1)-y(3)) - (z(1)-z(3))*(y(1)-y(2)))/((x(1)-x(2))*(y(1)-y(3)) - (x(1)-x(3))*(y(1)-y(2)));
     b = ((z(1) - z(3)) - a * (x(1) - x(3)))/(y(1)-y(3));
     c = z(1) - a * x(1) - b * y(1);
     plane = [a b -1 c]

     mask=abs(plane*[data; ones(1,size(data,2))]);    %求每个数据到拟合平面的距离
     total=sum(mask<sigma);              %计算数据距离平面小于一定阈值的数据的个数

     if total>pretotal            %找到符合拟合平面数据最多的拟合平面
         pretotal=total;
         bestplane=plane;          %找到最好的拟合平面
    end  
 end
 %显示符合最佳拟合的数据
mask=abs(bestplane*[data; ones(1,size(data,2))])<sigma;    
hold on;
k = 1;
for i=1:length(mask)
    if mask(i)
        inliers(1,k) = data(1,i);
        inliers(2,k) = data(2,i);
        plot3(data(1,i),data(2,i),data(3,i),'r+');
        k = k+1;
    end
end

 %%% 绘制最佳匹配平面
 bestParameter1 = bestplane(1);
 bestParameter2 = bestplane(2);
 bestParameter3 = bestplane(4);
 xAxis = min(inliers(1,:)):max(inliers(1,:));
 yAxis = min(inliers(2,:)):max(inliers(2,:));
 [x,y] = meshgrid(-30000:1000:30000);
 z = bestParameter1 * x + bestParameter2 * y + bestParameter3;
 mesh(x, y, z);
 title(['bestPlane:  z =  ',num2str(bestParameter1),'x + ',num2str(bestParameter2),'y + ',num2str(bestParameter3)]);

拟合结果如下：

然后我们在Mathematica上画出平面与椭球：

代码如下：

data = {{4272.199656, -9122.847094, 
    332.461913}, {8091.936548, -8215.105235, 
    8019.448934}, {9250.537919, -4377.145579, 
    13254.361230}, {7326.513290, 3801.632308, 16413.922920},
   {4775.406342, 7889.177050 , 15098.049929}};

P3 = ListPlot3D[data, ColorFunction -> Function[{x, y, z}, Hue[x]]]

P2 = ContourPlot3D[
  1.8204 x + 0.81444 y - z - 20.3705 == 0, {x, -30000, 
   30000}, {y, -30000, 30000}, {z, -30000, 30000}]

P1 = ContourPlot3D[
  x^2/11831.96520773810 + y^2/7748.536683771910 + 
    z^2/37674.40280711745 - 14504.08387417353 == 0, {x, -20000, 
   20000}, {y, -20000, 20000}, {z, -30000, 30000}]

Show[{P1, P2, P3}]

结果如下：

感觉有那么点意思了。

下面我们求出两个面相交的参数方程，使用Mathematica计算还是比较方便的。

列出以下三个方程：

1.8204 x + 0.81444 y - z - 20.3705 = 0

x^2/11831.96520773810 + y^2/7748.536683771910 + z^2/37674.40280711745 - 14504.08387417353 = 0

z = 23375.89994468299 Sin[t]

这里23375.89994468299 通过sqrt(14504.08387417353*37674.40280711745)求出来的，z是在这个范围内变动的。

使用Mathematica解算如下：

Solve[{1.8204 x + 0.81444 y - z - 20.3705 == 0, 
  x^2/11831.96520773810 + y^2/7748.536683771910 + 
    z^2/37674.40280711745 - 14504.08387417353 == 0, 
  z == 23375.89994468299 Sin[t]}, {x, y, z}]

结果：

{{x -> 2.5336*10^-43 (3.90483*10^43 + 4.48093*10^46 Sin[t] - 7.40177*10^18 Sqrt[
       5.65524*10^54 - 8.37291*10^51 Sin[t] - 1.04594*10^55 Sin[t]^2]), 
  y -> 1.2668*10^-42 (2.28816*10^42 + 2.62575*10^45 Sin[t] + 3.30882*10^18 Sqrt[
       5.65524*10^54 - 8.37291*10^51 Sin[t] - 1.04594*10^55 Sin[t]^2]), 
  z -> 23375.9 Sin[t]},
 {x -> 2.5336*10^-43 (3.90483*10^43 + 4.48093*10^46 Sin[t] + 7.40177*10^18 Sqrt[
       5.65524*10^54 - 8.37291*10^51 Sin[t] - 1.04594*10^55 Sin[t]^2]), 
  y -> 1.2668*10^-42 (2.28816*10^42 + 2.62575*10^45 Sin[t] -  3.30882*10^18 Sqrt[
       5.65524*10^54 - 8.37291*10^51 Sin[t] - 1.04594*10^55 Sin[t]^2]),
  z -> 23375.9 Sin[t]}}

求到这里，这个参数方程已经求出来的，下面我画出来验证一下。

matlab代码如下：

clear all;
close all;
clc;

xx=[];yy=[];zz=[];

for t=0:0.001:2*pi
    x=2.5336*10^(-43) *(3.90483*10^43 + 4.48093*10^46*sin(t) -  7.40177*10^18*sqrt(5.65524*10^54 - 8.37291*10^51* sin(t) - 1.04594*10^55*sin(t).^2));
    y=1.2668*10^(-42) *(2.28816*10^42 + 2.62575*10^45*sin(t) + 3.30882*10^18*sqrt(5.65524*10^54 - 8.37291*10^51 *sin(t) - 1.04594*10^55* sin(t).^2));
    z=23375.9* sin(t);

    if(sqrt(5.65524*10^54 - 8.37291*10^51 *sin(t) - 1.04594*10^55* sin(t).^2)>0)
        xx=[xx x];
        yy=[yy y];
        zz=[zz z];
    end
end

for t=0:0.001:2*pi
    x=2.5336*10^(-43) *(3.90483*10^43 + 4.48093*10^46*sin(t) +  7.40177*10^18*sqrt(5.65524*10^54 - 8.37291*10^51* sin(t) - 1.04594*10^55*sin(t).^2));
    y=1.2668*10^(-42) *(2.28816*10^42 + 2.62575*10^45*sin(t) - 3.30882*10^18*sqrt(5.65524*10^54 - 8.37291*10^51 *sin(t) - 1.04594*10^55* sin(t).^2));
    z=23375.9* sin(t);
    if(sqrt(5.65524*10^54 - 8.37291*10^51 *sin(t) - 1.04594*10^55* sin(t).^2)>0)
        xx=[xx x];
        yy=[yy y];
        zz=[zz z];
    end
end

plot3(xx,yy,zz,'.')

x=[4272.199656,8091.936548,9250.537919,7326.513290,4775.406342];
y=[-9122.847094,-8215.105235,-4377.145579,3801.632308,7889.177050];
z=[332.461913,8019.448934,13254.361230,16413.922920,15098.049929];
hold on
plot3(x,y,z,'ro')

结果如下：

可以看出，点基本都在椭圆周围，效果不错，下面多用几组原始点对比看看：

呵呵，这就比较尴尬了，好像不怎么一样哦。

毕竟，这里只用了5个点，拟合点数多一些效果应该会好些吧 : )

关注公众号：只关于编程（ID：allcoding），每天推送MATLAB学习最常见的问题，每天进步一点点，业精于勤荒于嬉。

打开微信扫一扫哦！

鲜花

握手

雷人

路过

鸡蛋

该文章已有0人参与评论

请发表评论

全部评论

专题导读

More+

10-27 六六分期app的软件客服如何联系？(六六分期

11-06 可心卡盟:win10系统火狐flash插件崩溃怎么

11-06 亲亲特价:怎么删除回收站图标

11-06 济南大学虚拟社区:鲁大师节能降温的具体办

11-06 xlueops.exe:无线网络安装向导

11-06 女斗合众国:win7系统cf与主机连接不稳定怎

11-06 0xc000022-[cf烟雾头]cf怎么调烟雾头

11-06 qizideyouhuo:应用程序无法正常启动0xc0000

11-06 ipz-185:win7系统vcf文件怎么打开

11-06 傻哥蹦迪:win10系统s4怎么打开usb调试

11-06 八神浩树gtaste:回收站清空了怎么恢复

11-06 妖尾之黑色守护:win10系统电脑没有1440x900

11-06 校园至尊魔王小说:win7系统浏览网页时字体

11-06 女斗合众国:win10系统访问共享文件夹提示请

11-06 tokyo hot n0654:恢复win7系统默认字体一招

11-06 雨酷仙境:设置win7系统转移临时文件夹腾出

11-06 阿穆纳伊之杖:win7系统开始菜单在右边还原

11-06 tunespotting:win10系统火狐flash插件总是

11-06 甘尔葛分析师：计谋网站seo关键词暴涨有什

11-06 蔡贵霖: 计谋网站seo关键词暴涨有什么秘密

11-06 博益网首页:ao3网页版进入不了解决方法

11-06 漏斗子专栏: 网站数据分析小白易懂精华篇

11-06 见证双虹怎么做:win7系统开启telnet命令的

11-06 颾狐蝶蜋:系统资源不足无法完成请求的服务

11-06 国光中学校歌:提交网站到alexa查询详细步骤

11-06 西安有情天:静态网页和动态网页的区别

11-06 红木雅尚斋:外部链接构造对网站的好处

11-06 前官礼遇：防止域名劫持–增强域安全性的10

11-06 密传二转答案: 中文分词算法有哪些

11-06 金泉家园邮编:百度快照劫持的表现及应对方

Delphi下的COM编程技术 - 1发布时间：2022-07-18

Delphi实现Ftp客户端下载（支持断点续传，多线程传输）发布时间：2022-07-18

剪的笔顺,诠释剪的笔画,认识剪的部首

六六分期app的软件客服如何联系？(六六分期

2023-10-27

florent37/ViewAnimator: A fluent Android

2022-08-15

florent37/Shrine-MaterialDesign2: implem

2022-08-17

CVE-2020-36276

2022-09-23

SimpleSoftwareIO/simple-sms: Send and re

2022-08-13

阅读排行榜

1 六六分期app的软件客服如何联系？(六六分期

六六分期app的软件客服如何联系？不知道吗？加qq群【895510560】即可！标题：六六分期

阅读：19166|2023-10-27

2 可心卡盟:win10系统火狐flash插件崩溃怎么

今天小编告诉大家如何处理win10系统火狐flash插件总是崩溃的问题，可能很多用户都不知

阅读：9981|2022-11-06

3 亲亲特价:怎么删除回收站图标

今天小编告诉大家如何对win10系统删除桌面回收站图标进行设置，可能很多用户都不知道

阅读：8321|2022-11-06

4 济南大学虚拟社区:鲁大师节能降温的具体办

今天小编告诉大家如何对win10系统电脑设置节能降温的设置方法，想必大家都遇到过需要

阅读：8690|2022-11-06

5 xlueops.exe:无线网络安装向导

我们在使用xp系统的过程中,经常需要对xp系统无线网络安装向导设置进行设置，可能很多

阅读：8634|2022-11-06

6 女斗合众国:win7系统cf与主机连接不稳定怎

今天小编告诉大家如何处理win7系统玩cf老是与主机连接不稳定的问题，可能很多用户都不

阅读：9650|2022-11-06

7 0xc000022-[cf烟雾头]cf怎么调烟雾头

电脑对日常生活的重要性小编就不多说了，可是一旦碰到win7系统设置cf烟雾头的问题，很

阅读：8617|2022-11-06

8 qizideyouhuo:应用程序无法正常启动0xc0000

我们在日常使用电脑的时候，有的小伙伴们可能在打开应用的时候会遇见提示应用程序无法

阅读：7994|2022-11-06

9 ipz-185:win7系统vcf文件怎么打开

今天小编告诉大家如何对win7系统打开vcf文件进行设置，可能很多用户都不知道怎么对win

阅读：8648|2022-11-06

10 傻哥蹦迪:win10系统s4怎么打开usb调试

今天小编告诉大家如何对win10系统s4开启USB调试模式进行设置，可能很多用户都不知道怎

阅读：7531|2022-11-06

客服电话

电子邮件

【Matlab&Mathematica】对三维空间上的点进行椭圆拟合

请发表评论

全部评论

上一篇：

下一篇：

bluegill/katana:

OpenNBS/OpenNoteBlockStudio: An open-sou

juleswhite/mobilecloud-15: Programming C

JunichiT/JTMaterialSwitch: A Customizabl

delphi 数据导出进度条自己生成

剪的笔顺,诠释剪的笔画,认识剪的部首

六六分期app的软件客服如何联系？(六六分期

florent37/ViewAnimator: A fluent Android

florent37/Shrine-MaterialDesign2: implem

CVE-2020-36276

SimpleSoftwareIO/simple-sms: Send and re

关于我们

产品与服务

解决方案

139-2527-9053

客服电话

电子邮件

【Matlab&amp;Mathematica】对三维空间上的点进行椭圆拟合

请发表评论

全部评论

上一篇：

下一篇：

bluegill/katana:

OpenNBS/OpenNoteBlockStudio: An open-sou

juleswhite/mobilecloud-15: Programming C

JunichiT/JTMaterialSwitch: A Customizabl

delphi 数据导出 进度条自己生成

剪的笔顺,诠释剪的笔画,认识剪的部首

六六分期app的软件客服如何联系？(六六分期

florent37/ViewAnimator: A fluent Android

florent37/Shrine-MaterialDesign2: implem

CVE-2020-36276

SimpleSoftwareIO/simple-sms: Send and re

关于我们

产品与服务

解决方案

139-2527-9053

【Matlab&Mathematica】对三维空间上的点进行椭圆拟合

delphi 数据导出进度条自己生成