2015年全国大学生数学建模竞赛优秀论文

基于非线性曲线拟合的经纬度测量方法

摘要

本文首先基于天体物理学知识，构造出地球上某处直杆的影长与时间的函数关系式；然后运用非线性曲线拟合的方法，求解缺省参数，再根据直杆影长的变化规律，推算出测量点的地理位置及所处的日期。

在问题一中，本文以北京时间为参考时间，对地球上某一点处直杆影长的影响因素进行分析，发现其与直杆所处纬度、太阳直射点处纬度、所处时刻及经度等因素有关，结合地理知识构造出影长与影响因素的函数关系式。在各项参数均已给定的情况下，即可作出题目所要求的影长-时间变化曲线。

对于问题二，本文由附件1给定的时刻及其影长，运用非线性曲线拟合的方法，利用问题一中建立的关系式，将时间与影长作为已知参数，利用lsqcurvefit函数拟合求解经纬度参数。联系实际，筛选出可能的4个位置，并认为海南省白沙黎族自治县是最有可能的地点。

问题三与问题二基本相似，本文仍然在附件所得的数据基础上进行lsqcurvefit非线性曲线拟合，得到经度、纬度以及赤纬的可行解，根据所求赤纬，通过查表可以得到可能的日期。由附件2得到3个可能的地点与6个可能的日期，并认为其中新疆维吾尔自治区喀什地区巴楚县是最有可能的地点，5月24日或7月20日是最有可能的日期；由附件3同样得到3个可能的地点与6个可能的日期，认为湖北省十堰市郧西县与陕西省商洛市山阳县均是可能的地点，可能的日期为2月6日或11月6日前后。

对于问题四，首先用MATLAB进行图像处理并得到等时间间隔的图片，然后经过筛选得到21张图片。经滤镜处理后，由所得帧的图像得到影长与杆长的比例关系，进而得到不同时刻下的影长。在日期已知的情况下，问题四应用非线性拟合函数fit得到可行解，筛选后得到最可能地点为内蒙古自治区乌兰察布市丰镇市；若未给日期条件，在本题上一问的基础上，将太阳赤纬设为未知，利用fit函数求出可行解，经筛选得到最可能的地点为内蒙古自治区乌兰察布市，日期为6月6日或7月8日，与准确日期相差无几。

本文通过误差分析，证明本文所得结果具有很高的可信度。

关键词：非线性曲线拟合测量经纬度

一、问题重述

如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。

1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用你们建立的模型画出2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线。

2.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。将你们的模型应用于附件1的影子顶点坐标数据，给出若干个可能的地点。

3.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将你们的模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干个可能的地点与日期。

4.附件4为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度为2米。请建立确定视频拍摄地点的数学模型，并应用你们的模型给出若干个可能的拍摄地点。

如果拍摄日期未知，你能否根据视频确定出拍摄地点与日期

二、问题分析

1.对于问题一，依据地理学知识，直杆的影长与直杆长度以及太阳高度角有关，而太阳高度角又与赤纬、直杆所处纬度、测量时刻等参数有关。考虑到附件所给时间为北京时间，还应考虑经度带来的时差影响。基于以上考虑，构造出直杆的影长与赤纬、所处经度和纬度以及时刻之间的函数关系式，进而将题目所给参数代入关系式中，即可作出在给定时间天安门广场上给定直杆的影长随时间的变化曲线。

2.对于问题二，若想得到可能的地点，关键就在于得到直杆所处位置的经纬度坐标。由所给附件可以得到影长随时间的变化关系，因此可以代入问题一种建立的模型，通过最小二乘法迭代的方法估计出参数的可能取值，即直杆可能位置的经纬坐标。

3.对于问题三，仅是在问题二的基础上增多了日期这一未知量。这一未知量可以通过赤纬与日期对应表查询，因此同第二问一样，可以由附件得到影长随时间的变化关系，利用最小二乘法迭代估计出经纬度、赤纬三个参数可能取值，从而给出可能的地点和日期。

