1、1基于非线性曲线拟合的经纬度测量方法摘要本文首先基于天体物理学知识，构造出地球上某处直杆的影长与时间的函数关系式；然后运用非线性曲线拟合的方法，求解缺省参数，再根据直杆影长的变化规律，推算出测量点的地理位置及所处的日期。在问题一中，本文以北京时间为参考时间，对地球上某一点处直杆影长的影响因素进行分析，发现其与直杆所处纬度、太阳直射点处纬度、所处时刻及经度等因素有关，结合地理知识构造出影长与影响因素的函数关系式。在各项参数均已给定的情况下，即可作出题目所要求的影长-时间变化曲线。对于问题二，本文由附件 1 给定的时刻及其影长，运用非线性曲线拟合的方法，利用问题一中建立的关系式，将时间与影长作为已
2、知参数，利用 lsqcurvefit 函数拟合求解经纬度参数。联系实际，筛选出可能的 4 个位置，并认为海南省白沙黎族自治县是最有可能的地点。问题三与问题二基本相似，本文仍然在附件所得的数据基础上进行 lsqcurvefit 非线性曲线拟合，得到经度、纬度以及赤纬的可行解，根据所求赤纬，通过查表可以得到可能的日期。由附件 2 得到3 个可能的地点与 6 个可能的日期，并认为其中新疆维吾尔自治区喀什地区巴楚县是最有可能的地点，5 月 24 日或 7 月 20 日是最有可能的日期；由附件 3 同样得到 3 个可能的地点与 6 个可能的日期，认为湖北省十堰市郧西县与陕西省商洛市山阳县均是可能的地点，
3、可能的日期为 2 月 6 日或11 月 6 日前后。对于问题四，首先用 MATLAB 进行图像处理并得到等时间间隔的图片，然后经过筛选得到 21张图片。经滤镜处理后，由所得帧的图像得到影长与杆长的比例关系，进而得到不同时刻下的影长。在日期已知的情况下，问题四应用非线性拟合函数 fit 得到可行解，筛选后得到最可能地点为内蒙古自治区乌兰察布市丰镇市；若未给日期条件，在本题上一问的基础上，将太阳赤纬设为未知，利用 fit 函数求出可行解，经筛选得到最可能的地点为内蒙古自治区乌兰察布市，日期为 6 月 6 日或7 月 8 日，与准确日期相差无几。本文通过误差分析，证明本文所得结果具有很高的可信度。关
4、键词： 非线性曲线拟合 测量 经纬度2一、 问题重述如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。1. 建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用你们建立的模型画出 2015 年 10 月 22 日北京时间 9:00-15:00 之间天安门广场（北纬 39 度 54 分 26 秒,东经 116度 23 分 29 秒）3 米高的直杆的太阳影子长度的变化曲线。2. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。将你们的模型应用于附件
5、1 的影子顶点坐标数据，给出若干个可能的地点。3. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将你们的模型分别应用于附件 2 和附件 3 的影子顶点坐标数据，给出若干个可能的地点与日期。4. 附件 4 为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度为 2米。请建立确定视频拍摄地点的数学模型，并应用你们的模型给出若干个可能的拍摄地点。如果拍摄日期未知，你能否根据视频确定出拍摄地点与日期？二、问题分析1. 对于问题一，依据地理学知识，直杆的影长与直杆长度以及太阳高度角有关，而太阳高度角又与赤纬、直杆所处纬度、测量时刻等参数有关。
6、考虑到附件所给时间为北京时间，还应考虑经度带来的时差影响。基于以上考虑，构造出直杆的影长与赤纬、所处经度和纬度以及时刻之间的函数关系式，进而将题目所给参数代入关系式中，即可作出在给定时间天安门广场上给定直杆的影长随时间的变化曲线。2. 对于问题二，若想得到可能的地点，关键就在于得到直杆所处位置的经纬度坐标。由所给附件可以得到影长随时间的变化关系，因此可以代入问题一种建立的模型，通过最小二乘法迭代的方法估计出参数的可能取值，即直杆可能位置的经纬坐标。3. 对于问题三，仅是在问题二的基础上增多了日期这一未知量。这一未知量可以通过赤纬与日期对应表查询，因此同第二问一样，可以由附件得到影长随时间的变化
7、关系，利用最小二乘法迭代估计出经纬度、赤纬三个参数可能取值，从而给出可能的地点和日期。4. 对于问题四，通过 MATLAB 导入视频，并等时间间距地 21 张截图。在经过滤镜处理后，可以由截图得到影长与杆长的比例关系，从而得到不同时间下的影长。在日期已知的情况下，问题可以利用问题二的模型；若日期未知，则可以利用问题三的模型，得到所求参量。3三、模型假设1. 地球是半径为 的均匀圆球，即忽略其表面地形及地球微小的椭球率带来的影响；R2. 忽略太阳光在穿过地球大气层时产生的折射；3. 照射到地球的太阳光为平行光；4. 黄赤交角为一固定值。四、 符号约定0太阳直射点纬度测量点的纬度l测量点的经度a测
8、量点的时角1e太阳直射方向向量2地心指向测量点的方向向量h 太阳高度0L直杆长影长t测量点地方时0测量点地方时对应的北京时间五、模型建立与求解（一）问题一根据假设 1，以如图所示的方式建立空间直角坐标系 1。其中，以地心为坐标原点，垂直于赤道平面指向正北的轴为 z 轴，沿赤道平面指向太阳方向的轴为 x 轴，与二者成左手关系的轴为 y 轴。假定此时太阳直射点为地球上一点 A，设其纬度为 （北纬为正，南纬为负，以下涉及纬度符号均0同理） 。显然，由于 x 轴指向太阳方向，因此直射点 A 位于 oxz 平面内。设此时刻太阳直射方向为。设测量点纬度为 ，其时角 2为 ，其与地心构成的向量为 。1ea
9、2e4 0e1xzyA ae2赤 道观 测 点R图一 模型推导示意图根据以上设定，有100(cos,sin)aeRR2 ,a可求得 与 之间的空间 的余弦值为1e2（1）1200coscossins| ae根据太阳高度角 的定义，有 ，因此可得hh（2）00sincoscsina由图二所示的示意图可得，影长、直杆长之间存在关系 （3）0tanhL5h直 杆 长L0影 长 L阳 光图二 影长、杆长与太阳高度的关系在上述推导中，时角 ，其中 为测量点地方时。考虑到附录中所给时间均为北(12)att京时间，因此此处统一将地方时转化为北京时间表示。设测量点经度为 ，为方便起见，不妨假设l测量点在北京所
10、在的东八区的中间经线（120E）以东，因此二者之间经度差为 120lA由于地球每自转 24h 便转过 360，因此近似可认为经度上每度引起的时差为 4min，由前可知测量点在 120E 以东，因此测量点地方时与东八区区时有如下关系 0(12)5lt因此（4）0(12)5alt其中 为测量点地方时 对应的北京时间。0tt联立式（2） 、 （3） 、 （4） ，可得地球上某一点处直杆的影长与北京时间的函数关系式 （5）200 0000 0(12)cosssin15ci12ltLlt 查表 3可知 2015 年 10 月 22 日太阳直射点 ，并将已知条件 m，0.83 3L6， ，使用 MATLA
11、B 即可作出题目给定时间内天安门广场上 3m 长的39.07216.394lE直杆的影长变化曲线，如图三所示。图三 问题一所求曲线（二）问题二由附件可以得到测量点直杆的影长与北京时间的对应数据，可以在问题一所得到的模型的基础上，将待求的经、纬度坐标以及杆长看作为模型的参数，调用 MATLAB 中 lsqcurvefit4函数，通过迭代的方法实现最小二乘拟合，得到参数的估计值。Lsqcurvefit 函数是进行最小二乘非线性曲线拟合的函数，设定初始迭代起点进行拟合求最优解。为了更快得到最优解，将迭代起点选择为期望值附近（但并未限定所求解的范围） ，即杆长期望在 0 至 3 米之间，经度位于-18
12、0至 180之间，纬度的正弦值位于-1 到 1 之间。求解完成后，基于以下原则，在 MATLAB 返回的可行解当中进行了初步的筛选：1. 所有参数均应为纯实数，即虚部应等于零；2. 出于实际考虑，要求杆长均大于 0.5m；3. 参数完全重复的或相当接近的，取其中拟合程度最佳的一组参数。由附件 1 得到的数据见下表（时间已化为十进制）：北京时间（时） 影长（m ）14.7 1.14962582614.75 1.182198976714.8 1.21529695514.85 1.24905105214.9 1.2831953414.95 1.31799314915 1.35336404915.05
13、 1.38938709115.1 1.42615285615.15 1.46339985315.2 1.50148162215.25 1.54023181715.3 1.57985331615.35 1.62014451515.4 1.66127061315.45 1.70329063315.5 1.7462059115.55 1.79005091515.6 1.83501427215.65 1.88087500115.7 1.927918447表一 由附件 1 计算得到的北京时间与对应影长具体程序见附录，在代入由附件 1 得到的数据之后，直接得到的可行解如下：杆长（m） 纬度 经度 出现次数
14、1.9937 18.98 109.35 302.1479 -3.0152 104.32 2.2523 -3.6179 102.6662.4734 23.369 101.98 12.6254 24.185 99.892 12.794 24.819 97.757 2.7953 24.823 97.7434表二 附件 1 数据第一次迭代后的结果由表格可以看到某些解出现的次数不止一次，这是因为不同点起点得到相同的结果，说明该结果可能为最优解。由表格知（98N，109.35 E）在可行解中出现的频率极高，说明其可信度相对于其他解要更高；而在（23.36924.823N，101.9897.74E）这一范围
15、中，尽管每组解出现的次数不高，8但是可以发现这些数据分布相当密集，因此可以看做是同一个地点；（3.01523.6179S，104.32102.66E ）范围之间的解同理，但是由于其出现次数较少，显然其可信度不及另两组解。为了增强以上结果的可信度，现改变迭代起点进行第二次拟合。在这一次拟合中，将迭代起点设置偏离期望值。得到的结果如下表所示：杆长（m） 纬度 经度 出现次数1.9937 18.98 109.3494 11.5114 -2.6700 114.7980 1表三 附件 1 数据第二次迭代后的结果（18.98N，109.35E）这一组解再次出现，证明了这一组解得可信度确实相当高。将以上得到
16、的可能解汇总，并确定具体的地点如下：编号 杆长（m） 纬度 经度 具体地点1 1.9937 18.98 109.3494 海南省白沙黎族自治县2 2.7945 24.821 97.75 云南省德宏傣族景颇族自治州盈江县3 2.2523 -3.6179 102.66 印度尼西亚 Sinar Gunung Tebat Karai Kepahiang Regency 4 1.5114 -2.6700 114.7980 印度尼西亚 Asia Baru KuripanBarito Kuala Regency表四 问题二可能的地点分别对得到的四组解进行误差可视化和再次拟合，图像如下（从上到下依次为编号1、
17、2、3、4 的解）：9图四 问题二所有可能解的拟合及误差误差图显示拟合值与真实值之间的误差均为 10-3 级别，这再一次证明了解的可信度。这是基于本模型得到的四个可能的地点。考虑到迭代过程中出现次数悬殊较大，有理由相信海南省白沙黎族自治县（18.98N ，109.3494E ）是最有可能的地点，而其他三个地点可能性相对较小。（三）问题三第三问在第二问的基础上又增添了日期这个未知参数。由于日期与赤纬存在一定的对应关系，同第二问思路相同，我们将赤纬作为另一个进 lsqcurvefit 拟合，此时初步筛选原则除去问题二中提10出的 3 条外，再增添一条：赤纬的范围需位于 23.4333S23.433
18、3N 之间。1. 对附件 2 的分析由附件 2 得到的数据见下表（时间已化为十进制）：北京时间（时） 影长（m ）12.6833 1.24725620512.7333 1.2227945912.7833 1.19892148612.8333 1.17542896412.8833 1.15243957312.9333 1.1299174712.9833 1.1078354813.0333 1.08625420613.0833 1.06508107213.1333 1.04444626513.1833 1.02426412613.2333 1.00464031413.2833 0.98549090
19、813.3333 0.96679049413.3833 0.94858473513.4333 0.93092788113.4833 0.9137517513.5333 0.89710905113.5833 0.88097376213.6333 0.86549225913.6833 0.850504468表五 由附件 2 计算得到的北京时间与对应影长将附件 2 提取的数据代入模型，直接得到的可行解如下：杆长（m） 纬度 经度 赤纬 出现次数2.0008 39.893 79.744 20.761 102.3484 10.953 97.605 -11.902 13.8662 7.3544 97.72
20、7 -6.2274 1表五 附件 2 数据第一次迭代后的结果同样，为了增强解得可信度，使迭代起点偏离期望住进行了第二次拟合，直接得到的可行解如下：杆长（m） 纬度 经度 赤纬 出现次数2.0008 39.8926 79.7443 20.7606 111表六 附件 2 数据第二次迭代后的结果由此可见，第二次得到的这组解可信度很高。将表四中三组解依次编号为 1、2、3，以赤纬进行查表确定日期，最终得到下列可能的时间和地点：编号 纬度 经度 具体地点 日期1 39.893 79.744 新疆维吾尔自治区喀什地区巴楚县 5 月 24 日或 7 月 20 日2 10.953 97.605 缅甸海 2 月
21、 18 日或 10 月 25 日3 7.3544 97.727 缅甸海 3 月 5 日或 10 月 9 日表七 问题三之附件 2 可能的地点和日期分别进行误差可视化，得到的图像如下（从上到下依次为编号 1、2、3 的解）：12图五 附件 2 所有可能解的拟合及误差从误差图及拟合图都可以看出，第 1 组解的误差远远小于第 2、第 3 组，其误差在 10-4 级别；而另两组解的误差达到了 10-2 甚至 10-1 级别，可信度不高。从实际意义上也可以看到，位于陆地的可能性显然高于位于海洋的可能性。因此基于以上分析，地点位于缅甸海、日期为 2 月 18 日或 10月 25 日以及 3 月 5 日或
22、10 月 9 日的情况仅存在理论上可能；而地点位于新疆维吾尔自治区喀什地区巴楚县，日期为 5 月 24 日或 7 月 20 日的可能性非常高。2. 对附件 3 的分析与上述过程完全类似，由附件 2 得到的数据见下表（时间已化为十进制）：北京时间（时） 影长（m ）13.15 3.53314218413.2 3.54676802913.25 3.56179764313.3 3.57810071513.35 3.59575078313.4 3.6149342813.45 3.63542598313.5 3.65721827213.55 3.68054111513.6 3.70516783613.6
23、5 3.73127802513.7 3.75891791113.75 3.78808788813.8 3.81870101513.85 3.85080961913.9 3.8845852213.95 3.91991182814 3.95687599214.05 3.9955347914.1 4.0357508351314.15 4.077863059表八 由附件 3 计算得到的北京时间与对应影长用附件 3 提取的数据运用 lsqcurvefit 函数，直接得到的可行解如下：杆长（m） 纬度 经度 赤纬 出现次数3.0356 -32.849 110.24 15.959 102.3484 32.8
24、49 110.24 -15.959 163.8662 7.3544 97.727 -6.2274 13.0371 33.06 110.27 -15.731 13.0378 33.35 110.29 -15.43 14.5644 -18.866 96.681 20.331 15.5373 -16.421 103.8 16.387 1表九 附件 3 数据第一次迭代后的结果更改迭代起点为偏离期望值的初始值，进行第二次拟合，得到的结果如下：杆长（m） 纬度 经度 赤纬 出现次数3.0356 32.8488 110.2450 -15.9593 33.0391 33.5224 110.3023 -15.2
25、465 1表十 附件 3 数据第二次迭代后的结果可以看出估计结果的出现的次数相差很大，可以认为出现次数仅为一次的估计值可信度过低因此此处仅保留出现次数相对较高的几次结果。对照赤纬表得到对应的日期，通过谷歌地图得到具体的地点，结果如下：编号 纬度 经度 具体地点 日期1 -32.849 110.24 印度洋 5 月 6 日或 8 月 9 日2 32.849 110.24 湖北省十堰市郧西县2 月 6 日或 11 月 6日143 33.35 110.29 陕西省商洛市山阳县2 月 7 日或 11 月 5日表十一 问题三之附件 3 可能的地点和日期分别进行误差可视化和再次拟合，得到的图像如下（从上到
26、下依次为编号 1、2、3 的解）：图六 附件 3 所有可能解的拟合及误差由误差图可以看到，得到的这几组地点与日期均与原数据拟合的相当好，理论上都具有很大的可能性；但是考虑到实际情况，仍然认为在陆地上的可能性要大于在海洋上的可能性，因此可以认15为可以认为第 2 组、第 3 组所得到的地点与日期可能性更大一些。（四）问题四为了获取视频中直杆的影长信息，首先调用 MATLAB 中的 VideoReader，将附件 4 导入到MATLAB 中，并通过程序实现每隔 1000 帧获取一帧图片，随后从中筛选出 61 张直杆的影子相对清晰的图片，然后筛选出 21 张等时间间距（2 分钟）的图片作为提取影长的
27、信息来源。 。由于摄像机在进行拍摄时使用的是透视角度，其特点是在平面上相互平行的直线在这一视角中延长线将会相交，因此从截取的图片中直接获得的杆长与影长的比例关系可能已经不是其真实的比例关系。因此，在进行杆长测量之前，先使用了 PhotoShop 软件对图片进行滤镜处理，垂直透视参数设置为-38，镜头校准为 49.1%，使其由透视视角转换为平面视角，进而便于测量影长与杆长之间的比例关系。对于完全相同的一幅图，滤镜处理前后图片的对比如图七所示。滤 镜 处 理 后图七 滤镜处理前后效果对比图为了尽可能减小人为因素给模型输出结果带来的误差，我们选择在 PhotoShop 软件中，以软件自动建立的坐标来
28、计算影长与杆长之间的像素长度，得到两者之间的比例关系，而非直接在图片上以刻度尺量取，这样做可以最大程度避免人为测量带来的随机误差。设软件上直杆顶端坐标为，底端坐标为 ，测得某时刻影子顶端坐标为 ，影子的实际长度为 ，根据比1(,)xy2(,)xy(,)xyL例关系有 221()()xyL则有（6）221()()y实际测得的数据如下表所示：北京时间 影子实际长度8.9128 2.22568.9461 2.2009168.9794 2.17189.0127 2.14799.046 2.11809.0793 2.08899.1126 2.06809.1459 2.04039.1792 2.01359
29、.2125 1.99119.2458 1.95309.2791 1.94119.3124 1.90759.3457 1.88589.379 1.85979.4123 1.83519.4456 1.80989.4789 1.78599.5122 1.75759.5455 1.73379.5788 1.7150表十二 视频中提取出不同时刻的实际影长1. 在日期已知的前提下由视频可以得到拍摄日期为 7 月 13 日，可得此时赤纬 =21.9167N。此时问题四类似于问0题二。将问题四的已知参数代入问题二的模型。由于所求参数均具有实际意义，因此对它们进行一定的限制，利用非线性拟合函数 fit 进行拟合
30、，得到的可能取值如下：编号 纬度 经度 出现次数1 -6.5203 126.5685 132 40.3735 113.2984 73 -90 140 1表十三 问题四日期已知情况下得到的可能解根据误差图误差图以及数据拟合图（图八）可以看到第 1、3 组数据的拟合效果不好，误差绝对值非常大，可信度均不高，应该被舍去。17图八 问题四各可能解的误差图及再次拟合同样由误差图可以看到，第 2 组数据（40.3735N ，113.2984E）这组解产生的误差相对很小，因此这组解的可信度是很高的。经查，此坐标位于内蒙古自治区乌兰察布市丰镇市。2. 在日期未知的前提下由于拍摄日期位置，因此无法得知赤纬，所以
31、基于上一问的基础上，增加太阳赤纬的未知数，进行非线性曲线拟合，得到的结果如下：纬度 经度 赤纬 出现次数40.5311 112.8863 22.4886 40.6507 112.5483 22.95456-5.3519 126.3530 22.9545 1418 -5.7684 126.4333 22.5879-90 0.39 100 1表十四 问题四日期未知情况下得到的可能解作出误差图和数据拟合对比图如图九所示：图九 问题四（日期未知情况下）各可能解的误差图及再次拟合根据作出的误差图和数据拟合对比图，排除掉误差较大的一组解，保留其余两组解，并查的其对应的具体地点和日期如下：纬度 经度 赤纬
32、具体地点 日期-5.6845 126.4174 22.6624 班达海 6 月 6 日或 7 月 7 日40.5578 112.8126 22.5879 内蒙古自治区乌兰察布市丰镇市 6 月 6 日或 7 月 8 日表十五 问题四（日期未知）的可能地点和日期尽管这两组解的拟合程度均相当好，但是考虑到实际情况仍然是位于陆地上的可能性高于位于19海洋上的可能性，因此测试地点更有可能位于内蒙古自治区乌兰察布市丰镇市，日期为 6 月 6 日或7 月 8 日。这一结果与日期已知的情况下得到的参数高度吻合，证明了这一算法的合理性。六、模型评价基于非线性曲线拟合函数 lsqcurvefit 和 fit 进行
33、系数估算来实现求解,运行速度快，容错性强，可信度高。而且给定不同的初值可能出现相同的解，根据解出现的频率可以评价解的可信度。但是由于拟合的误差，解可能在某个微小范围内波动，而且可能出现附属解。参考文献1肖志勇，刘宇翔，一种新的纬度测量方法，大学物理，第 29 卷第 9 期：52，2010。2徐宝菜，应振华，地球概论教程，北京：高等教育出版社，1983。3 Table of the Declination of the Sun， http:/ 年 9 月 13 日。4刘浩，韩晶，MATLAB R2014a 完全自学一本通，北京：电子工业出版社，2015。附录本文使用软件为 Matlab R201
34、3b、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);ta
35、nh=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:3c0=i i i
36、;c=lsqcurvefit(fun2,c0,xdata,ydata)；20endfunction a b=fun2s()%问题 2 求解函数 迭代起点集中%a 杆长 所求纬度 太阳赤纬 所求经度%b 杆长 所求纬度的 sin 太阳赤纬 sin 所求经度x=xlsread(1.xlsx);xdata=x(:,1);ydata=x(:,2);a=;options = optimset(Display,off);for i=-1:0.01:1;c0=3+3*i i i*180;c=lsqcurvefit(fun2,c0,xdata,ydata,options);a=a;c;endb=a;a(:,2
37、)=asin(a(:,2);a(:,2)=a(:,2)*180/pi;c=find(a(:,1)0.5a=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(1.xlsx);xdata=x(:
38、,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;endb=a;a(:,2)=asin(a(:,2);a(:,2)=a(:,2)*180/pi;c=find(a(:,1)0.5a=a(c,:);b=b(c,:);c=find(imag(a(:,2)=0);21a=a(c,:);b=b(c,:);c=find(imag(a(:,3)=0);a=a(c,:);b=b(c,:);c d=sort(a
39、(:,1);a=a(d,:);b=b(d,:);function fun2w(a)%问题 2 绘图函数x=xlsread(1.xlsx);xdata2=x(:,1);ydata2=x(:,2);t=14.7:0.05:15.7;y1=fun2(a,t);subplot(2,1,1);plot(xdata2,ydata2,b*-)hold onplot(t,y1,r-)title(原始数据与拟合数据的对比);xlabel(T/h);ylabel(H/m);legend(原始数据, 拟合数据);y2=y1-ydata2;subplot(2,1,2);plot(t,y2,r.-)title(误差图
40、);xlabel(T/h);ylabel(H/m);legend(误差);function fun2ww(a)%问题 2 绘图函数x=xlsread(1.xlsx);xdata2=x(:,1);ydata2=x(:,2);t=9:0.01:16;y1=fun2(a,t);plot(xdata2,ydata2,b*-)hold onplot(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/
41、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 杆长 所求纬度 太阳赤纬 所求经度22%b 杆长 所求纬度的 sin 太阳赤纬 sin 所求经度x=xlsread(4.xlsx);xdata=x(:,1);ydata=x(:,2);a=;options = optimset(Display,off);for i=-1:0.01:1;c0= i i ;c=lsqcurvefit(fu
42、n4,c0,xdata,ydata,options);a=a;c;endb=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,:);c d=sort(a(:,1);a=a(d,:);b=b(d,:);function a b=fun32s()%问题 3 附件 2 求解函数 迭代起点集中%a 杆长 所求纬度 太阳赤纬 所求经度%b 杆长 所求纬度的 sin 太阳赤纬 sin 所求经度x=xlsread(2.xl
43、sx);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;endb=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.5a=a(c,:);b=b(c,:);c=find(imag(a(:,2)=0);a=a(
44、c,:);b=b(c,:);c=find(a(:,3)-23.4393a=a(c,:);b=b(c,:);c=find(imag(a(:,4)=0);a=a(c,:);23b=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(2.xlsx);xdata=x(:,1);ydata=x(:,2);a=;options = optimset(Display,off)
45、;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;endb=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.5a=a(c,:);b=b(c,:);c=find(imag(a(:,2)=0);a=a(c,:);b=b(c,:);c=find(a(:,3)-23.4393a=a(c,:);b=b(c,:);c=find(ima
46、g(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(2.xlsx);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 onplot(t,y1,r-)title(原始数据与拟合数据的对比);xlabel(T/h);24ylabel(H/m);legend(原始数据, 拟合数据
47、);y2=y1-ydata2;subplot(2,1,2);plot(t,y2,r.-)title(误差图 );xlabel(T/h);ylabel(H/m);legend(误差);function fun32ww(a)%问题 3 附件 2 绘图函数x=xlsread(2.xlsx);xdata2=x(:,1);ydata2=x(:,2);t=9:0.01:16;y1=fun3(a,t);plot(xdata2,ydata2,b*-)hold onplot(t,y1,r-)xlabel(T/h);ylabel(H/m);function a b=fun33s()%问题 3 附件 3 求解函数
48、迭代起点分散%a 杆长 所求纬度 太阳赤纬 所求经度%b 杆长 所求纬度的 sin 太阳赤纬 sin 所求经度x=xlsread(3.xlsx);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;endb=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.5a=a(c,:);b=b(c,:);c=find(imag(a(:
