2024年2月27日发(作者：扬州小学期中数学试卷)
基于非线性曲线拟合的经纬度测量方法
摘要
本文首先基于天体物理学知识，构造出地球上某处直杆的影长与时间的函数关系式；然后运用非线性曲线拟合的方法，求解缺省参数，再根据直杆影长的变化规律，推算出测量点的地理位置及所处的日期。
在问题一中，本文以北京时间为参考时间，对地球上某一点处直杆影长的影响因素进行分析，发现其与直杆所处纬度、太阳直射点处纬度、所处时刻及经度等因素有关，结合地理知识构造出影长与影响因素的函数关系式。在各项参数均已给定的情况下，即可作出题目所要求的影长-时间变化曲线。
对于问题二，本文由附件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
一、 问题重述
如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。
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张截图。在经过滤镜处理后，可以由截图得到影长与杆长的比例关系，从而得到不同时间下的影长。在日期已知的情况下，问题可以利用问题二的模型；若日期未知，则可以利用问题三的模型，得到所求参量。
2
三、模型假设
1. 地球是半径为R的均匀圆球，即忽略其表面地形及地球微小的椭球率带来的影响；
2. 忽略太阳光在穿过地球大气层时产生的折射；
3. 照射到地球的太阳光为平行光；
4. 黄赤交角为一固定值。
四、 符号约定
0
l
太阳直射点纬度
测量点的纬度
测量点的经度
测量点的时角
太阳直射方向向量
地心指向测量点的方向向量
太阳高度
直杆长
影长
测量点地方时
测量点地方时对应的北京时间
a
e1
e2
h
L0
L
t
t0
五、模型建立与求解
（一）问题一
根据假设1，以如图所示的方式建立空间直角坐标系[1]。其中，以地心为坐标原点，垂直于赤道平面指向正北的轴为z轴，沿赤道平面指向太阳方向的轴为x轴，与二者成左手关系的轴为y轴。假定此时太阳直射点为地球上一点A，设其纬度为0（北纬为正，南纬为负，以下涉及纬度符号均同理）。显然，由于x轴指向太阳方向，因此直射点A位于oxz平面内。设此时刻太阳直射方向为e1。设测量点纬度为，其时角[2]为a，其与地心构成的向量为e2。
3
ze2观测点yφae1ARφ0x赤道
图一 模型推导示意图
根据以上设定，有
e1(Rcos0cosa,0,Rsin0)
e2(Rcoscosa,Rcossina,Rsin)
可求得e1与e2之间的空间的余弦值为
cose1e2|e||ecoscos0cosasinsin0
12|根据太阳高度角h的定义，有h2，因此可得
sinhcoscos0cosasinsin0
由图二所示的示意图可得，影长、直杆长之间存在关系
LL0tanh
4
1）2）3）
（
（
（
阳光直杆长L0h影长L图二 影长、杆长与太阳高度的关系
在上述推导中，时角a(t12)
12，其中t为测量点地方时。考虑到附录中所给时间均为北京时间，因此此处统一将地方时转化为北京时间表示。设测量点经度为l，为方便起见，不妨假设测量点在北京所在的东八区的中间经线（120°E）以东，因此二者之间经度差为
ll120
由于地球每自转24h便转过360°，因此近似可认为经度上每度引起的时差为4min，由前可知测量点在120°E以东，因此测量点地方时与东八区区时有如下关系
tt0因此
(l120)
15
at0(l120)12 （4）
1512其中t0为测量点地方时t对应的北京时间。
联立式（2）、（3）、（4），可得地球上某一点处直杆的影长与北京时间的函数关系式
(l120)12sinsin01coscos0cost01512 （5）
LL0(l120)coscos0cost012sinsin01512查表[3]可知2015年10月22日太阳直射点010.833，并将已知条件L03 m，5
2
39.9072，l116.3914E，使用MATLAB即可作出题目给定时间内天安门广场上3m长的直杆的影长变化曲线，如图三所示。
图三 问题一所求曲线
（二）问题二
由附件可以得到测量点直杆的影长与北京时间的对应数据，可以在问题一所得到的模型的基础上，将待求的经、纬度坐标以及杆长看作为模型的参数，调用MATLAB中lsqcurvefit[4]函数，通过迭代的方法实现最小二乘拟合，得到参数的估计值。
Lsqcurvefit函数是进行最小二乘非线性曲线拟合的函数，设定初始迭代起点进行拟合求最优解。
为了更快得到最优解，将迭代起点选择为期望值附近（但并未限定所求解的范围），即杆长期望在0至3米之间，经度位于-180°至180°之间，纬度的正弦值位于-1到1之间。
求解完成后，基于以下原则，在MATLAB返回的可行解当中进行了初步的筛选：
1. 所有参数均应为纯实数，即虚部应等于零；
2. 出于实际考虑，要求杆长均大于0.5m；
3. 参数完全重复的或相当接近的，取其中拟合程度最佳的一组参数。
由附件1得到的数据见下表（时间已化为十进制）：
北京时间（时）
14.7
14.75
14.8
6
影长（m）
1.149625826
1.182198976
1.215296955
14.85
14.9
14.95
15
15.05
15.1
15.15
15.2
15.25
15.3
15.35
15.4
15.45
15.5
15.55
15.6
15.65
15.7
1.249051052
1.28319534
1.317993149
1.353364049
1.389387091
1.426152856
1.463399853
1.501481622
1.540231817
1.579853316
1.620144515
1.661270613
1.703290633
1.74620591
1.790050915
1.835014272
1.880875001
1.927918447
表一 由附件1计算得到的北京时间与对应影长
具体程序见附录，在代入由附件1得到的数据之后，直接得到的可行解如下：
杆长（m）
1.9937
2.1479
…
2.2523
2.4734
2.6254
2.794
…
2.7953
纬度
18.98
-3.0152
…
-3.6179
23.369
24.185
24.819
…
24.823
经度
109.35
104.32
…
102.66
101.98
99.892
97.757
…
97.74
34
1
1
6
出现次数
30
表二 附件1数据第一次迭代后的结果
由表格可以看到某些解出现的次数不止一次，这是因为不同点起点得到相同的结果，说明该结果可能为最优解。
由表格知（98°N，109.35°E）在可行解中出现的频率极高，说明其可信度相对于其他解要更高；而在（23.369°~24.823°N，101.98°~97.74°E）这一范围中，尽管每组解出现的次数不高，但是可以发现这些数据分布相当密集，因此可以看做是同一个地点；（3.0152°~3.6179°S，7
104.32~102.66°E）范围之间的解同理，但是由于其出现次数较少，显然其可信度不及另两组解。
为了增强以上结果的可信度，现改变迭代起点进行第二次拟合。在这一次拟合中，将迭代起点设置偏离期望值。得到的结果如下表所示：
杆长（m）
1.9937
1.5114
纬度
18.98
经度
109.3494
出现次数
1
1 -2.6700 114.7980
表三 附件1数据第二次迭代后的结果
（18.98°N，109.35°E）这一组解再次出现，证明了这一组解得可信度确实相当高。将以上得到的可能解汇总，并确定具体的地点如下：
编号
1
2
3
4
杆长（m）
1.9937
2.7945
2.2523
1.5114
纬度
18.98
24.821
-3.6179
-2.6700
经度
109.3494
97.75
102.66
114.7980
具体地点
海南省白沙黎族自治县
云南省德宏傣族景颇族自治州盈江县
印度尼西亚Sinar Gunung Tebat Karai Kepahiang Regency
印度尼西亚Asia Baru KuripanBarito Kuala Regency
表四 问题二可能的地点
分别对得到的四组解进行误差可视化和再次拟合，图像如下（从上到下依次为编号1、2、3、4的解）：
8
图四 问题二所有可能解的拟合及误差
误差图显示拟合值与真实值之间的误差均为10-3级别，这再一次证明了解的可信度。这是基于本模型得到的四个可能的地点。考虑到迭代过程中出现次数悬殊较大，有理由相信海南省白沙黎族自治县（18.98°N，109.3494°E）是最有可能的地点，而其他三个地点可能性相对较小。
（三）问题三
第三问在第二问的基础上又增添了日期这个未知参数。由于日期与赤纬存在一定的对应关系，同第二问思路相同，我们将赤纬作为另一个进lsqcurvefit拟合，此时初步筛选原则除去问题二中提9
出的3条外，再增添一条：赤纬的范围需位于23.4333°S~23.4333°N之间。
1. 对附件2的分析
由附件2得到的数据见下表（时间已化为十进制）：
北京时间（时）
12.6833
12.7333
12.7833
12.8333
12.8833
12.9333
12.9833
13.0333
13.0833
13.1333
13.1833
13.2333
13.2833
13.3333
13.3833
13.4333
13.4833
13.5333
13.5833
13.6333
13.6833
影长（m）
1.247256205
1.22279459
1.198921486
1.175428964
1.152439573
1.12991747
1.10783548
1.086254206
1.065081072
1.044446265
1.024264126
1.004640314
0.985490908
0.966790494
0.948584735
0.930927881
0.91375175
0.897109051
0.880973762
0.865492259
0.850504468
表五 由附件2计算得到的北京时间与对应影长
将附件2提取的数据代入模型，直接得到的可行解如下：
杆长（m）
2.0008
2.3484
3.8662
纬度
39.893
10.953
7.3544
经度
79.744
97.605
97.727
赤纬
20.761
-11.902
-6.2274
出现次数
10
1
1
表五 附件2数据第一次迭代后的结果
同样，为了增强解得可信度，使迭代起点偏离期望住进行了第二次拟合，直接得到的可行解如下：
杆长（m）
2.0008
纬度 经度 赤纬 出现次数
1 39.8926 79.7443 20.7606
10
表六 附件2数据第二次迭代后的结果
由此可见，第二次得到的这组解可信度很高。将表四中三组解依次编号为1、2、3，以赤纬进行查表确定日期，最终得到下列可能的时间和地点：
编号
1
2
3
纬度 经度 具体地点 日期
5月24日或7月20日
2月18日或10月25日
3月5日或10月9日
39.893 79.744
新疆维吾尔自治区喀什地区巴楚县
10.953 97.605
7.3544 97.727
缅甸海
缅甸海
表七 问题三之附件2可能的地点和日期
分别进行误差可视化，得到的图像如下（从上到下依次为编号1、2、3的解）：
11
图五 附件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得到的数据见下表（时间已化为十进制）：
北京时间（时）
13.15
13.2
13.25
13.3
13.35
13.4
13.45
13.5
13.55
13.6
13.65
13.7
13.75
13.8
13.85
13.9
13.95
14
14.05
14.1
12
影长（m）
3.533142184
3.546768029
3.561797643
3.578100715
3.595750783
3.61493428
3.635425983
3.657218272
3.680541115
3.705167836
3.731278025
3.758917911
3.788087888
3.818701015
3.850809619
3.88458522
3.919911828
3.956875992
3.99553479
4.035750835
14.15 4.077863059
表八 由附件3计算得到的北京时间与对应影长
用附件3提取的数据运用lsqcurvefit函数，直接得到的可行解如下：
杆长（m）
3.0356
2.3484
3.8662
3.0371
3.0378
4.5644
5.5373
纬度
-32.849
32.849
7.3544
33.06
33.35
-18.866
-16.421
经度
110.24
110.24
97.727
110.27
110.29
96.681
103.8
赤纬
15.959
-15.959
-6.2274
-15.731
-15.43
20.331
16.387
出现次数
10
16
1
1
1
1
1
表九 附件3数据第一次迭代后的结果
更改迭代起点为偏离期望值的初始值，进行第二次拟合，得到的结果如下：
杆长（m）
3.0356
3.0391
纬度
32.8488
33.5224
经度
110.2450
110.3023
赤纬
-15.9593
-15.2465
出现次数
3
1
表十 附件3数据第二次迭代后的结果
可以看出估计结果的出现的次数相差很大，可以认为出现次数仅为一次的估计值可信度过低因此此处仅保留出现次数相对较高的几次结果。对照赤纬表得到对应的日期，通过谷歌地图得到具体的地点，结果如下：
编号
1
2
3
纬度 经度 具体地点
印度洋
日期
5月6日或8月9日
-32.849 110.24
32.849 110.24
湖北省十堰市郧西县 2月6日或11月6日
33.35 110.29
陕西省商洛市山阳县 2月7日或11月5日
13
表十一 问题三之附件3可能的地点和日期
分别进行误差可视化和再次拟合，得到的图像如下（从上到下依次为编号1、2、3的解）：
图六 附件3所有可能解的拟合及误差
由误差图可以看到，得到的这几组地点与日期均与原数据拟合的相当好，理论上都具有很大的可能性；但是考虑到实际情况，仍然认为在陆地上的可能性要大于在海洋上的可能性，因此可以认为可以认为第2组、第3组所得到的地点与日期可能性更大一些。
（四）问题四
14
为了获取视频中直杆的影长信息，首先调用MATLAB中的VideoReader，将附件4导入到MATLAB中，并通过程序实现每隔1000帧获取一帧图片，随后从中筛选出61张直杆的影子相对清晰的图片，然后筛选出21张等时间间距（2分钟）的图片作为提取影长的信息来源。。
由于摄像机在进行拍摄时使用的是透视角度，其特点是在平面上相互平行的直线在这一视角中延长线将会相交，因此从截取的图片中直接获得的杆长与影长的比例关系可能已经不是其真实的比例关系。因此，在进行杆长测量之前，先使用了PhotoShop软件对图片进行滤镜处理，垂直透视参数设置为-38，镜头校准为49.1%，使其由透视视角转换为平面视角，进而便于测量影长与杆长之间的比例关系。对于完全相同的一幅图，滤镜处理前后图片的对比如图七所示。
滤镜处理后图七 滤镜处理前后效果对比图
为了尽可能减小人为因素给模型输出结果带来的误差，我们选择在PhotoShop软件中，以软件自动建立的坐标来计算影长与杆长之间的像素长度，得到两者之间的比例关系，而非直接在图片上以刻度尺量取，这样做可以最大程度避免人为测量带来的随机误差。设软件上直杆顶端坐标为
(x1,y1)，底端坐标为(x2,y2)，测得某时刻影子顶端坐标为(x,y)，影子的实际长度为L，根据比例关系有
(xx2)2(yy2)2y2y1
L2则有
2(xx2)2(yy2)2
L （6）
y2y1实际测得的数据如下表所示：
北京时间
8.9128
8.9461
8.9794
9.0127
9.046
9.0793
9.1126
15
影子实际长度
2.2256
2.2009
2.1718
2.1479
2.1180
2.0889
2.0680
9.1459
9.1792
9.2125
9.2458
9.2791
9.3124
9.3457
9.379
9.4123
9.4456
9.4789
9.5122
9.5455
9.5788
2.0403
2.0135
1.9911
1.9530
1.9411
1.9075
1.8858
1.8597
1.8351
1.8098
1.7859
1.7575
1.7337
1.7150
表十二 视频中提取出不同时刻的实际影长
1. 在日期已知的前提下
由视频可以得到拍摄日期为7月13日，可得此时赤纬0=21.9167°N。此时问题四类似于问题二。将问题四的已知参数代入问题二的模型。由于所求参数均具有实际意义，因此对它们进行一定的限制，利用非线性拟合函数fit进行拟合，得到的可能取值如下：
编号
1
2
3
纬度
-6.5203
40.3735
-90
经度
126.5685
113.2984
140
出现次数
13
7
1
表十三 问题四日期已知情况下得到的可能解
根据误差图误差图以及数据拟合图（图八）可以看到第1、3组数据的拟合效果不好，误差绝对值非常大，可信度均不高，应该被舍去。
16
图八 问题四各可能解的误差图及再次拟合
同样由误差图可以看到，第2组数据（40.3735°N，113.2984°E）这组解产生的误差相对很小，因此这组解的可信度是很高的。经查，此坐标位于内蒙古自治区乌兰察布市丰镇市。
2. 在日期未知的前提下
由于拍摄日期位置，因此无法得知赤纬，所以基于上一问的基础上，增加太阳赤纬的未知数，进行非线性曲线拟合，得到的结果如下：
纬度
40.5311
…
40.6507
-5.3519
…
-5.7684
-90
经度
112.8863
…
112.5483
126.3530
…
126.4333
0.39
赤纬
22.4886
…
22.9545
22.9545
…
22.5879
100 1
14
6
出现次数
表十四 问题四日期未知情况下得到的可能解
作出误差图和数据拟合对比图如图九所示：
17
图九 问题四（日期未知情况下）各可能解的误差图及再次拟合
根据作出的误差图和数据拟合对比图，排除掉误差较大的一组解，保留其余两组解，并查的其对应的具体地点和日期如下：
纬度 经度 赤纬 具体地点
班达海
日期
6月6日或7月7日
-5.6845 126.4174 22.6624
40.5578 112.8126 22.5879
内蒙古自治区乌兰察布市丰镇市 6月6日或7月8日
表十五 问题四（日期未知）的可能地点和日期
尽管这两组解的拟合程度均相当好，但是考虑到实际情况仍然是位于陆地上的可能性高于位于海洋上的可能性，因此测试地点更有可能位于内蒙古自治区乌兰察布市丰镇市，日期为6月6日或7月8日。这一结果与日期已知的情况下得到的参数高度吻合，证明了这一算法的合理性。
六、模型评价
基于非线性曲线拟合函数lsqcurvefit和fit进行系数估算来实现求解,运行速度快，容错性强，可信度高。而且给定不同的初值可能出现相同的解，根据解出现的频率可以评价解的可信度。但是由于拟合的误差，解可能在某个微小范围内波动，而且可能出现附属解。
18
参考文献
[1]肖志勇，刘宇翔，一种新的纬度测量方法，大学物理，第29卷第9期：52，2010。
[2]徐宝菜，应振华，地球概论教程，北京：高等教育出版社，1983。
[3] Table of the Declination of the Sun，/~wsanford/exo/sundials/DEC_，2015年9月13日。
[4]刘浩，韩晶，MATLAB R2014a完全自学一本通，北京：电子工业出版社，2015。
附录
本文使用软件为Matlab R2013b、Adobe Photoshop CS6、Microsoft Office Excel 2013。其中Matlab
R2013b的程序代码如下（空行表示以下为另一段代码）：
function beijingplot()
%问题1绘图函数
t=9:0.01:15;
sinh=sin(39.907*pi/180)*sin(-(10+5/6)*pi/180)+cos(39.907*pi/180)*cos(-(10+5/6)*pi/180)*cos(((t-12)+(116.3914-120)*(1/15))*2*pi/24);
arc=asin(sinh);
arc=pi/2-asin(sinh);
tanh=tan(arc);
tanhl=3*tanh;
plot(t,tanhl);
function y=fun2(a,t)
%问题2函数
y=a(1)*(1 -(a(2)*sin(10.617*pi/180) + cos(((t-12)+(a(3)-120)*(1/15))*pi/12)*sqrt(1-a(2)^2)*cos(10.617*pi/180)).^2).^(1/2)./(a(2)*sin(10.617*pi/180) + cos(((t-12)+(a(3)-120)*(1/15))*pi/12)*sqrt(1-a(2)^2)*cos(10.617*pi/180));
for i=-3:0.1:3
c0=[i i i];
c=lsqcurvefit(\'fun2\',c0,xdata,ydata)；
end
function [a b]=fun2s()
%问题2求解函数 迭代起点集中
%a [杆长 所求纬度 太阳赤纬 所求经度]
%b [杆长 所求纬度的sin 太阳赤纬sin 所求经度]
x=xlsread(\'\');
xdata=x(:,1)\';
ydata=x(:,2)\';
a=[];
options = optimset(\'Display\',\'off\');
for i=-1:0.01:1;
19
c0=[3+3*i i i*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)>0.5&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,:);
[c d]=sort(a(:,1));
a=a(d,:);
b=b(d,:);
function [a b]=fun2ss()
%问题2求解函数 迭代起点偏离期望值
%a [杆长 所求纬度 太阳赤纬 所求经度]
%b [杆长 所求纬度的sin 太阳赤纬sin 所求经度]
x=xlsread(\'\');
xdata=x(:,1)\';
ydata=x(:,2)\';
a=[];
options = optimset(\'Display\',\'off\');
for i=-10:0.1:10;
c0=[i i i];
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)>0.5&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,:);
[c d]=sort(a(:,1));
a=a(d,:);
b=b(d,:);
function fun2w(a)
%问题2绘图函数
x=xlsread(\'\');
xdata2=x(:,1)\';
20
ydata2=x(:,2)\';
t=14.7:0.05:15.7;
y1=fun2(a,t);
subplot(2,1,1);
plot(xdata2,ydata2,\'b*-\')
hold on
plot(t,y1,\'r-\')
title(\'原始数据与拟合数据的对比\');
xlabel(\'T/h\');
ylabel(\'H/m\');
legend(\'原始数据\',\'拟合数据\');
y2=y1-ydata2;
subplot(2,1,2);
plot(t,y2,\'r.-\')
title(\'误差图\');
xlabel(\'T/h\');
ylabel(\'H/m\');
legend(\'误差\');
function fun2ww(a)
%问题2绘图函数
x=xlsread(\'\');
xdata2=x(:,1)\';
ydata2=x(:,2)\';
t=9:0.01:16;
y1=fun2(a,t);
plot(xdata2,ydata2,\'b*-\')
hold on
plot(t,y1,\'r-\')
function y=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 [a b]=fun32s()
%%问题3附件2求解函数 迭代起点集中在期望解附近
%a [杆长 所求纬度 太阳赤纬 所求经度]
%b [杆长 所求纬度的sin 太阳赤纬sin 所求经度]
x=xlsread(\'\');
xdata=x(:,1)\';
ydata=x(:,2)\';
a=[];
options = optimset(\'Display\',\'off\');
for i=-1:0.01:1;
c0=[ i i ];
c=lsqcurvefit(\'fun4\',c0,xdata,ydata,[],[],options);
a=[a;c];
end
b=a;
a(:,1)=asin(a(:,1));
a(:,1)=a(:,1)*180/pi;
21
c=find(imag(a(:,1))==0);
a=a(c,:);
b=b(c,:);
c=find(imag(a(:,2))==0);
a=a(c,:);
b=b(c,:);
[c d]=sort(a(:,1));
a=a(d,:);
b=b(d,:);
function [a b]=fun32s()
%问题3附件2求解函数 迭代起点集中
%a [杆长 所求纬度 太阳赤纬 所求经度]
%b [杆长 所求纬度的sin 太阳赤纬sin 所求经度]
x=xlsread(\'\');
xdata=x(:,1)\';
ydata=x(:,2)\';
a=[];
options = optimset(\'Display\',\'off\');
for i=-1:0.01:1;
c0=[2.5+2.5*i i i i*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)>0.5&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)<23.4393&a(:,3)>-23.4393&imag(a(:,3))==0);
a=a(c,:);
b=b(c,:);
c=find(imag(a(:,4))==0);
a=a(c,:);
b=b(c,:);
[c d]=sort(a(:,1));
a=a(d,:);
b=b(d,:);
function [a b]=fun32ss()
%问题3附件2求解函数 迭代起点分散
%a [杆长 所求纬度 太阳赤纬 所求经度]
%b [杆长 所求纬度的sin 太阳赤纬sin 所求经度]
x=xlsread(\'\');
xdata=x(:,1)\';
ydata=x(:,2)\';
22
a=[];
options = optimset(\'Display\',\'off\');
for i=-10:0.1:10;
c0=[2.5+2.5*i i i i*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)>0.5&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)<23.4393&a(:,3)>-23.4393&imag(a(:,3))==0);
a=a(c,:);
b=b(c,:);
c=find(imag(a(:,4))==0);
a=a(c,:);
b=b(c,:);
[c d]=sort(a(:,1));
a=a(d,:);
b=b(d,:);
function fun32w(a)
%问题3附件2绘图函数
x=xlsread(\'\');
xdata2=x(:,1)\';
ydata2=x(:,2)\';
t=12.6833:0.05:13.6833;
y1=fun3(a,t);
subplot(2,1,1);
plot(xdata2,ydata2,\'b*-\')
hold on
plot(t,y1,\'r-\')
title(\'原始数据与拟合数据的对比\');
xlabel(\'T/h\');
ylabel(\'H/m\');
legend(\'原始数据\',\'拟合数据\');
y2=y1-ydata2;
subplot(2,1,2);
plot(t,y2,\'r.-\')
title(\'误差图\');
xlabel(\'T/h\');
ylabel(\'H/m\');
legend(\'误差\');
function fun32ww(a)
%问题3附件2绘图函数
23
x=xlsread(\'\');
xdata2=x(:,1)\';
ydata2=x(:,2)\';
t=9:0.01:16;
y1=fun3(a,t);
plot(xdata2,ydata2,\'b*-\')
hold on
plot(t,y1,\'r-\')
xlabel(\'T/h\');
ylabel(\'H/m\');
function [a b]=fun33s()
%问题3附件3求解函数 迭代起点分散
%a [杆长 所求纬度 太阳赤纬 所求经度]
%b [杆长 所求纬度的sin 太阳赤纬sin 所求经度]
x=xlsread(\'\');
xdata=x(:,1)\';
ydata=x(:,2)\';
a=[];
options = optimset(\'Display\',\'off\');
for i=-10:0.1:10;
c0=[2.5+2.5*i i i i*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)>0.5&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)<23.4393&a(:,3)>-23.4393&imag(a(:,3))==0);
a=a(c,:);
b=b(c,:);
c=find(imag(a(:,4))==0);
a=a(c,:);
b=b(c,:);
[c d]=sort(a(:,1));
a=a(d,:);
b=b(d,:);
function [a b]=fun33s()
%问题3附件2求解函数 迭代起点在期望值附近
%a [杆长 所求纬度 太阳赤纬 所求经度]
%b [杆长 所求纬度的sin 太阳赤纬sin 所求经度]
x=xlsread(\'\');
xdata=x(:,1)\';
ydata=x(:,2)\';
24
a=[];
options = optimset(\'Display\',\'off\');
for i=-1:0.01:1;
c0=[2.5+2.5*i i -i i*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)>0.5&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)<23.4393&a(:,3)>-23.4393&imag(a(:,3))==0);
a=a(c,:);
b=b(c,:);
c=find(imag(a(:,4))==0);
a=a(c,:);
b=b(c,:);
[c d]=sort(a(:,1));
a=a(d,:);
b=b(d,:);
function fun33w(a)
%问题3附件3绘图函数
x=xlsread(\'\');
xdata2=x(:,1)\';
ydata2=x(:,2)\';
t=13.15:0.05:14.15;
y1=fun3(a,t);
subplot(2,1,1);
plot(xdata2,ydata2,\'b*-\')
hold on
plot(t,y1,\'r-\')
title(\'原始数据与拟合数据的对比\');
xlabel(\'T/h\');
ylabel(\'H/m\');
legend(\'原始数据\',\'拟合数据\');
y2=y1-ydata2;
subplot(2,1,2);
plot(t,y2,\'r.-\')
title(\'误差图\');
xlabel(\'T/h\');
ylabel(\'H/m\');
legend(\'误差\');
function fun32ww(a)
%问题3附件2绘图函数
25
x=xlsread(\'\');
xdata2=x(:,1)\';
ydata2=x(:,2)\';
t=9:0.01:16;
y1=fun3(a,t);
plot(xdata2,ydata2,\'b*-\')
hold on
plot(t,y1,\'r-\')
xlabel(\'T/h\');
ylabel(\'H/m\');
function fun4video()
%视频帧提取函数
obj=VideoReader(\'\');
numFrames = OfFrames;
for k = 1 :1000: numFrames
frame = read(obj,k);
imshow(frame);
imwrite(frame,strcat(num2str(k),\'.jpg\'),\'jpg\');
end
function [e f]=fun4f()
%第四题第二问求解函数
%e [所求纬度 所求经度]
%f [所求纬度sin 所求精度]
x=xlsread(\'\');
xdata2=x(:,1)\';
ydata2=x(:,2)\';
e=[];
s=[-1 -180];
for i=0:20
s=[-1+i*0.1 -180+i*18];
options=fitoptions(\'Method\',\'NonlinearLeastSquares\',\'Lower\',[-1 100],\'Upper\',[1
180],\'StartPoint\',s,\'Display\',\'off\');
z=fittype(\'2*(1 -(a*sin(21.9167*pi/180) + cos(((x-12)+(c-120)*(1/15))*pi/12)*sqrt(1-a^2)*cos(21.9167*pi/180)).^2).^(1/2)./(a*sin(21.9167*pi/180) + cos(((x-12)+(c-120)*(1/15))*pi/12)*sqrt(1-a^2)*cos(21.9167*pi/180))\',\'options\',options);
[a b c]=fit(xdata2\',ydata2\',z);
d=[a.a a.c];
e=[e;d];
end
f=e;
e(:,1)=asin(e(:,1));
e(:,1)=e(:,1)*180/pi;
function [e f]=fun42f()
%第四题第二问求解函数
%e [所求纬度 所求经度]
%f [所求纬度sin 所求精度]
26
x=xlsread(\'\');
xdata2=x(:,1)\';
ydata2=x(:,2)\';
e=[];
s=[-1 -3.9 -180];
for i=0:20
s=[-1+i*0.1 -3.9+i*0.039 -180+i*18];
options=fitoptions(\'Method\',\'NonlinearLeastSquares\',\'Lower\',[-1 0.38 100],\'Upper\',[1 0.39
180],\'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-12)+(c-120)*(1/15))*pi/12)*sqrt(1-a^2)*sqrt(1-b^2))\',\'options\',options);
[a b c]=fit(xdata2\',ydata2\',z);
d=[a.a a.b a.c];
e=[e;d];
end
f=e;
e(:,1)=asin(e(:,1));
e(:,1)=e(:,1)*180/pi;
function fun4w(a)
%问题4绘图函数
x=xlsread(\'\');
xdata2=x(:,1)\';
ydata2=x(:,2)\';
t=xdata2;
y1=fun4(a,t);
subplot(2,1,1);
plot(xdata2,ydata2,\'b*-\')
hold on
plot(t,y1,\'r-\')
title(\'原始数据与拟合数据的对比\');
xlabel(\'T/h\');
ylabel(\'H/m\');
legend(\'原始数据\',\'拟合数据\');
y2=y1-ydata2;
subplot(2,1,2);
plot(t,y2,\'r.-\')
title(\'误差图\');
xlabel(\'T/h\');
ylabel(\'H/m\');
legend(\'误差\');
function fun4ww(a)
%问题2绘图函数
x=xlsread(\'\');
xdata2=x(:,1)\';
ydata2=x(:,2)\';
t=9:0.01:15;
y1=fun4(a,t);
plot(xdata2,ydata2,\'b*-\')
hold on
27
plot(t,y1,\'r-\')
xlabel(\'T/h\')
ylabel(\'H/m\')
28
更多推荐
得到,可能,时间,太阳,影长,附件,所求,问题
发布评论