4.对于问题四，通过MATLAB导入视频，并等时间间距地21张截图。在经过滤镜处理后，可以由截图得到影长与杆长的比例关系，从而得到不同时间下的影长。在日期已知的情况下，问题可以利用问题二的模型；若日期未知，则可以利用问题三的模型，得到所求参量。

三、模型假设

1.地球是半径为R的均匀圆球，即忽略其表面地形及地球微小的椭球率带来的影响；

2.忽略太阳光在穿过地球大气层时产生的折射；

3.照射到地球的太阳光为平行光；

4.黄赤交角为一固定值。

四、 符号约定

太阳直射点纬度 测量点的纬度 测量点的经度 测量点的时角 太阳直射方向向量

地心指向测量点的方向向量

h 太阳高度 直杆长 影长 测量点地方时

测量点地方时对应的北京时间

五、模型建立与求解

（一）问题一

根据假设1，以如图所示的方式建立空间直角坐标系[1]。其中，以地心为坐标原点，垂直于赤道平面指向正北的轴为z 轴，沿赤道平面指向太阳方向的轴为x 轴，与二者成左手关系的轴为y 轴。假定此时太阳直射点为地球上一点A ，设其纬度为0?（北纬为正，南纬为负，以下涉及纬度符号均同理）。显然，由于x 轴指向太阳方向，因此直射点A 位于oxz 平面内。设此时刻太阳直射方向为1e r

。设测量点纬度为?，其时角[2]为a ?，其与地心构成的向量为2e r

。

图一模型推导示意图

根据以上设定，有

可求得1e r 与2e r

之间的空间θ的余弦值为

12

0012cos cos cos cos sin sin ||||

a e e e e θ??????==+r r

r r （1）

根据太阳高度角h 的定义，有2

h π

θ=

-，因此可得

00sinh cos cos cos sin sin a ?????=+（2）

由图二所示的示意图可得，影长、直杆长之间存在关系

tanh

L L =

（3） 图二影长、杆长与太阳高度的关系

在上述推导中，时角(12)12

a t π

?=-?

，其中t 为测量点地方时。考虑到附录中所给时间均为北

京时间，因此此处统一将地方时转化为北京时间表示。设测量点经度为l ，为方便起见，不妨假设测量点在北京所在的东八区的中间经线（120°E ）以东，因此二者之间经度差为

由于地球每自转24h 便转过360°，因此近似可认为经度上每度引起的时差为4min ，由前可知测量点在120°E 以东，因此测量点地方时与东八区区时有如下关系 因此

0(120)121512

a l t π???-=+-?????o （4）

其中0t 为测量点地方时t 对应的北京时间。

联立式（2）、（3）、（4），可得地球上某一点处直杆的影长与北京时间的函数关系式

0000

(120)

cos cos cos 12sin sin 1512L L l t π????=??-+-?+?????

???o

（5） 查表[3]可知2015年10月22日太阳直射点010.833?=-o

，并将已知条件03L =m ，

39.9072?=o ，116.3914l E =o ，使用MATLAB 即可作出题目给定时间内天安门广场上3m 长的直

杆的影长变化曲线，如图三所示。

图三问题一所求曲线

（二）问题二

由附件可以得到测量点直杆的影长与北京时间的对应数据，可以在问题一所得到的模型的基础上，将待求的经、纬度坐标以及杆长看作为模型的参数，调用MATLAB 中lsqcurvefit [4]函数，通过迭代的方法实现最小二乘拟合，得到参数的估计值。

Lsqcurvefit 函数是进行最小二乘非线性曲线拟合的函数，设定初始迭代起点进行拟合求最优解。 为了更快得到最优解，将迭代起点选择为期望值附近（但并未限定所求解的范围），即杆长期望在0至3米之间，经度位于-180°至180°之间，纬度的正弦值位于-1到1之间。

求解完成后，基于以下原则，在MATLAB 返回的可行解当中进行了初步的筛选： 1. 所有参数均应为纯实数，即虚部应等于零； 2. 出于实际考虑，要求杆长均大于；

3.参数完全重复的或相当接近的，取其中拟合程度最佳的一组参数。

由附件1得到的数据见下表（时间已化为十进制）：

北京时间（时）影长（m）

15

