2015年高教社杯数学建模论文.pdfVIP

1 基于最小二乘拟合和蒙特卡洛基于最小二乘拟合和蒙特卡洛法实现太阳影子定位法实现太阳影子定位 摘要摘要 本文基于太阳高度角的相关原理，建立影子长度与直杆长度、纬度、积日和 时角的函数关系，使用最小二乘拟合法、蒙特卡洛法、分类筛选法、误差修正法 来实现太阳影子定位。 对于问题一，建立多元分析模型来分析影子长度与纬度、积日、地方时间之 间的变化规律。 固定其中两个变量，分析第三个变量与影子长度的变化规律并画 出天安门广场直杆影子长度的变化曲线。 对于问题二，建立最小二乘拟合参数模型，使用附件 1 的数据，得出直杆的 地点可能为：北纬 19 度 16 分、东经 108 度 39 分（中国海南省） ，南纬 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 分、东经 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 度 22 分（蒙古色楞格） ；北纬 45 度 36 分、东经 110 度（蒙古赛音山达） ；北纬 48 度 21 分、东经 110 度（蒙古鄂嫩） 。 若日期未知，可应用问题三中的蒙特卡洛拟合模型得出若干可能地点与日 期。 关键词：关键词：太阳高度角 最小二乘拟合 蒙特卡洛法 拟合优度 误差修正模型 2 一一、问题重述、问题重述 如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面， 太阳影子 定位技术就是通过分析视频中物体的太阳影子变化， 确定视频拍摄的地点和日期 的一种方法。 1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律， 并应用你们建立的模型画出 2015年10月22 日北京时间 9:00-15:00 之间天安门 广场（北纬 39 度 54 分 26 秒,东经 116 度 23 分 29 秒）3 米高的直杆的太阳影子 长度的变化曲线。 2.根据某固定直杆在水平地面上的太阳影子顶点坐标数据， 建立数学模型确 定直杆所处的地点。将你们的模型应用于附件 1 的影子顶点坐标数据，给出若干 个可能的地点。 3. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型 确定直杆所处的地点和日期。 将你们的模型分别应用于附件 2 和附件 3 的影子顶 点坐标数据，给出若干个可能的地点与日期。 4附件 4 为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估 计出直杆的高度为 2 米。请建立确定视频拍摄地点的数学模型，并应用你们的模 型给出若干个可能的拍摄地点。 如果拍摄日期未知，你能否根据视频确定出拍摄地点与日期？ 二二、问题分析、问题分析 针对问题一，建立多元分析模型来分析影子长度L与纬度、积日V、地方 时间T之间的变化规律。 太阳高度角与赤纬角、 纬度、 时角t（t=|180 15 |T） 之 间 有 如 下 关 系 式 11 ； sin( )sin( )sin( )+cos( )cos( )cos( ),t （ 2 2 (284V) 23.45sin 365 ） 。由上可得L与、V、T之间的函数关系式 为： 2 0 1( , , ) ( ( , , )sin( ) ( , , ) LfT V LfT V fT V ，固定其中两个变量可得第三个变 量与太阳影子长度L的函数关系式， 分析第三个变量与影子长度L的变化规律并 画出天安门广场直杆影子长度的变化曲线。 3 针对问题二，建立最小二乘拟合参数模型。附件 1 中的时间是北京时间而并 非杆子所处的当地时间，设北京时间为 0 T， 0 TT -a（a为时差） 。在关系式 2 0 1( , , ) ( , , ) LfT V L fT V 中未知参数为 0 La，把影子长度L与北京时间 0 T导 入 matlab 中，用最小二乘法拟合参数,运用拟合优度 new R（ new R越接近 1，拟合 程度越好）检验得出拟合参数 0 La，。由时差a可算出测量地点的经度 (120) 4 a ，再结合纬度可得到直杆的地点。运用附件 1 中的数据结合该 模型可得出若干可能的地点。 针对问题三，建立蒙特卡洛随机拟合参数模型。由于变量增加，拟合出来的 参数不唯一，故本文采取蒙特卡洛随机拟合模型产生大量拟合参数，根据参数分 布 情 况并结合 拟合优度 挑选出符合实际情况 的参数。在 函数 关系式 2 0 1(,) (,) LfT V L fT V 中，地方时间 0 TT -a（ 0 T为北京时间） ，分析可知在关系 式中未知参数为 0 ,a, L ,V，把L与 0 T导入 matlab 中，运用最小二乘法结合拟 合优度检验得出大量拟合参数 0 ,a, L ,V， ，根据参数分布情况，挑选有代表性的 参数 0 ,a, L ,V，然后由时差 a 可算出测量地点的经度，再结合纬度可得到 直杆的地点。运用附件 2、3 中的数据结合该模型可得出若干可能的地点。 针对问题四，首先采用相似原理求出影子近似长度，运用问题 2 中的最小二 乘拟合模型求出可能的地点，然后建立误差修正模型修正结果。本问已知杆长L 为 2m,由相似原理，有 00 Ll Ll （l与 0 l分别为影像中的杆长和影长） ，使用截屏编 辑工具对附件 4 中的视频进行处理，每隔 4 分钟，输出一次l和 0 l，计算得到近 似影子长度 0 L。根据视频中的北京时间 0 T可以得到一组太阳影子长度 0 L和 0 T的 数据，使用最小二乘拟合模型，可求出拍摄该视频若干可能的地点。由于实际影 子长度与测量出的影子长度之间存在一定误差（存在偏角） ，此处使用修正模型 来减小结果误差。同理，若该视频没有日期，可先求测量出的太阳影子长度 0 L与 4 北京时间 0 T的数据， 再使用蒙特卡洛拟合参数模型求出拍摄视频若干可能的地点 与日期。 三、模型假设三、模型假设 1、假设杆子与地面垂直 2、假设地球是圆的 3、不考虑太阳光线在穿过大气层时的折射、太阳的视角、高山阻挡、海拔高度 等因素的影响，照射到地球的太阳光可以看成是平行光线 四、符号说明四、符号说明 符号 含义 单位 L 太阳影子长度 米 纬度 度 V 积日 天 a 时差 分钟 太阳高度角 度 0 L 杆长 米 赤纬角 度 t 时角 度 T 地方时 时 五、模型的建立与求解五、模型的建立与求解 5.15.1 问题一模型建立与求解问题一模型建立与求解- -多元分析模型多元分析模型 5.1.1 模型建立 本文建立了多元分析模型来分析影子长度L与纬度、积日V、地方时间T 之间的变化规律。 查阅文献 【 1 】知太阳影子长度 L与太阳高度角有关，它们有如下关 系: 0 tan( ) L L ,其中 0 L为杆长，即有 0 tan( ) L L ，具体如下图所示： 5 太阳高度角与赤纬角、纬度、时角t（t=|180 15 |T）之间有如下关系式: sin( )sin( )sin( )+cos( )cos( )cos( ),t 赤纬角的计算公式如下： 2 (284V) 23.45sin 365 , V为积日，即日期在年内的顺序号，例如，1 月 1 日其积日为 1，平年 12 月 31 日的积日为 365。将与纬度代入太阳高度角计算公式，推导出 sin( )( , , )fT V，可得影子长度L与纬度、积日V、时间T之间的函数关系 式为： 2 0 1( , , ) ( , , ) LfT V L fT V ，固定其中两个变量可得第三个变量与太阳影子长 度L的函数关系式， 分析第三个变量与L的变化规律并画出天安门广场 3 米长直 杆影子长度在 2015 年 10 月 22 日北京时间 9:00-15:00 之间的变化曲线。 5.1.25.1.2 模型求解：模型求解： 求出L与各个参数的函数关系式如下： 2 0 1 (sin( )sin( )cos( )cos( )cos(|180 15 |)2 (284V) ,23.45sin sin( )sin( )cos( )cos( )cos(|180 15 |)365 LT L T 固定其中两个参数，分析第三个参数与太阳影子长度之间的关系。 （1）取定时间T，纬度，由常识可知回归线左右太阳影子的情况略有不 同，此处以北半球为例，设杆长 0 3,L 取12T ，=12 度、23 度 26 分、40 度 来分析太阳影子长度L与积日V之间的变化规律,使用 matlab 软件(程序见附录) 0 L L 6 画出L与V的图像如下： 0100200300400 0 1 2 3 4 5 6 7 积 日 影子长度 N 23 26 N 40 N 12 结果分析如下： 由图像可知，在北纬 40 度和北纬 23 度 26 分，正午 12 点，随着积日V的增 加，影子长度先变短后变长，在170V 左右，北纬 40 度区域太阳影子长度达到 最小约 1 米， 此时北纬 40度地区处于夏季， 影子长度最短， 符合实际情况。170V 左右，太阳直射北回归线，在正午 12 点，太阳影子长度为 0。由常识可知在正 午 12 点，因为北纬 12 度处在南北回归线之间，故它的太阳影子长度随着积日 V 的增加，在正午 12 点会经历两次太阳影子长度为 0，上图很好的验证了这一点。 （2）取定时间T,积日V，使用 matlab 软件画出L与的图像如下： 此处，T=12，即为正午 12 点，积日V=80、180，横坐标纬度已换算成弧度。 杆长 0 L=3. 7 00.511.5 0 5 10 15 20 25 30 35 40 45 50 北 半 球 纬 度 影子长度 积 日 80 积 日 180 结果分析如下： 由图像（这里纬度的单位已换算成弧度）可知，积日 V=80 时，在正午 12 点，随着纬度的增加，太阳影子长度从 0 开始增加，当达到 90时，太阳 影子达到最长。由常识可知 V=80 左右，太阳直射赤道，随着纬度增加，太阳影 子长度逐渐增加，上图符合实际情况。 积日 V=180 时，在正午 12 点，随着纬度的增加，太阳影子长度先减小， 后增加。在 V=180 左右，太阳直射北回归线，故当纬度从 0 开始增加时，太阳影 子长度先减小，在北回归线处为 0，后增加，在 90处达到最大，上图符合实际 情况。 （3）取定积日V,纬度.此处积日V=295，纬度为北纬 39 度 54 分 26 秒、东 经 116 度 23 分 29 秒（天安门广场） ，经过计算得到太阳影子长度L与时间T之间的 函数关系式如下： 2 3 1 (0.0450.7688cos(|15180|) 0.0450.7688cos(|15180|) T L T ， 使用 matlab 软件（程序见附录）画出天安门广场2015 年 10 月 22 日北京时 间 9:00-15:00 之间太阳影子长度L与时间T的变化曲线如下： 8 9101112131415 2 2.5 3 3.5 4 4.5 时 间 影子长度 结果分析如下：由上图分析可知，固定积日V和纬度,随着时间的增加，太阳 影子长度先减小，在 12 点左右达到最小 2.2 米，随后随着时间的增加而增加， 与实际情况相符。 5.25.2 问题二模型问题二模型- -最小二乘拟合参数模型：最小二乘拟合参数模型： 5.2.1 模型求解 由第一问，有如下函数关系式：sin( )sin( )sin( )+cos( )cos( )cos( ),t 2 (284V) 23.45sin, 365 sin( )( , , )fT V ， 2 0 1( , , ) ( , , ) LfT V L fT V 附件 1 中所给的时间是北京时间而并非杆子所处的当地时间， ，由测量日期 可 知 积 日V ，设 北 京 时 间 为0 T， 且 有 0 TT -a。 分 析 可 知 在 关 系 式 2 0 1( , , ) ( , , ) LfT V L fT V 中未知参数为 0 La，将附件 1 中影子顶点坐标换算成 影子长度L，把L与 0 T导入 matlab 中，运用最小二乘法求出拟合参数 9 0 La，由时差 a 可算出测量地点的经度（120 4 a ） ，由经纬度, 可 得到直杆的地点，最后进行拟合优度分析。其原理如下： 对非线性方程，计算残差平方和 2 () c QLL （ C L为用拟合参数计算出的 影子长度）以及 2 RL。拟合优度 new R=1 Q R ， new R越接近 1，说明拟合参 数 0 La，越可靠， new R越小，说明拟合参数越不可靠。根据 new R的大小判断拟 合出的参数是否可靠。 5.2.2 模型求解： 使用 matlab 软件计算得到参数 0 La，如下（程序见附录） ： 参数 a 0 L 拟合值 1 N 1916 -45.4min 2.0365m 拟合值 2 S 3 39 -69.27min 2.2527m 拟合优度结果分析：使用 matlab（程序见附录）计算得到拟合值 1、2 的拟合优 度 1 0.9997 new R、 new2 R0.9995， 12 , newnew RR都接近于 1，说明上述拟合值 1、拟 合值 2 可靠。 由此得到两个可能的地点如下表： 地点 经度 纬度 所属地区 地点 1 E: 108 39 N: 1916 中国海南省 地点 2 E: 102 41 S: 3 39 南苏门答腊 5.35.3 问题三问题三模型模型- -蒙特卡洛随机蒙特卡洛随机拟合参数模型：拟合参数模型： 5.3.1 模型建立 由于参数增加，拟合出的参数不唯一，故本问建立蒙特卡洛拟合模型产生大 量拟合参数， 根据参数分布情况选出拟合优度很高且符合实际情况的参数。 现由 10 附件 2、 3 已知直杆太阳影子顶点坐标和北京时间 0 T,可换算出太阳影子长度L和 地方时间T( 0 TT -a,a为时差） 。分析可知在关系式 2 0 1( , , ) ( , , ) LfT V L fT V 中未 知参数为 0 ,a, L ,V，把L与 0 T导入 matlab 中，运用最小二乘法得到大量拟合 参数 0 ,a, L ,V，根据参数分布规律并结合拟合优度（ new R越接近 1，说明拟合参 数 0 ,a,L ,V越可靠， new R越小，说明拟合参数越不可靠。）挑选符合实际情况的 参数， 此时可由时差a算出测量地点的经度， 再结合纬度可得到直杆的地点， 由积日V可算出日期。 5.3.2 模型求解： 使用 matlab（程序见附录）求出大量拟合参数 0 ,a, L , ,V根据参数分布情况 得到三类具有代表性的拟合参数及相应的拟合优度 new R见下表: 附件 2 杆长 0 L 纬度 时差a 积日V 拟合优度 new R 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 杆长 0 L 纬度 时差a 积日V 拟合优度 new R 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.987218 通过计算得到附件 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:119 39 N:46 25 内蒙古自治区 12.20 或 12.22 E:118 26 N:36 57 山东临沂 12.21 E:118 04 N:3314 江苏宿迁 12.20 或 12.22 5.45.4 问题四问题四模型的建立与求解模型的建立与求解- -相似原理及误差相似原理及误差修正模型修正模型 5.4.1 相似原理求解影子长度并进行参数拟合 （1）问题分析 首先假设视频中的时间是北京时间，附件 4 是 2m 一根直杆在太阳下的影子 变化的视频，由相似原理，有 00 Ll Ll ，l与 0 l分别为直杆和影子的像素，使用截 屏编辑工具，每隔 4 分钟，输出一次影像中直杆和太阳影子的像素比例，计算得 到影长 0 L，由视频中的北京时间可以得到一组太阳影子长度 0 L随着北京时间变 化的数据，使用最小二乘拟合参数模型，求出拍摄该视频可能的地点。但是，上 述方法可能存在较大误差，究其原因，是因为并没有考虑到相机是否正对直杆和 影子所成平面，现构建误差修正模型来修正测量出的影子长度可能存在的误差， 从而修正拍摄地点。 （2）问题求解 12 把求出的影子长度导入 matlab（程序见附录）求解出可能的地点为： 纬度 经度 地点 北纬 39 度 31 分 东经 109 度 04 分 内蒙古巴彦淖尔市 5.4.2 误差修正模型的建立与求解： （1）模型建立： 假设由视频测量出的各时间段影子长度 C L与实际影子长度L相差一个修正 常数b，即有 C LbL 2 0 1( , , ) , ( , , ) LfT V fT V 得到 2 0 1( , , ) ( , , ) C LfT V Lb fT V ，将 相似原理模型中的 C L与北京时间 0 T导入 matlab 中， 使用最小二乘法拟合出未知 参数, , ,a b求得相对精确的若干可能地点。 若该视频没有日期，分析上述关系式 2 0 1( , , ) ( , , ) C LfT V Lb fT V ，未知参数 为, , ,a bV,在 matlab 中拟合出上述未知参数，求出拍摄地点与日期。 （2）模型求解： 在 matlab 中求解出参数, , ,a b计算得到若干修正后的地点如下表： 纬度 经度 地点 北纬 49 度 36 分 东经 104 度 22 分 蒙古色楞格 北纬 45 度 36 分 东经 110 度 蒙古赛音山达 北纬 48 度 21 分 东经 110 度 蒙古鄂嫩 5.4.3 结果分析 因为拍摄时相机不一定正对直杆与影子所成平面， 所以建立了上述修正模型 对地点进行修正，假设中由视频测量出的各时间段影子长度 C L与实际影子长度 L相差一个修正常数b，但实际上b是随着时间的变化而变化的，故上述误差修 正模型只能修正部分误差。 5.4.4 思考题 若拍摄日期未知，同样可使用蒙特卡洛拟合参数模型拟合出未知参数，但此 时拟合的未知参数有 4 个，会存在较大的误差，此时应该做大量的参数拟合，得 13 出拟合优度高且符合实际的参数，然后分类筛选出具有代表性的参数，确定若干 可能的地点。 六、模型优缺点六、模型优缺点 优点：优点：问题 1 建立了多元分析模型，画出影子长度与各个因素之间的图像，简洁 明了的刻画了太阳影子长度与各个因素之间的变化规律。问题 2、3 使用拟合优 度分析出拟合出的参数的合理性。 问题 4 建立了误差修正模型来修正测量出的太 阳影子长度存在的误差。 缺点：缺点：当拟合的参数较多时，最小二乘非线性拟合模型拟合出的参数会存在较大 的误差，问题 4 中误差修正模型拟合所得结果并不太好，究其原因可能是由于测 量出的太阳影子数据误差较大。 七、参考文献七、参考文献 1 陈晓勇 郑科科 对建筑日照计算中太阳赤纬角公式的探讨J 浙江建筑 2011 年 9 月 2 百度百科 赤纬角 -0w6U6M0w6U6M RXzzoc8wvFkGZDnSqltesqsKQdeqXYiV9GTVX0jxkRZ2vtneQKRXzzoc8wvFkGZDnSqltesqsKQdeqXYiV9GTVX0jxkRZ2vtneQK 2015/9/142015/9/14 3司守奎 孙玺菁 数学建模算法与应用 国防工业出版社 2011 年 93 P 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.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) （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*sin(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/18 0)*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)*pi)*t2).2)./(sin(12/18 15 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-(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:A22); 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.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);%把大量的拟合值导入 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*(xdata+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,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) . 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-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+canshu(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 h(p,j)=a(j); end xdata=xlsread(附件 4 数据.xlsx,Sheet1,A2:A12); y=2*sqrt(1-