1、1 基于最小二乘拟合和蒙特卡洛 法实现太阳影子定位 摘要 本文 基于 太阳高度角的 相关 原理， 建立影子长度与 直杆长度、纬度、积日和时角的函数关系， 使用 最小二乘拟合法、蒙特卡洛法、 分类筛选法、 误差修正法来实现太阳影子定位。 对 于问题一，建立 多元分析模型来分析影子长度 与 纬度、积日 、 地方时间之间的变化规律。 固定其中两个变量 ，分析第三个变量 与影子长度 的变化规律并画出天安门广场直杆影子长度的变化曲线。 对于问题二， 建立 最小二乘拟合参数模型 ，使用附件 1 的数据，得出 直杆的地点可能为 ： 北纬 19 度 16 分 、 东经 108 度 39 分（中国海南省 ） ，
2、南纬 3 度 39分、 东经 102 度 41 分（ 南 苏门答腊 ）。 对于问题三， 建立 蒙特卡洛拟合参数模型 ， 得出大量拟合参数，通过拟合优度检验和分类筛选 ， 得出 附件 2 中直杆的可能地点 与日期 为 ：北纬 34 度 57 分、东经 118 度 53 分（ 江苏连云港 ）、 12 月 20 日或 12 月 22 日 ； 北纬 37 度 14 分、东经 119 度 31 分（渤海）、 12 月 21 日 ； 北纬 33 度 14 分，东经 115 度 49 分（安徽阜阳） 、 12 月 20 日或 12 月 22 日。 附件 3 中可能的地点与日期为：北纬 46 度 25分、东经
3、 119 度 39 分（内蒙古自治区） 、 12 月 20 日或 12 月 22 日 ； 北纬 36 度 57 分、东经 118 度 26 分（山东临沂） 、 12 月 21 日 ； 北纬 33 度 14 分，东经 118 度 4 分（江苏宿迁）、 12 月 20 日或 12 月 22 日。 对于问题四 ， 首先 使用 相似原理 计算 太阳影子长度， 结合问题 二 模型 得出若干可能地点为： 北纬 39 度 31 分、 东经 109 度 04 分（内蒙古鄂尔多斯市） 。 然后 建立 误差 修正 模型 对结果进行修正 。得到修正后的若干可能地点为： 北纬 49 度 36 分、东经 104 度 2
4、2 分（ 蒙古色楞格 ） ； 北纬 45 度 36 分、东经 110度 （ 蒙古赛音山达 ） ； 北纬 48 度 21 分 、 东经 110 度 （ 蒙古鄂嫩） 。 若日期未知， 可应用问题三中 的蒙特卡洛拟合模型得出 若干可能地点与日期。 关键词： 太阳高度角 最小二乘拟合 蒙特卡洛法 拟合优度 误 差修正模型 2 一 、问题重述 如何确定视频的拍摄地点和拍摄 日期 是视频数据分析的重要方面，太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的 地点和日期的一种方法。 1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用你们建立的模型画出 2015年 1
5、0月 22日北京时间 9:00-15:00之间天安门广场（北纬 39 度 54 分 26 秒 ,东经 116 度 23 分 29 秒） 3 米高的直杆的太阳影 子长度的变化曲线。 2.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。将你们的模型应用于附件 1 的影子顶点坐标数据，给出若干个可能的地点。 3. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将你们的模型分别应用于附件 2 和附件 3 的影子顶点坐标数据，给出若干个可能的地点与日期。 4附件 4 为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直
6、杆的高度为 2 米。请建立确定视频拍摄地点的数学模型，并应用你们的模型给出若干个可能的拍摄 地点。 如果拍摄日期未知，你能否根据视频确定出拍摄地点与日期？ 二 、问题分析 针对问题一，建立 多元分析模型来分析影子长度 L 与纬度 、积日 V 、 地方时间 T 之间的变化规律。太阳高度角 与赤纬 角 、纬度 、时角 t（ t=|180 15 |T ）之 间 有 如 下 关 系 式 1 ； sin( ) sin( ) sin( ) + cos( ) cos( ) cos( ),t（ 22 ( 2 8 4 V )2 3 . 4 5 s i n365 ）。 由上可得 L 与 、 V 、 T 之间的函数
7、关系式为： 20 1 ( , , ) ( ( , , ) s in ( ) )( , , )L f T VL f T Vf T V ，固定其中两个变量可得第三个变量与 太阳影子长度 L 的函数关系式， 分析第三个变量与 影子长度 L 的变化规律并画出天安门广场直杆影子长度的变 化曲线。 3 针对问题二， 建立 最小二乘拟合参数模型 。附件 1 中的时间是北京时间而并非杆子所处的当地时间，设 北京时间为 0T ， 0T T-a （ a 为时差） 。 在关系式20 1 ( , , )( , , )L f T VL f T V中未知参数为 0La， ， 把 影子长度 L 与 北京时间 0T 导入 m
8、atlab 中， 用最小二乘法拟合 参数 ,运用 拟合 优 度 newR （ newR 越接近 1，拟合程度越好） 检验得 出拟合参数 0La， 。由时差 a 可算出测量地点的经度( 120 )4a， 再 结合 纬度 可得到直杆 的地点 。 运用附件 1 中的数据结合该模型可得出若干可能的地点。 针对问题三， 建立 蒙特卡洛随机拟合参数模型。由于变量增加， 拟合出来的参数不唯一，故 本 文 采取 蒙特卡洛随机拟合模型产生大量拟合参数， 根据参数分布情况 并结合 拟合优度 挑选出 符合实际情况 的参数。在 函数 关系式20 1 ( , , )( , , )L f T VL f T V中，地方时间
9、 0T T-a （ 0T 为北 京时间 ） ， 分析可知在关系式中未知参数为 0,a, L ,V ， 把 L 与 0T 导入 matlab 中，运用最小二乘法 结合 拟合优度检验 得 出大量 拟合参数 0,a, L ,V ， ， 根据参数分布情况 ，挑选有代表性的参数 0,a, L ,V ， 然后 由 时差 a 可算出测量地点的经度 ， 再 结合 纬度 可得到直杆 的地点 。 运用附件 2、 3 中的数据结合该模型可得出若干可能的地点。 针对问题四， 首先采用 相似原理 求出影子近 似长度 ，运用问题 2 中的最小二乘拟合模型求出可能的地点，然后建立 误差修正模型修正结果。 本问已知 杆长 L
10、为 2m,由相似原理，有00LlLl （ l 与 0l 分别为 影像中的 杆长和影长 ） ，使用 截屏 编辑工具 对附件 4 中的视频进行处理 ， 每隔 4 分钟，输出一次 l 和 0l ， 计算得到 近似 影 子 长 度 0L 。根据 视频中的 北京 时间 0T 可以得到一组太阳影子长度 0L 和 0T 的数据 ，使用 最小二乘拟合 模型，可求出拍摄该视频 若干 可能的 地点。 由于实际影子长度与测量出的影子长度之间存在一定误差 （存在偏角） ，此处使用修正模型来减小结果误差 。 同理， 若该视频没有日期， 可先 求 测量出的 太阳影子长度 0L 与4 北京 时间 0T 的数据 ，再使用 蒙
11、特卡洛拟合参数 模型求出拍摄视频 若干 可能的 地点与日期。 三、模型假设 1、假设杆子与地面垂直 2、假设地球是圆的 3、 不考虑太阳光线在穿过大气层时的折射、太阳的视角、高山阻挡、海拔高度等因素 的影响 ，照射到地球的太阳光可以看成是平行光线 四、符号说明 符号 含义 单位 L 太阳 影子长度 米 纬度 度 V 积日 天 a 时差 分钟 太阳高度角 度 0L 杆长 米 赤纬角 度 t 时角 度 T 地方时 时 五、模型的建立与求解 5.1 问题一模型建立与求解 -多元分析模型 5.1.1 模型建立 本文建立了多元分析模型来分析影子长度 L 与纬度 、积日 V 、 地方 时间 T之间的变化规
12、律。 查阅文献 【 1 】 知 太阳影子长 度 L 与 太阳高度角 有关，它们有如下关系 : 0tan( ) LL ,其中 0L 为杆长，即有 0tan( )LL ， 具体 如下图所示： 5 太阳高度角 与赤纬角 、纬度 、时角 t （ t=|180 15 |T ）之间有如下关系式 : sin( ) sin( ) sin( ) +cos( ) cos( ) cos( ),t 赤纬角的计算公式如下： 2 ( 2 8 4 V )2 3 . 4 5 s i n 365 , V 为积日，即日期在年内的顺序号， 例如， 1 月 1 日其积日为 1，平年 12月 31 日的积日为 365。将 与 纬度 代
13、入 太阳 高度角计算公式，推导出sin( ) ( , , )f T V ，可得影子长度 L 与 纬度 、积日 V 、 时间 T 之间的函数关系式为： 20 1 ( , , )( , , )L f T VL f T V，固定其中两个变量可得第三个变量与太阳影子长度 L 的函数关系式， 分析第三个变量与 L 的变化规律并画出天安门广场 3 米长 直杆影子长度 在 2015 年 10 月 22 日北京时间 9:00-15:00 之间 的变化曲线。 5.1.2 模型求解： 求出 L 与各个参数的函数关系式 如下：20 1 ( s i n ( ) s i n ( ) c o s ( ) c o s (
14、) c o s ( | 1 8 0 1 5 |) ) 2 ( 2 8 4 V ), 2 3 . 4 5 s i ns i n ( ) s i n ( ) c o s ( ) c o s ( ) c o s ( | 1 8 0 1 5 |) 3 6 5LTL T 固定其中两个参数，分析第三个参数与太阳影子长度之间的关系。 （ 1） 取定时间 T ，纬度 ，由常识可知回归线左右太阳影子的情况略有不同，此处以北 半球为例， 设杆长 0 3,L 取 12T ， =12 度 、 23 度 26 分 、 40 度来分析太阳影子长度 L 与积日 V 之间的变化规律 ,使用 matlab 软件 (程序见附录
15、) 0LL6 画出 L 与 V 的图像如下： 0 100 200 300 40001234567积日影子长度N 2 3 2 6 N 4 0 N 1 2 结果分析如下： 由图像可知，在北纬 40 度 和北纬 23 度 26 分 ，正午 12 点，随着积日 V 的增加，影子长度先变短后变长，在 170V 左右，北纬 40 度 区域太阳影子长度达到最小 约 1米，此时北纬 40度 地区处于夏季，影子长度最短，符合实际情况。 170V左右，太阳直射北回归线，在正午 12 点，太阳影子长度为 0。由常识可知在正午 12 点，因为北纬 12 度 处在南北回归线之间，故它的太阳影子长度随着积日 V的增加，在
16、正午 12 点会经历两次太阳影子长度为 0，上图很好的验证了这一点。 （ 2）取定时间 T ,积日 V ，使用 matlab 软件画出 L 与 的图像如下： 此处， T =12，即为正午 12 点，积日 V =80、 180，横坐标纬度已 换算成弧度。杆长 0L =3. 7 0 0 . 5 1 1 . 505101520253035404550北半球纬度影子长度积日 8 0积日 1 8 0结果分析如下： 由图像 （ 这里 纬度的单位已换算成弧度） 可知，积日 V=80 时，在正午 12点，随着纬度 的增加，太阳影子长度从 0 开始增加，当 达到 90时，太阳影子达到最长。由常识可知 V=80
17、左右，太阳直射赤道，随着纬度增加，太阳影子长度逐渐增加，上图符合实际情况。 积日 V=180 时，在正午 12 点，随着纬度 的增加，太阳影子长度先减小，后增加。在 V=180 左右，太阳直射北回归线，故当纬度从 0 开始增加时，太阳影子长度先减小，在北回归线处为 0，后增加，在 90处达到最大，上图符合实际情况。 （ 3）取定积日 V ,纬度 .此处 积日 V =295， 纬度为 北纬 39 度 54 分 26 秒 、东经 116 度 23 分 29 秒（天安门广场） ， 经过计算得到 太阳 影子长度 L 与时间 T 之间的函数 关系式如下： 23 1 ( 0 . 0 4 5 0 . 7 6
18、 8 8 c o s ( | 1 5 1 8 0 |) )0 . 0 4 5 0 . 7 6 8 8 c o s ( | 1 5 1 8 0 |)TL T ， 使用 matlab 软件 （程序见附录）画出 天安门广场 2015 年 10 月 22 日北京时间 9:00-15:00 之间 太阳影子长度 L 与 时间 T 的变化曲线如下：8 9 10 11 12 13 14 1522 . 533 . 544 . 5时间影子长度结果分析如下：由上图分析可知， 固定积日 V 和纬度 ,随着时间的增加，太阳影子长度先减小，在 12 点左右达到最小 2.2 米，随后随着时间的增加而增加，与实际情况相符。
19、5.2 问题二模型 -最小二乘拟合参数模型： 5.2.1 模型求解 由第一问，有如下函数关系式： sin( ) sin( ) sin( ) +cos( ) cos( ) cos( ),t 2 ( 2 8 4 V )2 3 . 4 5 s i n ,365 sin( ) ( , , )f T V ， 20 1 ( , , )( , , )L f T VL f T V 附件 1 中所给的时间是北京时间而并非杆子所处的当地时间， ，由测量日期可知积日 V ， 设 北 京 时 间 为 0T ， 且有 0T T-a 。 分 析 可 知 在 关 系 式20 1 ( , , )( , , )L f T VL
20、 f T V中未知参数为 0La， ， 将附件 1 中影子顶点坐标换算成影子长度 L ， 把 L 与 0T 导入 matlab 中，运用最小二乘法 求出拟合参数9 0La， ， 由时差 a 可算出测量地点的 经度 （ 120 4a） ， 由经纬度 ,可得到直杆 的地点 ，最后进行拟合优度分析。 其原理如下： 对非线性方程 ， 计算残差平方和 2()cQ L L （ CL 为 用拟合参数 计算出的影子长度 ） 以及 2RL 。 拟合 优 度 newR =1 QR， newR 越接近 1，说明拟合参数 0La， 越可靠， newR 越小，说明拟合参数越不可靠。根据 newR 的大小判断拟合出的参数
21、是否可靠。 5.2.2 模型求解： 使用 matlab 软件计算得到参数 0La， 如下（程序见附录）： 参数 a 0L 拟合值 1 N 1916 -45.4min 2.0365m 拟合值 2 S 339 -69.27min 2.2527m 拟合优度结果分析：使用 matlab（程序见附录）计算得到 拟合值 1、 2 的 拟合优度 1 0.9997newR 、 new 2R 0.9995 ， 12,new newRR都接近于 1，说明上述拟合值 1、 拟合值 2 可靠。 由此得到两个可能的地点如下表： 地点 经度 纬度 所属地区 地点 1 E: 10839 N: 1916 中国海南 省 地点
22、2 E: 10241 S: 339 南 苏门答腊 5.3 问题三 模型 -蒙特卡洛随机 拟合参数模型： 5.3.1 模型建立 由于 参数 增加， 拟合出的 参数不唯一，故 本问 建立 蒙特卡洛拟合模型产生大量拟合参数， 根据参数分布情况 选出拟合优度很高 且符合实际情况的参数 。 现由10 附件 2、 3 已知 直杆 太阳影 子顶点坐标和北京时间 0T ,可换算出太阳影子长度 L 和地方时间 T ( 0T T-a ,a 为时差） 。 分析可知在关系式 20 1 ( , , )( , , )L f T VL f T V中未知参数为 0,a, L ,V ， 把 L 与 0T 导入 matlab 中
23、，运用最小二乘法 得到 大量 拟合参数 0,a, L ,V ， 根据参数分布 规律 并结合 拟合优度 （ newR 越接近 1，说明拟合参数 0,a,L,V 越可靠， newR 越小，说明拟合参数越不可靠。 ） 挑选 符合实际情况 的参数 ，此时可由时差 a 算出测量地点的经度 ， 再结合 纬度 可得到直杆 的地点 ，由积日 V 可算出日期。 5.3.2 模型求解： 使用 matlab（程序见附录） 求出 大量 拟合参数 0,a, L , ,V 根据参数分布情况得到三类具有代表性的拟合参数及 相应的 拟合优度 newR 见下表 : 附件 2 杆长 0L 纬度 时差 a 积日 V 拟合优度 ne
24、wR 0.600888 0.649033 -1.98926 -10.3739 0.91827 0.659672 0.61037 -4.44488 -10.5632 0.919039 0.715595 0.577863 -10.0269 -11.2387 0.920741 附件 3 杆长 0L 纬度 时差 a 积日 V 拟合优度 newR 1.158541 0.807488 -1.39226 -9.07736 0.983326 1.969888 0.61298 -6.28522 -10.0661 0.988077 2.10947 0.580352 -5.31753 -10.9593 0.9872
25、18 通过计算得到附件 2 与附件 3 中直杆可能的地点与日期如下表： 11 附件 2 经度 纬度 地点 日期 E: 31119 N: 1437 渤海 12.21 E: 53118 N: 5734 江苏连云港 12.20 或 12.22 E: 49115 N: 1433 安徽阜阳 12.20 或 12.22 附件 3 经度 纬度 地点 日期 E: 11939 N: 4625 内蒙古自治区 12.20 或 12.22 E: 11826 N: 3657 山东临沂 12.21 E:11804 N: 3314 江苏宿迁 12.20 或 12.22 5.4 问题四 模型的建立与求解 -相似原理及误差 修
26、正模型 5.4.1 相似原理 求解影子长度 并进行参数拟合 （ 1）问题分析 首先假设 视频中的 时间是北京时间 ， 附件 4 是 2m 一根直杆在太阳下的影子变化的视频 ，由相似原理，有00LlLl ， l 与 0l 分别为直杆和影子 的像素，使用截屏 编辑工具， 每隔 4 分钟，输出一次影像中直杆和 太阳影子 的像素比例，计算得到 影长 0L ，由视频中的 北京 时间可以得到一组太阳影子长度 0L 随着 北京 时间变化的数据 ， 使用最小二乘拟合参数 模型，求出拍摄该视频可能的 地点。 但是，上述方法可能存在较大误差，究其原因，是因为 并没有考虑到 相机是否正对直杆和影子所成平面 ， 现
27、构建误差修正模型来修正测量出的影子长度 可能 存在的误差，从而修正拍摄地点。 （ 2）问题求解 12 把求出的影子长度导入 matlab（ 程序 见附录）求解出可能的地点为： 纬度 经度 地点 北纬 39 度 31 分 东经 109 度 04 分 内蒙古巴彦淖尔市 5.4.2 误差修正模型的建立与求解： （ 1）模型建立： 假设由视频测量出的各时间段影子长度 CL 与实际影子长度 L 相差一个 修正常数 b ，即有 CL b L 20 1 ( , , ) ,( , , )L f T Vf T V得到 20 1 ( , , )( , , )C L f T VLbf T V，将相似原理模型中的 C
28、L 与北京时间 0T 导入 matlab 中，使用最小二乘法拟合出未知参数 , , ,ab 求得相对精确的若干可能地点 。 若该视频没有日期， 分析上述关系式 20 1 ( , , )( , , )C L f T VLbf T V，未知参数为 , , ,ab V ,在 matlab 中 拟合出上述未知参数，求出拍摄地点与日期。 （ 2） 模型求解： 在 matlab 中求解出参数 , , ,ab 计算得到若干修正后的地点如下表： 纬度 经度 地点 北纬 49 度 36 分 东经 104 度 22 分 蒙古色楞格 北纬 45 度 36 分 东经 110 度 蒙古赛音山达 北纬 48 度 21 分
29、 东经 110 度 蒙古鄂嫩 5.4.3 结果分析 因为拍摄时相机不一定正对直杆与影子所成平面 ，所以建立了上述修正模型对地点进行修正 ，假设 中 由视频测量出的各时间段影子长度 CL 与实际影子长度L 相差一个修正 常数 b ，但实际上 b 是随着时间的变化而变化的，故上述 误差 修正模型只能修正 部分误差 。 5.4.4 思考题 若拍摄日期未知，同样可 使 用蒙特卡洛拟合参数模型拟合出未知参数，但此时拟合的未知参数有 4 个，会存在较大的误差， 此时应该 做 大量的 参数 拟合，得13 出 拟合优度高且 符合实际的参数， 然后 分类筛选出具有代表性的参数，确定若干可能的地点。 六、模型优缺
30、点 优点： 问题 1 建立了多元分析模型，画出影子长度与各个因素之间的图像，简洁明 了的刻画了太阳影子长度与各个因素之间的变化规律。 问题 2、 3 使用拟合优度 分析出 拟合出的参数的合理性。 问题 4 建立了误差修正模型来修正测量出的太阳影子长度存在的误差。 缺点： 当拟合的参数 较多时，最小二乘非线性拟合模型拟合出的参数会存在较大的误差，问题 4 中 误差修正模型 拟合 所得结果并不太好，究其原因 可能 是由于 测量出的太阳影子数据误差较大。 七、参考文献 1 陈晓勇 郑科科 对建筑日照计算中太阳赤纬角公式的探讨 J 浙江建筑 2011 年 9 月 2 百度百科 赤纬角 http:/ 2
31、015/9/14 3司守奎 孙玺菁 数学建 模算法与应用 国防工业出版社 2011 年 93P 4王昌明 可照时数和太阳高度角计算公式的简化证明 J 山东气象 1989（ 02） 5 汪和平 影子与季节、纬度的关系 J 数学园地 2008 年第 9 期 6 刘彦刚 四季与赤纬角和时角 J 太阳能 1991 年 03 期 7唐家德 基于 MATLAB 的非线性曲线拟合 J计算机与现代化 2008 年第 6 期 八、附录 问题一代码： （ 1） 取定时间和积日，分析纬度的与影子长度的变化规律图 clc,clear ; x=0:0.1:pi/2; v=80;%赤道与南回归线之间。 t1=sin(23
32、.45*sin(2*pi*(284+v)./365)*pi/180); t2=sqrt(1-t1.2); 14 y=3*sqrt(1-(sin(x)*t1+cos(x)*t2).2)./(sin(x)*t1+cos(x)*t2; plot(x,y,r) hold on x=0:0.1:pi/2; v=180;%赤道与北回归线之间。 t1=sin(23.45*sin(2*pi*(284+v)./365)*pi/180); t2=sqrt(1-t1.2); y=3*sqrt(1-(sin(x)*t1+cos(x)*t2).2)./(sin(x)*t1+cos(x)*t2; plot(x,y) （
33、2）取定时间和纬度，分析积日的与影子长度的变化规律图 clc;clear; v=1:365;%取定纬度为北回归线，时间为 12时。 t1=sin(23.45*sin(2*pi*(284+v)/365)*pi/180); t2=sqrt(1-t1.2); y=3*sqrt(1-(sin(23.26/180)*pi)*t1+cos(23.26/180)*pi)*t2).2)./(sin(23.26/180)*pi)*t1+cos(23.26/180)*pi)*t2); plot(v,y,r) hold on v=1:365;%北回归线以外任一点的纬度，时间为 12时。 t1=sin(23.45*s
34、in(2*pi*(284+v)/365)*pi/180); t2=sqrt(1-t1.2); y=3*sqrt(1-(sin(40/180)*pi)*t1+cos(40/180)*pi)*t2).2)./(sin(40/180)*pi)*t1+cos(40/180)*pi)*t2); plot(v,y,k) hold on v=1:365;%取赤道与北回归线之间任一点的纬度，时间为 12时。 t1=sin(23.45*sin(2*pi*(284+v)/365)*pi/180); t2=sqrt(1-t1.2); y=3*sqrt(1-(sin(12/180)*pi)*t1+cos(12/180
35、)*pi)*t2).2)./(sin(12/1815 0)*pi)*t1+cos(12/180)*pi)*t2); plot(v,y) （ 3） 北京一竿子影子随时间的变化规律图） clc;clear; x=9:0.1:15; t=0.045+0.7688*cos(15*x-180)*pi/180) y=(3*sqrt(1-t.2)./t; plot(x,y) 问题二 拟合并优度检验 代码： clc,clear a=; h=; for p=1:500%下面一段程序做数据拟合和优度检验（采用蒙特卡洛方法，产生500 组拟合值）。 f=(canshu,xdata)canshu(1)*sqrt(1-
36、(sin(canshu(2)*sin(10.511*pi/180)+. cos(canshu(2)*cos(10.511*pi/180)*cos(180-0.25*(xdata+canshu(3)*pi/ . 180).2)./(sin(canshu(2)*sin(10.511*pi/180)+cos(canshu(2)*cos . (10.511*pi/180)*cos(180-0.25*(xdata+canshu(3)*pi/180); y0=xlsread(附件 1 数据 .xlsx,Sheet1,B2:B22); x0=xlsread(附件 1 数据 .xlsx,Sheet1,A2:A
37、22); canshu=lsqcurvefit(f,rand(3,1),x0,y0); a=canshu;%下面一段程序做优度检验。 L=a(1);%杆长。 b=a(2);%纬度。 t=a(3);%时间差，用于确定经度。 for j=1:3 h(p,j)=a(j); 16 end xdata=xlsread(附件 1 数据 .xlsx,Sheet1,A2:A22); y=L*sqrt(1-(sin(b)*sin(10.511*pi/180)+cos(b)*cos(10.511*pi/180)*cos .(180-0.25*(xdata+t)*pi/180).2)./(sin(b)*sin(10
38、.511*pi/180)+ . cos(b)*cos(10.511*pi/180)*cos(180-0.25*(xdata+t)*pi/180); plot(xdata,y,*);%拟合值的散点图。 y0=xlsread(附件 1 数据 .xlsx,Sheet1,B2:B22);%实测值。 k=; d=; for i=1:21 k(i)=(y(i)-y0(i)2; d(i)=y0(i)2; end k; sum(k); sum(d); Renew=1-sqrt(sum(k)/sum(d);%得出优度值 Renew. h(p,4)=Renew; xlswrite(输出 1.xlsx,h);%把大
39、量的拟合值导入 Excel 表，进行分类，得出规律并挑选。 End 问题二 拟合并优度检验 代码 ： clc,clear h=; a=; for i=1:500%下面一段程序做数据拟合和优度检验（采用蒙特卡洛方法，产生500 组拟合值）。 17 f=(canshu,xdata)canshu(1)*sqrt(1-(sin(canshu(2)*sin(23.45*pi*sin .(2*pi*(284+canshu(4)/365)/180)+cos(canshu(2)*cos(23.45*pi*sin . (2*pi*(284+canshu(4)/365)/180)*cos(180-0.25*(xd
40、ata+canshu(3) . *pi/180).2)./(sin(canshu(2)*sin(23.45*pi*sin(2*pi*(284+canshu(4) . )/365)/180)+cos(canshu(2)*cos(23.45*pi*sin(2*pi* . (284+canshu(4)/365)/180)*cos(180-0.25*(xdata+canshu(3)*pi/180); y0=xlsread(附件 3 数据 .xlsx,Sheet1,B2:B12); x0=xlsread(附件 3 数据 .xlsx,Sheet1,A2:A12); canshu=lsqcurvefit(f
41、,rand(4,1),x0,y0); a=canshu; for j=1:4%下面一段程序做优度检验。 h(i,j)=a(j); end L=a(1);%杆长。 b=a(2);%纬度。 t=a(3);%时间差，用于确定经度。 T=a(4);%积日。 xdata=xlsread(附件 3 数据 .xlsx,Sheet1,A2:A12); y=L*sqrt(1-(sin(b)*sin(23.45*pi*sin(2*pi*(284+T)/365)/180)+cos(b) .*cos(23.45*pi*sin(2*pi*(284-T)/365)/180)*cos(180-0.25*(xdata+t)
42、.18 *pi/180).2)./(sin(b)*sin(23.45*pi*sin(2*pi*(284+T)/365)/180)+ . cos(b)*cos(23.45*pi*sin(2*pi*(284+T)/365)/180)*cos(180-0.25* . (xdata+t)*pi/180); plot(xdata,y,*);%拟合值的散点图。 y0=xlsread(附件 3 数据 .xlsx,Sheet1,B2:B12);%实测值。 k=; d=; for p=1:11 k(p)=(y(p)-y0(p)2; d(p)=y0(p)2; end sum(k); sum(d); Renew=1
43、-sqrt(sum(k)/sum(d);%得出优度值 Renew. h(i,5)=Renew; xlswrite(输出 3.xlsx,h);%把大量的拟合值导入 Excel 表，进行分类，得出规律并挑选。 End 问题四拟合并优度检验代码： clc,clear for p=1:500%下面一段程序做数据拟合和优度检验（采用蒙特卡洛方法，产生500 组拟合值）。 f=(canshu,xdata)2*sqrt(1-(sin(canshu(1)*sin(23.21*pi/180)+cos . (canshu(1)*cos(23.21*pi/180)*cos(180-0.25*(xdata+cansh
44、u(2) . *pi/180).2)./(sin(canshu(1)*sin(23.21*pi/180)+cos(canshu(1) . 19 *cos(23.21*pi/180)*cos(180-0.25*(xdata+canshu(2)*pi/180); y0=xlsread(附件 4 数据 .xlsx,Sheet1,B2:B12); x0=xlsread(附件 4 数据 .xlsx,Sheet1,A2:A12); canshu=lsqcurvefit(f,rand(2,1),x0,y0); a=canshu;%下面一段程序做优度检验。 L=a(1); b=a(2); for j=1:2
45、h(p,j)=a(j); end xdata=xlsread(附件 4 数据 .xlsx,Sheet1,A2:A12); y=2*sqrt(1-(sin(L)*sin(23.21*pi/180)+cos(L)*cos(23.21*pi/180)* . cos(180-0.25*(xdata+b)*pi/180).2)./(sin(L)*sin(23.21*pi/ . 180)+cos(L)*cos(23.21*pi/180)*cos(180-0.25*(xdata+b)*pi/180);%拟合值的散点图。 y0=xlsread(附件 4.1 数据 .xlsx,Sheet1,B2:B12);%实
46、测值。 k=; d=; for i=1:11 k(i)=(y(i)-y0(i)2; d(i)=y0(i)2; end k; sum(k); sum(d); Renew=1-sqrt(sum(k)/sum(d);%得出优度值 Renew. h(p,3)=Renew; xlswrite(输出 4.xlsx,h);%把大量的拟合值导入 Excel 表，进行分类，得出规律并挑选。 20 End 问题四误差修正代码： clc,clear for p=1:500%下面一段程序做数据拟合和优度检验（采用蒙特卡洛方法，产生500 组拟合值）。 f=(canshu,xdata)2*sqrt(1-(sin(can
47、shu(1)*sin(23.21*pi/180)+cos . (canshu(1)*cos(23.21*pi/180)*cos(180-0.25*(xdata+canshu(2) . *pi/180).2)./(sin(canshu(1)*sin(23.21*pi/180)+cos(canshu(1) . *cos(23.21*pi/180)*cos(180-0.25*(xdata+canshu(2)*pi/180)+canshu(3); y0=xlsread(附件 4 数据 .xlsx,Sheet1,B2:B12); x0=xlsread(附件 4 数据 .xlsx,Sheet1,A2:A12); canshu=lsqcurvefit(f,rand(3,1),x0,y0); a=canshu;%下面一段程序做优度检验。 L=a(1);%纬度。 b=a(2);%时差，用来求经度。 M=a(3);%修正系数。 for j=1:3 h(p,j)=a(j); end xdata=xlsread(附件 4 数据 .xl