表一由附件1计算得到的北京时间与对应影长具体程序见附录，在代入由附件1得到的数据之后，直接得到的可行解如下：杆长（m）纬度经度出现次数

30

6 ………

1

1

34

………

表二附件1数据第一次迭代后的结果

由表格可以看到某些解出现的次数不止一次，这是因为不同点起点得到相同的结果，说明该结果可能为最优解。

由表格知（98°N，°E）在可行解中出现的频率极高，说明其可信度相对于其他解要更高；而在（°~°N，°~°E）这一范围中，尽管每组解出现的次数不高，但是可以发现这些数据分布相当密集，因此可以看做是同一个地点；（°~°S，~°E）范围之间的解同理，但是由于其出现次数较少，显然其可信度不及另两组解。

为了增强以上结果的可信度，现改变迭代起点进行第二次拟合。在这一次拟合中，将迭代起点设置偏离期望值。得到的结果如下表所示：

杆长（m）纬度经度出现次数

1

1

表三附件1数据第二次迭代后的结果

（°N，°E）这一组解再次出现，证明了这一组解得可信度确实相当高。将以上得到的可能解汇总，并确定具体的地点如下：

编号

杆长

（m）

纬度经度具体地点

1 海南省白沙黎族自治县

2 云南省德宏傣族景颇族自治州盈江县

3 印度尼西亚SinarGunungTebatKaraiKepahiangRegency

4 印度尼西亚AsiaBaruKuripanBaritoKualaRegency

表四问题二可能的地点

分别对得到的四组解进行误差可视化和再次拟合，图像如下（从上到下依次为编号1、2、3、4的解）：

图四问题二所有可能解的拟合及误差

误差图显示拟合值与真实值之间的误差均为10-3级别，这再一次证明了解的可信度。这是基于本模型得到的四个可能的地点。考虑到迭代过程中出现次数悬殊较大，有理由相信海南省白沙黎族自治县（°N，°E）是最有可能的地点，而其他三个地点可能性相对较小。

（三）问题三

第三问在第二问的基础上又增添了日期这个未知参数。由于日期与赤纬存在一定的对应关系，同第二问思路相同，我们将赤纬作为另一个进lsqcurvefit拟合，此时初步筛选原则除去问题二中提出的3条外，再增添一条：赤纬的范围需位于°S~°N之间。

1.对附件2的分析

由附件2得到的数据见下表（时间已化为十进制）：

北京时间（时）影长（m）

表五由附件2计算得到的北京时间与对应影长

将附件2提取的数据代入模型，直接得到的可行解如下：

杆长（m）纬度经度赤纬出现次数

10

1

1

表五附件2数据第一次迭代后的结果

同样，为了增强解得可信度，使迭代起点偏离期望住进行了第二次拟合，直接得到的可行解如下：

杆长（m）纬度经度赤纬出现次数

1

表六附件2数据第二次迭代后的结果

由此可见，第二次得到的这组解可信度很高。将表四中三组解依次编号为1、2、3，以赤纬进行查表确定日期，最终得到下列可能的时间和地点：

编号纬度经度具体地点日期

1 新疆维吾尔自治区喀什地区巴楚县5月24日或7月20日

2 缅甸海2月18日或10月25日

3 缅甸海3月5日或10月9日

表七问题三之附件2可能的地点和日期

分别进行误差可视化，得到的图像如下（从上到下依次为编号1、2、3的解）：

图五附件2所有可能解的拟合及误差

从误差图及拟合图都可以看出，第1组解的误差远远小于第2、第3组，其误差在10-4级别；而另两组解的误差达到了10-2甚至10-1级别，可信度不高。从实际意义上也可以看到，位于陆地的可能性显然高于位于海洋的可能性。因此基于以上分析，地点位于缅甸海、日期为2月18日或10月25日以及3月5日或10月9日的情况仅存在理论上可能；而地点位于新疆维吾尔自治区喀什地区巴楚县，日期为5月24日或7月20日的可能性非常高。

2.对附件3的分析

与上述过程完全类似，由附件2得到的数据见下表（时间已化为十进制）：

北京时间（时）影长（m）

14

