1、 1 太阳影子定位 摘要 太阳影子定位技术就是通过分析物体的太阳影子长度变化，来确定物体所在的时间和地理位置。本文通过分析有关太阳影子各因素之间的关系，采用几何关系及 MATLAB软件编程等方法，对 所给 问题 分别给出了数学模型及处理方案。 对于问题一，根据题目 要求 ，首先确定影响影子长度的各个因素 （竿长，纬度，时间，日期） ，然后再根据几何知识确定它们之间的数学关系，建立相关的数学模型。再运用 MATLAB进行编程及绘出影长与各个变化因素的变化曲线图。 对于问题二，根据题目可知，在时间点 、 日期、 影子坐标已知的条件下，需要求 出所测点的地理位置，即经纬度。 我们根据问题一的相关结论
2、，做出 合理的假设。 根据附件 1中所给点求出影长， 用 MATLAB软件拟合出 所求点的影长与当地时间的关系曲线，确定各个影长所对应的当地时间 ，找到对应的北京时间。得到所求地与北京的时间差，即可用时 间差和经度的关系求得当地的经度。 在 此问题求解 中，我们 运用相关公式校准 坐标系，分析各个公式之间的相互转换，计算出题目所求地点的纬度。从而，确定当地的位置。 对于问题三，给定时间与影子的坐标，确定日期及地理位置。经度的确定与问题二中求得经度的方法一样。对于纬度的求解，则是运 用相关因素之间的公式，转换变化得出日期与纬度之间的关系。 再用 MATLAB软件 对变量（纬度） 进行穷举， 得到
3、最优解， 得出所求纬度， 确定 具体的地理位置 。 对于问题四，用 MATLAB软件分析视频，将视频处理成图片。同 样 ，用时间差来求出经度，并用公式算得纬度，以此来确定所测物的位置 及日期 。 最后，我们对于所建立的数学模型的优缺点做出了评价 关键词： matlab 影子变化 经纬度 三角变换 坐标矫正 2 一 问题重述 1 1 背景 确定视频的拍摄地点和拍摄 日期 是视频数据分析的重要方面，太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的 地点和日期 的一种方法。 1 2 需要解决的问题 问题一 ：建立影子长度变化 与各个参数的数学模型，并应用此模型解答给定条件（ 20
4、15年 10 月 22日北京时间 9:00-15:00之间天安门广场（北纬 39度 54分 26秒 ,东经 116 度 23分 29秒） 3米高的直杆）的影子变化曲线。 问题二： 根 据固定直杆在水平地面上的影子端点坐标，建立数学模型确定所处地点。依据 此数学模型确定所给 的影子顶点坐标数据的若干个可能地点。 问题三： 在 前一问基础上，若不给出日期，建立模型求解地点以 及日期。依据此模型确定所给条件的地点及日期。 问题四： 给出一根杆的影子视频，并且已知杆长 2m，建立模型确定拍摄地点。 拓 展： 建立模型确定拍摄地点及 日期。 二 问题分析 这属于 竿 影日照数学问题 。把 竿 顶影子端点
5、 坐标 移动的轨迹 ， 近似看成 太阳的 运行轨迹。运用 时角， 赤纬角，太阳高度角以及太阳方位角的角度关系进行推导。 2 1 问题一 影子的长度变化和时间 、 纬度、 日期的变化存在某种数学关系 ，因为有三个变量，在已知物体长度的情况下， 我们要建立影子和单个因素的关系，所以 必须先固定其中两个变量进而得出影子长度和另一变量的函 数关系 。在计算给定条件影子长度变化规律时，只需带入相应数学模型即可得到要求所需。 2 2 问题二 当给定影子端点坐标时，由于建立的坐标系 是随机的，所以必须先矫正成标准的坐标系， 由于 大致可认为影子的端点坐标连线是太阳的运动轨迹 反 方向，在此情况下就可以矫正坐
6、标系。要确定地点，一定要知道经纬度， 计算经度时，可3 以通过时差和经度差的关系进行求解 。计算纬度时，可以利用第一问的数学模型反过来求 解 。 2 3 问题三 本题在上一题的基础上， 要 求解日期 。同样，此题坐标是也不是标准的，也要进行矫正。 经度的算法可以参照问题二， 但是纬度的计算时没有日期，使得模型里有两个变量，不能直接求解。由于只有一个方程，所以 不 能直接求解两个未知量，但是两个未知量都是有 范围 的， 所以在这个范围内 对这两个变量同时进行穷举 ，进而解的最优解 。 2 4 问题四 本题给了一个视频，必须对该视频通过 matlab 进行处理，得到每一帧的图像，进而建立坐标系 得
7、到 影子的坐标。得到坐标后，可以利用问题二和问题三的模型进行求解。 三 模型假 设 1.假设题目中所给的数据全都真实可靠 。 2.假设影子定点坐标连线方向及为太阳运动方向的反方向。 3.假设大气折射对影子形成无影响。 4.假设太阳方位 角，太阳高度角，赤纬角是准确的。 5.假设北京正午 12： 00时影子最短。 6.假设地点都在北半球。 四 符号说明 表 1 符号说明 符号 说明 单位 时角 度 st 太阳时 小时 t 时差 小时 tt 北京时间 时 赤纬角 度 4 图 1 n 日期序号 太阳高度角 度 经度 度 当地纬度 度 A 太阳方位角 度 H 竿长 米 L 影长 米 五 模型的建立与解
8、决 5.1 问题一： 5 1 1 模型的准备 依据题目要求，要求建立影长和各个参数的变换规律，这个参数基本有纬度，时间，日期，竿长这几个变量，而影长和竿长的关系研究没什么意义，再次忽略，只研究影长和其他变量之间的关系。 5 1 2 模型的建立 依据题意，竿长已知，若已知太阳高度角，根据三角关系即可得到影子的长度。 依据竿影日照原理（如图一）可 以得到竿长，影长以及太阳高度角的三角关系 tanLH （ 1） 其中， H是竿的长度， L是影子的长度，是太阳高度角， 已知 太阳高度角的近似计算公式： coscoscossinsinsin （ 2） 其中， 是当地纬度， 是赤纬角， 是时角。 依据公式
9、 （ 2）可以 看出，要 求得 太阳高度角，必须 求得 赤纬角 和 时角 。 同时，我们知道太阳时角是以正午 12 点为零度开始计算，每小时 15度，上午为正，下午为负。因此可得时角公式： 5 ）12（15 st （ 3） 其中 ts为 太阳时。 赤纬角就是当时日期太阳直射的纬度。有计算公式： )3 65 ）n(2 8 42（s in2 3. 45 （ 4） 其中 n 为日期序号（如： 1 月 1 日 n=1）。 通过公式（ 2）（ 3）（ 4）的联立可 以 得到 太阳高度角的方程： c o s c o s c o s s i n s i n （a r c s i n ） （ 5） 进而我们可
10、得 tan的值，根据公式（ 1）可得影长模型： tanHL （ 6） 根据以上条件，可得数学模型 t a nHLc o s c o s c o s s i n s i n s i n)3 6 5）n2 ( 2 8 4（s i n2 3 . 4 5）12（15 st（ 7） 5.1.3 模型的求解 模型中及为影长 L 和时间 t，纬度 ，以及日期 n 的函数关系。 当其 中两个自变量确定后，就可建立影长 L 和另外一个自变量的模型。 （ 1） 影长 L 和时间 t 的模型 给定日期 n和纬度 ，模型就变成了影长 L和时间 t的一元函数， 应用 Matlab即可得到影长的变化曲线。 在此，我们验证
11、了赤道上 1月 1日的 影长 变化（如图 2） 由图 可以看出，当 1 月 1 日时， 9:00 到 15:00 的曲线为开口向上的抛物线，在早上 9:00 时，由于太阳 直射 南半球，所以影子长，到了当地正午 12:00 时影子最短，下午又开始增长，符合实际，模型基本成立。 6 图 3 赤道正午 12： 00 一年中影长的变化 （ 2）影长 L 和日期 n 的模型 当我们确定时间 t以及纬度 时 ， 模型就变成了影长 L和 日期 n的一元函数，应用 Matlab即可得到影长的变化曲线。 在此，我们验证了赤道正午 12： 00 一年中影长的变化（如图 3） 由图可以看出，一年内影子从第一天 开
12、始递减，当大致达到 3月份（北京六图 2 赤道 1 月 1 日的影长变化 0 50 100 150 200 250 300 350 4000 . 20 . 40 . 60 . 811 . 21 . 41 . 6日期 n影长L7 月份）时，影子最短，接近 0 m，而求一年内有两次接近于 0 m 的时候（即北京春分和秋分）。基本符合实际，模型基本成立。 （ 3） 影长 L 和纬度 的模型 当我们确定时间 t以及 日期 n时， 模 型就变成了影长 L和 纬度 的一元函数，应用 Matlab即可得到影长的变化曲线。 在此，我们验证了 1月 1日 正午 12： 00各个纬度 影长的变化（如图 4） 在图
13、 中验证北京 1 月 1 日正午 12： 00 的影长，基本符合实际，所以可认为建模基本成立。 （ 4） 本题要求 模型。 已知北京地理坐标是 北纬 39 度 54 分 26 秒 ,东经 116 度 23 分 29 秒 ，日期序号为第 295 天，竿长 3 m， 可以将这些数据直接带入 影长 L 和时间 ts 的模型 ，即可得以下曲线，如图 5。 图 4 1 月 1 日 正午 12： 00 各个纬度影长的变化 8 5.2 问题二： 5 2 1 模 型的准备 模型建立之前，我们分析数据得到所给影子顶点坐标并非 以 标准的 东西南北方向 坐标系下的坐标，所以我们必须进行矫正，把坐标 系修正成正南正
14、北的坐标系 。而后确定时差来确定经度，进而得到纬度。 5 2 2 模型的建立 （ 1）坐标系 矫正 由于 坐标 系建立 方式 不确定 ，所以必须进行 坐标系 矫正，得到 新 坐标系。 我们把给定的数据通过 matlab拟合成曲线，得到斜率 k，并通过三角关系式 ： k arctan （ 8） 得到偏转角度 （ 角如图 6） 。 进而通过坐标旋转公式 yxc o s s i n s i n c o s y x （ 9） 获得新的 坐标 （ x1,y1） 。 注 ： 矫正 坐标 系以东西方向为 x 轴，南北方向为 y轴。 图 5 北京 10 月 22 日 影长 变化 9 10 11 12 13 1
15、4 153 . 544 . 555 . 566 . 57时间 t影长L图 6 9 （ 2） 时差的求解 通过所给坐标在 matlab中进行拟合，得到一条 影长 L关于时间 ts的 抛物线方程： L=ats2-bts+c (10) 其中 L为影 长， t 为时间。 解出最低点坐标 t0，利用北京时间 12： 00时影子最短，利用比例关系 )12()( 0 ts tttt （ 11） 其 中 ts为 时间 ， t为 时差 ， t0为最低点时间， t 北 为对应的北京时间。 计算出时差 t 。 (3) 经度的求解 已 知两地经度相 差 1度，时间相差 4分钟， 所以可列出： 4北 st (12) 其
16、中 为当地经度， 北 为北京经度， ts为时差。 通过公式（ 12）解得经度 （ 4）纬度的求解 太阳方位角 就是太阳在方位上的角度，它和坐标有以下关系： tanAy 11 x （ 13） 其中（ x1， y1）为矫正还得坐标， A为太阳方位角。 得到太阳方位角： 11yarctanA x 又有太阳方 位角和 时角，赤纬角的关系： cos cossinsinA （ 14） 再结合公式（ 2），可得方程组： coscoscossinsinsincoscossinsinAtanAxy 11 （ 15） 10 即可解得纬度。 5 2 3 模型的求解 （ 1） 坐标系 矫正 通过 matlab 把 附
17、录 1数据进行拟合得到如图 7曲线，并获得斜率 k=0.147。 通过 公式（ 8）可得 =0.1459 通过 坐标旋 转 公式（ 9），把所有 坐标 矫正 ，得到以下 矫正 坐标（部分） 表 2 矫正 坐标（部分） x1 y1 x1 y1 x1 y1 1.097820723 0.341217234 1.272291262 0.344065235 1.461264662 0.345184945 1.131680076 0.341898561 1.308821889 0.344353761 1.501074729 0.345092317 1.16603411 0.342507147 1.3459
18、60682 0.344653933 1.541705381 0.345081175 1.200996308 0.34312738 1.383806575 0.344951205 1.461264662 0.345184945 1.236339702 0.343590484 1.4221326 0.3450768 1.583014038 0.34486926 （ 2） 时差 及经度 的求解 将所给坐标在 matlab 中进行拟合，得到 影长 L关于时间 ts的 抛物线方程 L=0.1483ts2-3.4541ts+20.5245 （ 9）， 拟合曲线如图 6， 解得最低点时间坐标为 ts0=11
19、.5984，同时所给数据第一个坐标对应的时间坐标为 ts=13.8760。进而解得 到时差 t=0.4224 。因为已知北京经度大约为 116 ，通过方程（ 11）可求出该点经度 =110.04 。 图 7 11 （ 3）纬度的求解 由于 模型（ 15）求解 比较复杂，程序不能直接解 出 纬度角 ，于是我们在 纬度区间 0 /2之间间隔 0.001取值依次带入方程，并求解出最优解。解得结果如表 3： 表 3 纬度 纬度 24.579 24.6944 24.8090 1.2032 1.3178 1.43239 1.60428 25.3820 25.49 2.0626 2.23453 2.4064
20、2 26.069 26.241 2.97938 26.5279 26.6998 3.609 3.8388 27.2727 27.4446 从表中可以看出，纬度有两个大致范围 1.2 3.8 ， 24.5 27.4 。差距不是很大，取平均值，纬度为 2.5 和 25.9 。 （ 4）检验 把求得的纬度带入问题一的影长和时间模型中，检验所给数据中第一个时间 14 27时（ 14： 16）的影长数据 与图像中的值对比。如图 9，图 10 图 8 影长 拟合 12 结果在误差范围 内 ，可认为此点符合要求 （ 经度 110.04E ，纬度 25.9 N 地点： 广西壮族自治区桂林市龙胜各族自治县 ）
21、图 9 检验 纬度 25.9 图 10 检验 纬度 2.5 13 结果在误差范围内 ，可认为此点符合要求（经度 110.04E ，纬度 2.5 N 地点：马来西亚北海岸）。 5.3 问题三： 5 3 1 模型的准备 和问题二一样，建立模型时， 分析数据得到所给影子顶点坐标并非标准 的正南正北方向建立的坐标 系 ，所以我们必须进行矫正，把坐标 矫正成南北向的坐标系 。而后确定时差来确定经度，进而得到纬度。 5 3 2 模型的建立 （ 1） 坐标 矫正 同问题二一样， 坐标系 的 建立不确定，所以必须 进行标准矫正，得到 矫正 坐标系。通过附件 2里的坐标，我们可以通过 matlab 程序 获得斜
22、率 k结合公式 获得偏转角 。并 通过 公式 （ 9）进行坐标矫正 获得新的 矫正 坐标 （ x2， y2）。 （ 2） 时差的求解 通过所给坐标在 matlab中进行拟合，得到一条 影长 L关于时间 ts的 抛物线方程 (10),解出最低点坐标 ts0，利用北京时间 12： 00时影子最短，利用比例关系 式 （ 11） 计算出时差 t。 (3) 经度的求解 已知经度相差 1度，时间相差 4分钟，所以可列出： 经度求解公式 ( ) 并 通过 该公式 解得 地点 经度 （ 4）纬度的求解 本题在求解纬度的过程和问题二大致相同，但是在未知数上多了一个日期的自变量，大大加大了解题难度。但是仍可用第二
23、问的模型进行求解，模型为： coscoscossinsinsincoscossinsinAtanAxy 111x(16) 5 3 3 模型的求解 1.附录 2 的其求解 （ 1） 坐标系矫正 14 图 11 通过 matlab把附录 2里的 数据进行拟合得到如图 11曲线，并获得斜率 k=0.5372。 通过 公式（ 8）可得 =0.4930 通过 坐 标旋转公式（ 9），把所有 坐标 矫正，得到以下矫正坐标（部分） 表 4 矫正坐标（部分） x1 y1 x1 y1 x1 y1 -1.00624245 0.736969584 -0.85204970 0.742108341 -0.7043148
24、5 0.743678412 -0.97479728 0.738239028 -0.8219652 0.742746611 -0.67550308 0.743638181 -0.94371104 0.739474265 -0.79219901 0.743215255 -0.64687405 0.743469091 -0.91285487 0.740492552 -0.76266282 0.743466949 -0.61847509 0.743083051 -0.88226954 0.741429305 -0.73339747 0.74363711 -0.59025887 0.742568152
25、（ 2） 时差及经度的求解 将所给坐标在 matlab 中进行拟合，得到 影长 L关于时间 ts的 抛物线方程 L=0.0902ts2-3.1845ts+26.4555 拟合曲线如图 12， 解得最低点时间坐标为 ts0=16.2242，同时所给数据第一个坐标对应的时间坐标为 ts=14.7429。进而解得到时差 t=2.1812 。因为已知北京经度大约为 116 ，通过方程（ 11）可求出该点经度 = 83.6616 。 15 （ 3）纬度的求解 由于模型（ 15） 有两个变量， 求解比较复杂，程序不能直接解出纬度角 和具体日期， 于是我们在纬度区间 0 /2之间间隔 0.001取值 ，日期
26、区间 1365之间间隔 1取值 依次带入方程，并求解出最优解。解得结果如表 3： 表 5 纬度 纬度 1.7307 5.9340 6.6758 8.6538 9.1483 9.3956 10.631 12.115 18.049 19.532 20.769 25.219 26.208 28.681 34.862 36.593 36.840 37.087 39.065 39.807 40.549 表 6 日期序号 日期序号 294 291 202 143 64 167 178 73 262 256 121 240 110 109 133 140 172 209 153 164 181 由于 数据庞
27、大，要一一验证，应用上题验证方法进行验证，可得纬度大致为 36.8N，经度为 83.66 E，日期序号在 160天左右。 位置大致为 新疆维吾尔自治区和图 12 16 田地区民丰县 。 2.附录 3 的其求解 和附录 2 的解法完全相同，得到的相应数据如下： 斜率 k=-0.0435,旋转角 =-0.0435,进而到矫正坐标如表 7 表 7 矫正坐标（部分） x1 y1 x1 y1 x1 y1 1.112001542 3.383443038 1.390551382 3.369242962 1.680759226 3.364065307 1.16694705 3.37984892 1.44750
28、3596 3.36753845 1.740578163 3.364189627 1.222274789 3.376671815 1.505029155 3.366059837 1.801066001 3.364744006 1.2778892 3.373807471 1.563032501 3.365102491 1.862222741 3.365628537 1.333981397 3.371264583 1.621609191 3.364371044 1.924143941 3.366847568 拟合得到的影长 曲线如图 13， 方程为： L=0.2964ts2-7.8768ts+5.8
29、064 解得最低点时间坐标为 ts0=13.2854，同时所给数据第一个坐标对应的时间坐标为 ts=13.7297。进而解得到时差 t=0.7056。因为已知北京经度大约为 116 ，通过方程（ 11）可求出该点经度 =105.8 。 图 13 17 同样一个用附录 2的解题步骤解得附录 3中的纬度及日期如下表 8 表 8纬度及日期 纬度 4.203 28.18 45.741 30.412 48.4615 19.2857 50.4395 5.934 16.5659 52.41 46.73 52.664 21.510 30.1648 46.2362 11.1263 7.664 7.664835
30、23.98 12.11 27.692 日期 313 254 193 254 184 59 174 330 47 176 221 191 54 73 114 22 341 343 290 22 283 由于 数据庞大，要一一验证，应用上题验证方法进行验证，可得纬度大致为 30.2N，经度为 105.8 E，日期序号在 180天左右。位置大致为 重庆市潼南县 。 5.4 问题四 ： 5 4 1 模型的准备 此题给了一个 40多 分钟的视频，可以将视频应用 matlab将视频的帧按图片截取出来， 并应用 matlab对图片进行批量处理得到竿顶坐标以及影子的像素长度 。最后应用第一问，以及第三问的模型
31、进行求解 5 4 2 模型的建立 （ 1） 根据公式（ 1），可得太阳高度角 LHtanarc 又 根据公式（ 2）可得 c o s c o s c o s s i n s i n ）t a na r c（s i n LH （ 17） 进而建立了影 长 L和 纬度 的关系，只要 影长固定，即可求的对应的纬度。 在计算经度时，模型和上面模型建立过程相同。先拟合出影长和时间的图像，得到 影长 L 关于时间 ts的 抛物线方程 （ 10）。然后求解时间差 t，然后根据经度与时间的关系（ 12）得到该地经度。 模型为： 4)12()(北0sstttttt(18) （ 2） 经度的求解和上面的建立过 程
32、一致，先求时差，再求经度。 18 求解纬度的时，仍然可以利用模型（ 17）进行求解，但是要注意有两个自变量，求解 时要利用 matlab 进行求解。 5 4 3 模型的求解 （ 1）依据题意可知竿长 2 m， matlab进行数据处理得到一系列影长如表 9 表 9 影长（部分数据） 影长 2.3811 2.3782 2.3752 2.3574 2.3516 2.3486 2.3456 2.3280 2.3220 2.3043 2.3043 2.2983 2.2836 2.2777 2.2571 2.2512 2.2334 2.2363 2.2274 2.2157 2.2068 2.1979 2
33、.1920 2.1743 2.1625 2.1565 2.1535 通过 matlab可以拟合出影长和时间的关系 L=0.0815ts2-2.3284ts+16.6742 (19) 解得影长最短是 14.292。进而求的时差为 t=2.2726，经度为 82.29 。 求解纬度时，根据方程（ 17）可以拟合出影长 L和纬度的图像，如图 14。 由图可知同一个影长对应两个纬度，在图上取点验证可得纬度为 25.21和 64.17。 地点为 （ 25.21 N， 82.29E 印度）及（ 64.17 N， 82.29E 俄罗斯） 图 14 影长 和 纬度 19 （ 3） 因为时间和影长不变，所以经度
34、和（ 1）中相等 ，为 82.29 。 在模型（ 17）中， )365 ）n2 ( 2 8 4（sin2 3 . 45 ， 有一个变量日期 n，同时，方程中纬度 为自变量，此模型在该问中有两个未知量，所以不能解出具体数值。 欲求的日期和纬度，我们采用第三问的迭代法，在纬度范围 0 /2 之间以间隔 0.001 取值，在日期 1365 之间间隔为 一取值，利用 matlab 循环求解，得到最优解。得到的纬度和日期如表 10 表 10 纬度和日期（部分数据） 纬度 69.230 68.241 67.252 66.016 65.027 64.285 63.296 62.307 61.565 68.9
35、83 67.994 67.005 65.769 64.780 64.038 63.049 62.060 61.318 68.736 67.747 66.758 65.521 64.537 63.543 62.802 61.813 61.071 日期 172 172 172 172 172 172 172 172 172 172 172 172 172 172 172 172 172 172 172 172 172 172 172 172 172 172 172 从表中可看出纬度为 68，日期序号为 172。 地点（ 68N ， 82.29E 俄罗斯）。 六 模型的评价 6.1 模型的 优点 （
36、 1）本模型结构调理清晰，运算简单，可以正确的反映影子和各个变量之间的关系 。 （ 2） 本模型考虑 的因素较多，考虑比较详细，只要对相关数据进行修改，就可求出对应的参数指标 。 （ 3）本模型可以很好的利用影子长度变化确定物体的所在地的位置以及日期，有很好的利用价值。 6.2 模型的 缺点 本模型没有很好的考虑大气折射问题，而且利用了一些理论值，存在些许误差。 20 七 参考文献 【 1】寿纪麟，数学建模 -方法与范例方法与范例，西安交通大学出版 (1993) 【 2】郑鹏飞，林大钧，刘小羊，吴志庭，基于影子轨迹线反求采光效果的技术研究，华东理工大学学报 (自然科学版 )（ 2010） 【
37、3】薛定宇，陈阳泉，高等数学问题的 MATLAB 求解，北京清华大学出版社（ 2008） 【 4】尚月强，杨一都， Matlab 及其在数学建模中的应用，贵州师范大学学报 (自然科学版 )（ 2005） 【 5】郑鹏飞，基于影子轨迹线反求采光效果的技术研究 21 附录 本文 通过 matlab 建立数学模型 ，以下为 代码； 附录一： 求解影长和各个因素的关系 (1)影长和时间关系 n=input(输入日期序号 n=); H=input(输入杆长 H=); w=input(输入纬度 w=); t=9:0.01:15; %时间区间 q=24.45*sin(2*pi*(284+n)/365); %
38、q=赤纬 e=15*(t-12); %e=时角 h=sin(w)*sin(q/180*pi)+cos(w)*cos(q/180*pi)*cos(e/180*pi);%太阳高度角 root=sqrt(1-power(h,2); u=h./root; l=H./u; %影长 plot(t,l,g); (2)影长和日期的关系 t=input(输入时间 t=); H=input(输入杆长 H=); w=input(输入纬度 w=); n=1:365 ; %日期区间 q=24.45*sin(2*pi*(284+n)/365); %q=赤纬 e=15*(t-12); %e=时角 h=sin(w)*sin(
39、q/180*pi)+cos(w)*cos(q/180*pi)*cos(e/180*pi);%太阳高度角 root=sqrt(1-power(h,2); u=h./root; l=H./u; %影长 plot(n,l,g); （ 3） 影长和纬度的关系 n=input(输入日期序号 n=); H=input(输入杆长 H=); t=input(输入时间 t=); w=0:0.01:(pi/2) ; %纬度区间 q=24.45*sin(2*pi*(284+n)/365); %q=赤纬 e=15*(t-12); %e=时角 h=sin(w)*sin(q/180*pi)+cos(w)*cos(q/18
40、0*pi)*cos(e/180*pi);%太阳高度角 root=sqrt(1-power(h,2); u=h./root; l=H./u; %影长 z=w/pi*180; plot(z,l,g); 附录二： 求附件 1中数据的地理位置 p = polyfit(x,y,1);%求原始数据图线斜率 22 th = atan(p(1,1); %th旋转角 xy = cos(th),sin(th);-sin(th),cos(th);*x;y;%矫正坐标 x1 = xy(1,:);%校正后横坐标 y1 = xy(2,:);%校正后纵坐标 l = sqrt(x1.*x1+y1.*y1);%影长 %拟合影长
41、与时间关系方程 xt=13.7:0.05:14.7; yt=l; n=2; p1 = polyfit(xt,yt,n);%二次方程的系数 xl=7:0.1:16;%取时间段 xl funl = p1(1,1)*power(xl,2)+p1(1,2)*xl+p1(1,3);%当地影长与时间关系拟合方程 p2 = polyfit(xl,funl,2); xx = min(xl):0.1:max(xl); yy = polyval(p2,xl); Y = vpa(poly2sym(p2,x),6);%转化为多项式格式 Y1 = diff(Y);%求一阶导数 Xmin = eval(solve(Y1)
42、;%一阶导数等于零，极值点 Ymin = polyval(p2,Xmin); plot(xl,funl,g)%做出拟合方程的图像并标出极值点 time1 = solve(0.1483*time12-3.4541*time1+20.5245=1.149625826,time1); jingdu = 116.38-(14.7-(12+(time1(2)-Xmin)*15%经度 shicha = 14.7-(12+(time1(2)-Xmin) for i = 1:1:21; ss(i) = funwei(x(i),y(i),t(i); hh(i) = funH(x(i),y(i),t(i); en
43、d sum1 = ss(find(ss20); sum2 = ss(find(ss20); sum3 = hh(:); weidu1 = mean(sum1(:) weidu2 = mean(sum2(:) H = mean(sum3(:) 附录三： (1)求解附件 2， 3中数据的地理位置 p = polyfit(x2,y2,1);%求原始数据图线斜率 th = atan(p(1,1); %th旋转角 xy = cos(th),sin(th);-sin(th),cos(th);*x2;y2;%矫正坐标 x1 = xy(1,:);%校正后横坐标 y1 = xy(2,:);%校正后纵坐标 l =
44、 sqrt(x1.*x1+y1.*y1);%影长 %拟合影长与时间关系方程 23 xt = 13.7:0.05:14.7; yt = l; n = 2; p1 = polyfit(xt,yt,n);%二次方程的系数 xl = 8:0.1:24;%取时间段 xl funl = p1(1,1)*power(xl,2)+p1(1,2)*xl+p1(1,3);%当地影长与时间关系拟合方程 p2 = polyfit(xl,funl,2); xx = min(xl):0.1:max(xl); yy = polyval(p2,xl); Y = vpa(poly2sym(p2,x),6);%转化为多项式格式
45、Y1 = diff(Y);%求一阶导 Xmin = eval(solve(Y1);%一阶导数等于零，极值点 Ymin = polyval(p2,Xmin); plot(xl,funl,g);%做出拟合方程的图像并标出极值点 time1(:) = solve(0.09902*time12-3.17*time1+26.46=1.2473,time1); jingdu = 116.38-(12.7-(12-(Xmin-time1(1)*15;%经度 shicha = 12.7-(12-(Xmin-time1(1);%时差 for i = 1:1:21; NN(i) = funN(x2(i),y2(i
46、),t2(i); ww(i) = funwei2(x2(i),y2(i),t2(i); end (i),t2(i); end 附录四 ： 求解视频中的地理位置 p1 = polyfit(t4,l4,2);%二次方程的系数 xl=0:0.1:24;%取时间段 xl funl = p1(1,1)*power(xl,2)+p1(1,2)*xl+p1(1,3);%当地影长与时间关系拟合方程 p2 = polyfit(xl,funl,2); xx = min(xl):0.1:max(xl); yy = polyval(p2,xl); Y = vpa(poly2sym(p2,x),6);%转化为多项式格式
47、 Y1 = diff(Y);%求一阶导数 Xmin = eval(solve(Y1);%一阶导数等于零，极值点 Ymin = polyval(p2,Xmin); plot(xl,funl,o,xx,yy,Xmin,Ymin,*)%做出拟合方程的图像并标出极值点 time1 = solve(0.0815*time12-2.3284*time1+16.6742=2.381190669,time1); jingdu = 116.38-(t4(1)-(12-(Xmin-time1(1)*15;%经度 shicha = t4(1)-(12-(Xmin-time1(1); for i=1:1:61; NNN(i)=funN3(l4(i),t4(i); www(i)=funwei3(l4(i),t4(i); end 24 sum = www(:); weidu = mean(sum(:); 附录五 ： 求解纬度 function I = funwei( x,y,t ) %输入影长坐标，求出对应纬度 % 此处显示详细说明 th = at