表八由附件3计算得到的北京时间与对应影长

用附件3提取的数据运用lsqcurvefit函数，直接得到的可行解如下：

杆长（m）纬度经度赤纬出现次数

10

16

1

1

1

1

1

表九附件3数据第一次迭代后的结果

更改迭代起点为偏离期望值的初始值，进行第二次拟合，得到的结果如下：

杆长（m）纬度经度赤纬出现次数

3

1

表十附件3数据第二次迭代后的结果

可以看出估计结果的出现的次数相差很大，可以认为出现次数仅为一次的估计值可信度过低因此此处仅保留出现次数相对较高的几次结果。对照赤纬表得到对应的日期，通过谷歌地图得到具体的地点，结果如下：

编号纬度经度具体地点日期

1 印度洋5月6日或8月9日

2 湖北省十堰市郧西县2月6日或11月6日

3 陕西省商洛市山阳县2月7日或11月5日

表十一问题三之附件3可能的地点和日期

分别进行误差可视化和再次拟合，得到的图像如下（从上到下依次为编号1、2、3的解）：

图六附件3所有可能解的拟合及误差

由误差图可以看到，得到的这几组地点与日期均与原数据拟合的相当好，理论上都具有很大的可能性；但是考虑到实际情况，仍然认为在陆地上的可能性要大于在海洋上的可能性，因此可以认为可以认为第2组、第3组所得到的地点与日期可能性更大一些。 （四）问题四

为了获取视频中直杆的影长信息，首先调用MATLAB 中的VideoReader ，将附件4导入到MATLAB 中，并通过程序实现每隔1000帧获取一帧图片，随后从中筛选出61张直杆的影子相对清晰的图片，然后筛选出21张等时间间距（2分钟）的图片作为提取影长的信息来源。。

由于摄像机在进行拍摄时使用的是透视角度，其特点是在平面上相互平行的直线在这一视角中延长线将会相交，因此从截取的图片中直接获得的杆长与影长的比例关系可能已经不是其真实的比例关系。因此，在进行杆长测量之前，先使用了PhotoShop 软件对图片进行滤镜处理，垂直透视参数设置为-38，镜头校准为%，使其由透视视角转换为平面视角，进而便于测量影长与杆长之间的比例关系。对于完全相同的一幅图，滤镜处理前后图片的对比如图七所示。

图七滤镜处理前后效果对比图

为了尽可能减小人为因素给模型输出结果带来的误差，我们选择在PhotoShop 软件中，以软件自动建立的坐标来计算影长与杆长之间的像素长度，得到两者之间的比例关系，而非直接在图片上以刻度尺量取，这样做可以最大程度避免人为测量带来的随机误差。设软件上直杆顶端坐标为

11(,)x y ，底端坐标为22(,)x y ，测得某时刻影子顶端坐标为(,)x y ，影子的实际长度为L ，根据比

例关系有 则有

21

L =

（6） 实际测得的数据如下表所示：

北京时间

影子实际长度

表十二视频中提取出不同时刻的实际影长

1.在日期已知的前提下

=°N。此时问题四类似于问题二。将由视频可以得到拍摄日期为7月13日，可得此时赤纬

问题四的已知参数代入问题二的模型。由于所求参数均具有实际意义，因此对它们进行一定的限制，利用非线性拟合函数fit进行拟合，得到的可能取值如下：

编号纬度经度出现次数

1 13

2 7

3 -90 140 1

表十三问题四日期已知情况下得到的可能解

根据误差图误差图以及数据拟合图（图八）可以看到第1、3组数据的拟合效果不好，误差绝对值非常大，可信度均不高，应该被舍去。

图八问题四各可能解的误差图及再次拟合

同样由误差图可以看到，第2组数据（°N，°E）这组解产生的误差相对很小，因此这组解的可信度是很高的。经查，此坐标位于内蒙古自治区乌兰察布市丰镇市。

2.在日期未知的前提下

由于拍摄日期位置，因此无法得知赤纬，所以基于上一问的基础上，增加太阳赤纬的未知数，进行非线性曲线拟合，得到的结果如下：

纬度经度赤纬出现次数

6

………

………

14

-90 100 1

表十四问题四日期未知情况下得到的可能解

作出误差图和数据拟合对比图如图九所示：

图九问题四（日期未知情况下）各可能解的误差图及再次拟合根据作出的误差图和数据拟合对比图，排除掉误差较大的一组解，保留其余两组解，并查的其对应的具体地点和日期如下：

纬度经度赤纬具体地点日期

班达海6月6日或7月7日

内蒙古自治区乌兰察布市丰镇市6月6日或7月8日

表十五问题四（日期未知）的可能地点和日期

尽管这两组解的拟合程度均相当好，但是考虑到实际情况仍然是位于陆地上的可能性高于位于海洋上的可能性，因此测试地点更有可能位于内蒙古自治区乌兰察布市丰镇市，日期为6月6日或7月8日。这一结果与日期已知的情况下得到的参数高度吻合，证明了这一算法的合理性。

六、模型评价

基于非线性曲线拟合函数lsqcurvefit和fit进行系数估算来实现求解,运行速度快，容错性强，可信度高。而且给定不同的初值可能出现相同的解，根据解出现的频率可以评价解的可信度。但是由于拟合的误差，解可能在某个微小范围内波动，而且可能出现附属解。

参考文献

[1]肖志勇，刘宇翔，一种新的纬度测量方法，大学物理，第29卷第9期：52，2010。

[2]徐宝菜，应振华，地球概论教程，北京：高等教育出版社，1983。

[3]TableoftheDeclinationoftheSun，日。

[4]刘浩，韩晶，MATLABR2014a完全自学一本通，北京：电子工业出版社，2015。

附录

本文使用软件为MatlabR2013b、AdobePhotoshopCS6、MicrosoftOfficeExcel2013。其中MatlabR2013b 的程序代码如下（空行表示以下为另一段代码）：

functionbeijingplot()

%问题1绘图函数

t=9::15;

sinh=sin*pi/180)*sin(-(10+5/6)*pi/180)+cos*pi/180)*cos(-(10+5/6)*pi/180)*cos(((t-12)+*(1/15))*2*pi/24); arc=asin(sinh);

arc=pi/2-asin(sinh);

tanh=tan(arc);

tanhl=3*tanh;

plot(t,tanhl);

functiony=fun2(a,t)

%问题2函数

y=a(1)*(1-(a(2)*sin*pi/180)+cos(((t-12)+(a(3)-120)*(1/15))*pi/12)*sqrt(1-a(2)^2)*cos*pi/180)).^2).^(1/2)./( a(2)*sin*pi/180)+cos(((t-12)+(a(3)-120)*(1/15))*pi/12)*sqrt(1-a(2)^2)*cos*pi/180));

fori=-3::3

c0=[iii];

c=lsqcurvefit('fun2',c0,xdata,ydata)；

end

function[ab]=fun2s()

%问题2求解函数迭代起点集中

%a[杆长所求纬度太阳赤纬所求经度]

%b[杆长所求纬度的sin太阳赤纬sin所求经度]

x=xlsread('');

xdata=x(:,1)';

ydata=x(:,2)';

a=[];

options=optimset('Display','off');

fori=-1::1;

c0=[3+3*iii*180];

c=lsqcurvefit('fun2',c0,xdata,ydata,[],[],options);

a=[a;c];

end

b=a;

a(:,2)=asin(a(:,2));

a(:,2)=a(:,2)*180/pi;

c=find(a(:,1)>&imag(a(:,1))==0);

a=a(c,:);

b=b(c,:);

c=find(imag(a(:,2))==0);

a=a(c,:);

b=b(c,:);

c=find(imag(a(:,3))==0);

a=a(c,:);

b=b(c,:);

[cd]=sort(a(:,1));

a=a(d,:);

b=b(d,:);

function[ab]=fun2ss()

%问题2求解函数迭代起点偏离期望值

%a[杆长所求纬度太阳赤纬所求经度]

%b[杆长所求纬度的sin太阳赤纬sin所求经度] x=xlsread('');

xdata=x(:,1)';

ydata=x(:,2)';

a=[];

options=optimset('Display','off');

fori=-10::10;

c0=[iii];

c=lsqcurvefit('fun2',c0,xdata,ydata,[],[],options); a=[a;c];

end

b=a;

a(:,2)=asin(a(:,2));

a(:,2)=a(:,2)*180/pi;

c=find(a(:,1)>&imag(a(:,1))==0);

a=a(c,:);

b=b(c,:);

c=find(imag(a(:,2))==0);

a=a(c,:);

b=b(c,:);

c=find(imag(a(:,3))==0);

a=a(c,:);

b=b(c,:);

[cd]=sort(a(:,1));

a=a(d,:);

b=b(d,:);

functionfun2w(a)

%问题2绘图函数

x=xlsread('');

xdata2=x(:,1)';

ydata2=x(:,2)';

t=::;

y1=fun2(a,t);

subplot(2,1,1);

plot(xdata2,ydata2,'b*-')

holdon

plot(t,y1,'r-')

title('原始数据与拟合数据的对比');

xlabel('T/h');

ylabel('H/m');

legend('原始数据','拟合数据');

y2=y1-ydata2;

subplot(2,1,2);

plot(t,y2,'')

title('误差图');

xlabel('T/h');

ylabel('H/m');

legend('误差');

functionfun2ww(a)

%问题2绘图函数

x=xlsread('');

xdata2=x(:,1)';

ydata2=x(:,2)';

t=9::16;

y1=fun2(a,t);

plot(xdata2,ydata2,'b*-')

holdon

plot(t,y1,'r-')

functiony=fun3(a,t)

%问题3函数

y=a(1)*(1-(a(2)*a(3)+cos(((t-12)+(a(4)-120)*(1/15))*pi/12)*sqrt(1-a(2)^2)*sqrt(1-a(3)^2)).^2).^(1/2)./(a(2)* a(3)+cos(((t-12)+(a(4)-120)*(1/15))*pi/12)*sqrt(1-a(2)^2)*sqrt(1-a(3)^2));

function[ab]=fun32s()

%%问题3附件2求解函数迭代起点集中在期望解附近

%a[杆长所求纬度太阳赤纬所求经度]

%b[杆长所求纬度的sin太阳赤纬sin所求经度]

x=xlsread('');

xdata=x(:,1)';

ydata=x(:,2)';

a=[];

options=optimset('Display','off');

fori=-1::1;

c0=[ii];

c=lsqcurvefit('fun4',c0,xdata,ydata,[],[],options);

a=[a;c];

end

b=a;

a(:,1)=asin(a(:,1));

a(:,1)=a(:,1)*180/pi;

c=find(imag(a(:,1))==0);

a=a(c,:);

b=b(c,:);

c=find(imag(a(:,2))==0);

a=a(c,:);

b=b(c,:);

[cd]=sort(a(:,1));

a=a(d,:);

b=b(d,:);

function[ab]=fun32s()

%问题3附件2求解函数迭代起点集中

%a[杆长所求纬度太阳赤纬所求经度]

%b[杆长所求纬度的sin太阳赤纬sin所求经度]

x=xlsread('');

xdata=x(:,1)';

ydata=x(:,2)';

a=[];

options=optimset('Display','off');

fori=-1::1;

c0=[+*iiii*180];

c=lsqcurvefit('fun3',c0,xdata,ydata,[],[],options);

a=[a;c];

end

b=a;

a(:,2)=asin(a(:,2));

a(:,2)=a(:,2)*180/pi;

a(:,3)=asin(a(:,3));

a(:,3)=a(:,3)*180/pi;

c=find(a(:,1)>&imag(a(:,1))==0);

a=a(c,:);

b=b(c,:);

c=find(imag(a(:,2))==0);

a=a(c,:);

b=b(c,:);

c=find(a(:,3)<&a(:,3)>&imag(a(:,3))==0);

a=a(c,:);

b=b(c,:);

c=find(imag(a(:,4))==0);

a=a(c,:);

b=b(c,:);

[cd]=sort(a(:,1));

a=a(d,:);

b=b(d,:);

function[ab]=fun32ss()

%问题3附件2求解函数迭代起点分散

%a[杆长所求纬度太阳赤纬所求经度]

%b[杆长所求纬度的sin太阳赤纬sin所求经度] x=xlsread('');

xdata=x(:,1)';

ydata=x(:,2)';

a=[];

options=optimset('Display','off');

fori=-10::10;

c0=[+*iiii*180];

c=lsqcurvefit('fun3',c0,xdata,ydata,[],[],options); a=[a;c];

end

b=a;

a(:,2)=asin(a(:,2));

a(:,2)=a(:,2)*180/pi;

a(:,3)=asin(a(:,3));

a(:,3)=a(:,3)*180/pi;

c=find(a(:,1)>&imag(a(:,1))==0);

a=a(c,:);

b=b(c,:);

c=find(imag(a(:,2))==0);

a=a(c,:);

b=b(c,:);

c=find(a(:,3)<&a(:,3)>&imag(a(:,3))==0);

a=a(c,:);

b=b(c,:);

c=find(imag(a(:,4))==0);

a=a(c,:);

b=b(c,:);

[cd]=sort(a(:,1));

a=a(d,:);

b=b(d,:);

functionfun32w(a)

%问题3附件2绘图函数

x=xlsread('');

xdata2=x(:,1)';

ydata2=x(:,2)';

t=::;

y1=fun3(a,t);

subplot(2,1,1);

plot(xdata2,ydata2,'b*-')

holdon

plot(t,y1,'r-')

title('原始数据与拟合数据的对比');

xlabel('T/h');

ylabel('H/m');

legend('原始数据','拟合数据');

y2=y1-ydata2;

subplot(2,1,2);

plot(t,y2,'')

title('误差图');

xlabel('T/h');

ylabel('H/m');

legend('误差');

functionfun32ww(a)

%问题3附件2绘图函数

x=xlsread('');

xdata2=x(:,1)';

ydata2=x(:,2)';

t=9::16;

y1=fun3(a,t);

plot(xdata2,ydata2,'b*-')

holdon

plot(t,y1,'r-')

xlabel('T/h');

ylabel('H/m');

function[ab]=fun33s()

%问题3附件3求解函数迭代起点分散

%a[杆长所求纬度太阳赤纬所求经度]

%b[杆长所求纬度的sin太阳赤纬sin所求经度] x=xlsread('');

xdata=x(:,1)';

ydata=x(:,2)';

a=[];

options=optimset('Display','off');

fori=-10::10;

c0=[+*iiii*180];

c=lsqcurvefit('fun3',c0,xdata,ydata,[],[],options); a=[a;c];

end

b=a;

a(:,2)=asin(a(:,2));

a(:,2)=a(:,2)*180/pi;

a(:,3)=asin(a(:,3));

a(:,3)=a(:,3)*180/pi;

c=find(a(:,1)>&imag(a(:,1))==0);

a=a(c,:);

b=b(c,:);

c=find(imag(a(:,2))==0);

a=a(c,:);

b=b(c,:);

c=find(a(:,3)<&a(:,3)>&imag(a(:,3))==0);

a=a(c,:);

b=b(c,:);

c=find(imag(a(:,4))==0);

a=a(c,:);

b=b(c,:);

[cd]=sort(a(:,1));

a=a(d,:);

b=b(d,:);

function[ab]=fun33s()

%问题3附件2求解函数迭代起点在期望值附近%a[杆长所求纬度太阳赤纬所求经度]

%b[杆长所求纬度的sin太阳赤纬sin所求经度] x=xlsread('');

xdata=x(:,1)';

ydata=x(:,2)';

a=[];

options=optimset('Display','off');

fori=-1::1;

c0=[+*ii-ii*180];

c=lsqcurvefit('fun3',c0,xdata,ydata,[],[],options); a=[a;c];

end

b=a;

a(:,2)=asin(a(:,2));

a(:,2)=a(:,2)*180/pi;

a(:,3)=asin(a(:,3));

a(:,3)=a(:,3)*180/pi;

c=find(a(:,1)>&imag(a(:,1))==0);

a=a(c,:);

b=b(c,:);

c=find(imag(a(:,2))==0);

a=a(c,:);

b=b(c,:);

c=find(a(:,3)<&a(:,3)>&imag(a(:,3))==0);

a=a(c,:);

b=b(c,:);

c=find(imag(a(:,4))==0);

a=a(c,:);

b=b(c,:);

[cd]=sort(a(:,1));

a=a(d,:);

b=b(d,:);

functionfun33w(a)

%问题3附件3绘图函数

x=xlsread('');

xdata2=x(:,1)';

ydata2=x(:,2)';

t=::;

y1=fun3(a,t);

subplot(2,1,1);

plot(xdata2,ydata2,'b*-')

holdon

plot(t,y1,'r-')

title('原始数据与拟合数据的对比'); xlabel('T/h');

ylabel('H/m');

legend('原始数据','拟合数据');

y2=y1-ydata2;

subplot(2,1,2);

plot(t,y2,'')

title('误差图');

xlabel('T/h');

ylabel('H/m');

legend('误差');

functionfun32ww(a)

%问题3附件2绘图函数

x=xlsread('');

xdata2=x(:,1)';

ydata2=x(:,2)';

t=9::16;

y1=fun3(a,t);

plot(xdata2,ydata2,'b*-')

holdon

plot(t,y1,'r-')

xlabel('T/h');

ylabel('H/m');

functionfun4video()

%视频帧提取函数

obj=VideoReader('');

numFrames=;

fork=1:1000:numFrames

frame=read(obj,k);

imshow(frame);

imwrite(frame,strcat(num2str(k),'.jpg'),'jpg'); end

function[ef]=fun4f()

%第四题第二问求解函数

%e[所求纬度所求经度]

%f[所求纬度sin所求精度]

x=xlsread('');

xdata2=x(:,1)';

ydata2=x(:,2)';

e=[];

s=[-1-180];

fori=0:20

s=[-1+i*+i*18];

options=fitoptions('Method','NonlinearLeastSquares','Lower',[-1100],'Upper',[1180],'StartPoint',s,'Display', 'off');

z=fittype('2*(1-(a*sin*pi/180)+cos(((x-12)+(c-120)*(1/15))*pi/12)*sqrt(1-a^2)*cos*pi/180)).^2).^(1/2)./(a*si n*pi/180)+cos(((x-12)+(c-120)*(1/15))*pi/12)*sqrt(1-a^2)*cos*pi/180))','options',options);

[abc]=fit(xdata2',ydata2',z);

d=[第四题第二问求解函数

%e[所求纬度所求经度]

%f[所求纬度sin所求精度]

x=xlsread('');

xdata2=x(:,1)';

ydata2=x(:,2)';

e=[];

s=[];

fori=0:20

s=[-1+i*'Method','NonlinearLeastSquares','Lower',[],'Upper',[],'StartPoint',s,'Display','off');

z=fittype('2*(1-(a*b+cos(((x-12)+(c-120)*(1/15))*pi/12)*sqrt(1-a^2)*sqrt(1-b^2)).^2).^(1/2)./(a*b+cos(((x-1 2)+(c-120)*(1/15))*pi/12)*sqrt(1-a^2)*sqrt(1-b^2))','options',options);

[abc]=fit(xdata2',ydata2',z);

d=[问题4绘图函数

x=xlsread('');

xdata2=x(:,1)';

ydata2=x(:,2)';

t=xdata2;

y1=fun4(a,t);

subplot(2,1,1);

plot(xdata2,ydata2,'b*-')

holdon

plot(t,y1,'r-')

title('原始数据与拟合数据的对比');

xlabel('T/h');

ylabel('H/m');

legend('原始数据','拟合数据');

y2=y1-ydata2;

subplot(2,1,2);

plot(t,y2,'')

title('误差图');

xlabel('T/h');

ylabel('H/m');

legend('误差');

functionfun4ww(a)

%问题2绘图函数

x=xlsread('');

xdata2=x(:,1)';

ydata2=x(:,2)';

t=9::15;

y1=fun4(a,t);

plot(xdata2,ydata2,'b*-')

holdon

plot(t,y1,'r-')

xlabel('T/h')

ylabel('H/m')