2015年全国大学生数学建模竞赛A题国赛一等奖 23页

太阳影子定位 摘要 视频的拍摄地点和拍摄日期的确定是视频数据分析的重要方面，通过太阳影子定位 技术分析物体影子的变化，是对视频拍摄地点和时间确定的重要方法。本文应用地理中 太阳高度角、方位角、赤纬角等知识建立影子长度变化模型，通过最小二乘法拟合建立 地点和日期的优化模型。 针对问题一，利用太阳高度角，得到影长与杆长、太阳高度角的关系 )tan(hs Lk  ， 所涉及参数主要有：地理纬度、经度、日期、时间、杆高。根据地理领域相关知识，计 算太阳高度角，建立影子长度变化模型。将题目所给数据代入模型，应用 matlab 计算 得到影子长度变化曲线（图 3）。通过控制变量法分析影长与各参数的变化规律。 针对问题二，基于问题一的影长模型，分析其他条件一定而杆高不确定时的影长与 经纬度的关系。建立优化模型，通过最小二乘法拟合影长和影子端点坐标，求解参数。 得到结果：纬度 ''96.6'1619 ，经度 ''00.60'35108 ，海南省直辖县级行政单位东方市。 针对问题三，与问题二相比缺少日期参数，建立优化模型，用拟合方法确定未知地 点与日期。对于附件 2 得到 3 个可能地点，分别为纬度 ''07.4'3237 ，经度 ’‘00.60'1778 ， 新疆维吾尔自治区和田地区皮山县，2015 年 7 月 1 号；纬度 '11.12'3739 ‘ ，经度 ''00.36'1279 ，新疆维吾尔自治区喀什地区巴楚县，2015 年 7 月 14 号；纬度 '04.56'2635 ‘ ， 经度 '00.36'2477 ‘ ，喀什米尔，2015 年 6 月 23 号。对于附件 3 得到 3 个可能地点，纬 度 ''06.14'0847 ，经度 ''00.60'05112 ，蒙古，2015 年 2 月 6 号；纬度 ''60.44'2827 ，经 度 ''00.60'47110 ，湖南省邵阳市隆回县，2015 年 1 月 12 号；纬度 ''76.24'4866 ，经度 ''0'24113 ，俄罗斯，2015 年 3 月 3 日。 针对问题四，要求确定视频拍摄地点和日期，建立优化模型。利用图像处理技术， 对视频图像分帧，并转化为灰度矩阵，二值处理后得到二值图像，界定阀值提取影长。 用 matlab 软件计算影子端点坐标。将影子端点坐标值进行斜坐标系与直角坐标系的变 换，和时间值一起代入模型求解。当日期不确定时，用问题三模型求解。日期确定时用 问题二模型求解。参数结果：纬度 ''52.4'5442 ，经度 ''62.48'17112 ，内蒙古自治区锡 林郭勒盟苏尼特右旗； 。 关键字：最小二乘法拟合 太阳高度角 控制变量法 图像处理 灰度矩阵 一、问题重述 随着视频监控技术的应用与发展，视频数据的分析越来越重要，如何确定视频的拍 摄地点和拍摄日期是视频数据分析的重要方面。太阳影子定位技术就是通过分析视频中 物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。 基于太阳影子定位技术，解决以下问题： 1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用 模型画出 2015 年 10 月 22 日北京时间 9:00-15:00 之间天安门广场（北纬 39 度 54 分 26 秒,东经 116 度 23 分 29 秒）3 米高的直杆的太阳影子长度的变化曲线。 2.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆 所处的地点。由已知的影子顶点坐标数据，根据模型给出若干个可能的地点。 3.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆 所处的地点和日期。将模型应用于已知的的影子顶点坐标数据，给出若干个可能的地点 与日期。 4．根据一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的 高度为 2 米。请建立确定视频拍摄地点的数学模型，并应用模型给出若干个可能的拍摄 地点。若日期未知，能否根据视频确定出拍摄地点与日期。 二、问题分析 2.1 概论 根据太阳影子定位技术确定视频中物体位置和日期，这是一个涉及图像处理的最优 化问题。问题涉及地理、物理、数学等多学科知识，问题难点在于太阳影子定位过程涉 及多个参数的彼此关联，通过拟合确定变量。问题特点在于对涉及视频图像的处理，提 取关键性信息。 2.2 问题一 问题一要求建立影子长度变化的数学模型，通过查阅太阳影子定位技术的相关资 料，我们了解到影子的长度变化与物体高度、经度、纬度、时间、日期等参数密切相关。 因此解决问题一的关键在于确定与影子长度变化有关的参数以及影长和参数的关系，建 立影子长度关于各参数的关系模型。以此分析影子长度关于各参数的变化规律，并代入 实例计算。 2.3 问题二 问题二给出杆的影子顶点坐标，要求确定杆的地理位置。由于杆的高度未知，单独 应用问题一的模型难以解决问题二。因此，我们考虑在杆高度未知情况下，可通过太阳 方位角相关知识，以影长与太阳方位角和影子顶点坐标的关系，结合所给实际数据以杆 的经纬度为所求量，进行拟合，分析杆的地理位置。 2.4 问题三 问题三要求根据直杆在太阳下的影子长度端点坐标数据，确定直杆所处的地点和日 期。地点、日期、杆高均未知，与问题二类似，所以考虑应用问题二的拟合算法，建立 优化模型，确定未知的各个参数值。 2.5 问题四 问题四要求确定视频的拍摄日期与地点，日期未知时可以根据视频中直杆的影子变 化用问题三的模型求解，日期已知时可用问题二的模型求解。但由于所需影子信息都存 储在视频图像中，首先需要对视频图像进行处理，提取每帧画面中影子的长度、时间等 信息。之后用优化模型进行求解。 三、模型假设 1.假设题目所给数据坐标系以 x 轴正方向为正东，y 轴正方向为正北； 2.假设题目所给数据合理可靠； 3.假设物体影子长度不受温度、风等因素影响； 4.假设问题四视频中时间为当地时间； 5.假设光的传播不受空气折射率的影响。 四、符号说明与名词解释 K：物体影子长度； L：被照射物体的高度； hs ：太阳高度角，太阳高度角是指某地太阳光线与该地作垂直于地心的地表切线的夹 角；  ：太阳赤纬，太阳赤纬又称赤纬角，是地球赤道平面与太阳和地球中心的连线之间 的 夹角；  ：时角，时角是指地平线在地球赤道平面上的投影与当地时间 12 点时、地中心连 线在赤道平面上的投影之间的夹角； As :太阳方位角，即太阳所在的方位，指太阳光线在地平面上的投影与当地经线的夹 角，可近似地看作是竖立在地面上的直线在阳光下的阴影与正南方的夹角； x ：影长在正东方向 x 轴的分解坐标； y ：影长在正北方向 y 轴的分解坐标; N ：积日，所谓积日，就是日期在年内的顺序号；  ：地理纬度；  ：日角； 五、模型的建立与求解 5.1 问题一 影子长度变化模型 由于太阳光对地球表面的照射，地球表面物体会形成与照射太阳光同方向的影子， 随着日期、时间的变化，影子长度会发生相应的改变。由于物体地理位置的不同即经纬 度的不同，相同形状、大小的物体影子长度也会有所不同。根据太阳光照射地球表面物 体形成影子的原理，我们建立影子长度变化模型，分析影子长度关于各个参数的变化规 律。 5.1.1 影子长度变化模型的建立 我们以一根直杆形成的影子为例。太阳光照射直杆形成影子的情形可抽象成一个直 角三角形，如下图所示。 图 1 直杆在阳光下的投影 根据上图可知，杆的影长与杆高度的关系满足， )tan(hs Lk  。根据太阳高度角的定 义，角hs 即为该时刻的太阳高度角（本问计算过程所有的角度均以弧度值计算）。 图 2 太阳高度角、赤纬角、地理纬度的示意图 1.太阳高度角的计算 根据太阳高度角的定义，角hs 即为该时刻的太阳高度角。太阳高度角随着地方时和 赤纬的变化而变化。由参考文献[1]知太阳高度角的计算公式如下：  coscoscossinsinsin hs （1） 其中 为地理纬度， 为太阳赤纬， 为时角。 （1）由参考文献[1]知太阳赤纬的计算   0.0201cos30.3656cos20.758cos- 0.1712sin3-0.1149sin2sin2567.233723.0   （2） 式中 为日角， 02422.365/2 NNtt  ， （3） N 为积日，平年 12 月 31 日的积日为 365，闰年则为 366。 ]4/)1985[()1985(2422.06764.790  年份年份 INTN （4） （INT 表示向下取整数部分） （2）时角的计算 北京时间是中国采用北京时区作为标所在的东八准时间，是东经 120°地方的地方 时间。由于地球上每相差一个经度就会有 4 分钟时间差，当某地当地时间为 h 时，北京 时间为 H，设当地经度为 E(以东经度为正值，西经度为负值）。 由参考文献[1]知： )120(4 EHh  （5） 时角计算公式： )12()24/360(  h （6） 2.影子长度的计算 由图 1 可知,影子长度 K 的计算公式为： )]coscoscossin(sintan[)tan(  arc L hs Lk （7） L 为直杆的高度。 5.1.2 模型的求解 应用影子长度的变化模型，我们对 2015 年 10 月 22 日北京时间 9:00-15:00 之间天 安门广场（北纬 39 度 54 分 26 秒,东经 116 度 23 分 29 秒）3 米高的直杆的太阳影子长 度的变化情况进行分析。 1.太阳赤纬的计算 由题目所给数据知，日期为 2015 年 10 月 22 日。根据闰年判断方法可判断 2015 年 为平年，即 2015 年共 365 天，其二月份共 28 天，则 10 月 22 日的积日 N=295。 ]4/)19852015[()19852015(2422.06764.790  INTN =79.9424， 故 t=215.0576，日角  3.6977°=3.6977 /180=0.0645rad， 太阳赤纬 rad0258.047989.1 0.0201cos30.3656cos20.758cos- 0.1712sin3-0.1149sin2sin2567.233723.0      2.时角的计算 该地地理经度值 E= '''2923116 ，以度为单 E= 3600/296023116 ）（  =116.3914°。 同理，纬度 3600/26605439 ）（  =39.9072°，转化为弧度表示  =39.9072 /180=10.6962rad 根据地球经度差与时间差的关系，当北京时间 H 为 9:00,当地时间(以时表示) h=9-4/60 (120-116.3914)=8 .7594 换算成具体时间为 8:45:33。 同理，北京时间 H 为 15:00 时，当地时间(以时表示)h= 14.7594，换算成具体时间为 14:45:33。综上可知时角 )12()24/360(  h ，h 为 8:45:33-14:45:33。 3.影长的计算 将相关参数代入公式（7），得到 )])]12(15cos[cos0.02586962.10cos0258.0sin6962.10(sintan[ 3 )tan(  harchs Lk (h 为 8:45:33-14:45:33) 通过 matlab 求解可得到如下图所示的影子长度变化曲线。 图 3 影子长度变化曲线 主要时间点的影子长度如下表所示： 表 1 2015 年 10 月 22 日部分时间点的影子长度 时间 9:00 9:30 10:00 10:30 11:00 11:30 12:00 影子长度 7.4172 6.0471 5.1386 4.5175 4.0983 3.8372 3.7123 时间 12:30 13:00 13:30 14:00 14:30 15:00 影子长度 3.7147 3.8444 4.1110 4.5369 5.1670 6.0889 基于影子长度变化模型，对影子长度与各参数变化规律的分析： 1.影长与时间 以当地时间标定当地正午，以当地正午时间为轴，轴两侧时间点的影长变化曲线基 本成左右对称结构，且从太阳升起时刻到正午时刻影长递减，正午时刻至日落时刻影长 递增。如图 3 所示。 2.影长与杆长 根据公式（7）可知，在其他变量一定的情况下，影长 k 与杆高 L 正相关，即在同 一天的同一时刻同一地点，物体的高度值越大，影子长度越大。 3.影长与纬度 影长与纬度的变化关系如下图所示： 图 4 影长随纬度变化关系 在经度、日期、时间一定的情况下，影长随纬度变化规律为：影长的变化趋势关于赤道 对称，纬度逐渐接近赤道时，影长先增大后减小，在北纬 57°和南纬 37°左右取得极 大值。 5.2 问题二 物体位置确定模型 问题二给出了影子的顶点坐标，杆的高度未知，无法完全应用第一问的模型求解。 因此，我们在问题一的太阳高度角分析的基础上，提出与之有关的太阳方位角的分析方 法，分析影长与太阳方位角和影子顶点坐标的关系。以物体所在经纬度为未知量，利用 最小二乘法拟合实际数据坐标和影子长度变化函数中的坐标值，进而确定物体经纬度参 数。 5.2.1 模型的建立 本文选定以 x 轴正方向为正东，y 轴正方向为正北，建立二维直角坐标系。 图 5 太阳方位角与太阳高度角 根据太阳方位角的计算公式：   coscos sinsinsincos hs hsAs  （8） 将杆的影子长度在 x,y 轴分解，得到如下图： 图 6 影长的正向分解 由图 6 可得， Asky Askx cos sin   （9） 影长是与参数有关的变量，结合已知的同一地点同一时刻的相同杆的影长坐标值（ ii yx , ） 运用最小二乘法拟合影长和坐标值建立问题二的优化模型，不仅考虑影长也考虑了角 度。 问题二的优化模型为：    n i iii LhNEyyLhNExxLhNEkk 1 222 ])),,,,(()),,,,(()),,,,([(min  （10） 由上式可得到确定的经度值 E、纬度值 ，即可确定所求的地点。 5.2.2 模型的求解 问题二给出了若干组 2015 年 4 月 18 日 14:42-15:42 时间段的影子顶点坐标值 （ ii yx , ），将之与问题一模型中的坐标（ yx, ）拟合，通过 matlab 软件求参数，具体 过程如下所示。 图 7 matlab 中 cftool 计算参数 图 8 参数输入 求得经纬度的参数值:纬度 ''96.6'1619 ，经度 ''00.60'35108 ，地点：海南省直辖县 级行政单位东方市。 图 9 海南省经纬度定位 5.3 问题三 日期和地点确定模型 问题三在问题二的基础上加大难度，增加未知量日期，要求确定地点和日期两个未 知参数，应用问题二的拟合方法就可以确定未知地点和日期。 问题三优化模型如下：           Asky Askx LhNEyyLhNExxLhNEkk n i iii cos sin ])),,,,(()),,,,(()),,,,([(min 1 222  （11） 将所给数据，某地某日北京时间 13:09-14:09 时间段以及某地某日北京时间 12:41-13:41 时间段的坐标值（ ii yx , ），代入上述模型用 matlab 软件求解： （1）附件 2 数据求得 3 个可能地点： 1 纬度 ''07.4'3237 ，经度 ’‘00.60'1778 ，新疆维吾尔自治区和田地区皮山县，2015 年 7 月 1 号； 2 纬度 '11.12'3739 ‘ ，经度 ''00.36'1279 ，新疆维吾尔自治区喀什地区巴楚县，2015 年 7 月 14 号； 3 纬度 '04.56'2635 ‘ ，经度 '00.36'2477 ‘ ，喀什米尔，2015 年 6 月 23 号。 （2）附件 3 数据求得 3 个可能地点： 1 纬度 ''53.7'1436 ，经度 ''32.54'11111 ，山西省临汾市蒲县，2015 年 1 月 12 号 2 纬度 ''23.28'5846 ，经度 ''27.43'06111 ，蒙古，2015 年 2 月 9 号 3 纬度 ''73.25'5340 ，经度 ''74.40'12111 ，呼和浩特，2015 年 3 月 9 号 图 10 呼和浩特经纬度定位 5.4 问题四 视频拍摄地点和日期确定模型 问题四的数据为视频图像，我们需要从视频中提取关键信息，即影子长度值和对应 的时间值，将影子长度转化为坐标（x,y）并建立优化模型。若日期未知，则应用问题 三的日期和地点确定模型；若日期已知，应用问题二的优化模型求解。首先将视频分帧， 转化为灰度矩阵，二值处理得到二值图像。通过界定阀值提取影长，用 matlab 软件选 择每帧画面的影子端点坐标值和杆底端坐标并记录时间值。 5.4.1 模型的建立 视频是由帧画面快速变化形成的，一段时长的视频图像可分为若干帧画面，每一帧 画面看作一张图片，对应一个时刻。我们将视频分帧计算，每帧画面中的信息为固定值， 可提取得到对应某一时刻的关键信息。 1.帧画面的选择 由于人眼的视觉暂留效应，为使视频画面平滑连续，连续图像的变化达每秒 24 帧 以上。我们提取数据进行计算求解时，不必提取视频中全部帧画面，合理提取使得帧画 面间时间间隔满足计算精确度即可。如本题视频时长大概 40 分钟，我们每隔两分钟提 取一帧画面的数据，共提取 21 帧画面。 2.画面像素点与实际长度单位米的比例转换 画面中的图像以像素点表示，实际场景中影长以米为单位表示，对画面中的杆长进 行像素点计算得到像素点值 P,则有画面中每一个像素点与长度单位米的比例关系： q=L/p,表示每一像素点对应实际长度 q 米。 3.影子端点坐标的坐标系转换 提取的帧画面中，可以得到影子端点的坐标值以及对应时间值，该坐标值为二维 图片坐标系下的坐标，我们需将其转换为斜坐标系，之后将图片上的斜坐标转换为实际 情景下的坐标。 （1）平面二维坐标转换为斜坐标 在平面二维坐标系中，选定过杆底端的并与 x 轴在实际情景中垂直的线作为斜坐标 系 y 轴。y 轴与 x 轴在平面图片中存在夹角α，根据夹角可求得斜坐标系下坐标值 x 和 y。 图 11 画面中平面二维坐标系 图 12 斜坐标系 （2）斜坐标转换为实际情景直角坐标 转换后的实际情景直角坐标系: 图 13 实际情景直角坐标系 坐标值转换：   sin tan 斜 斜 斜 yy xyx   （12） 4.完成对视频中影子端点坐标的提取与转换，得到（ ii yx , ）。 建立优化模型：           Asky Askx LhNEyyLhNExxLhNEkk n i iii cos sin ])),,,,(()),,,,(()),,,,([(min 1 222  （13） 确定视频拍摄地点和日期。 5.4.2 模型的求解 根据问题四提供的视频图像，视频时长大概 40 分钟，我们每隔两分钟提取一帧画 面的数据，共提取 21 帧画面。 通过对图像中杆长像素点的计算 p=0.08，以及所给出的实际杆长 L=2m,得到像素点 与实际长度的比值 q=p/L=0.04,即每一个像素点对应实际长度 0.04 米。 根据视频图像内的信息显示（图 12），在当地时间上午 8:54:15，杆影在杆的右侧， 表明此时太阳在杆的左侧方向。 图 14 问题四视频截图 在上午 8:54:15 的时候，太阳位于杆的东侧，因此可判断图片中向左方向为东。 以 x 轴正向为东，y 轴正向为北，则视频中的坐标系方向如下图所示： 图 15 视频中杆所在坐标系 运用 matlab 软件分帧扫描视频，得到所取的 21 帧画面中影子端点坐标值和时间值如下： 表 2 视频画面中影子端点坐标值和时间值的变化 时间 8:54 8:56 8:58 9:00 9:02 9:04 9:06 9:08 9:10 9:12 9:14 坐标 x/m 0.3303 0.3383 0.3215 0.3199 0.3134 0.3085 0.3045 0.3005 0.2964 0.2899 0.2883 坐标 y/m 0.0195 0.0195 0.0195 0.0195 0.0173 0.0151 0.0130 0.0130 0.0108 0.0086 0.0086 时间 9:16 9:18 9:20 9:22 9:24 9:26 9:28 9:30 9:32 9:34 坐标 x/m 0.2835 0.2819 0.2770 0.2747 0.2699 0.2650 0.2609 0.2569 0.2529 0.2489 坐标 y/m 0.0086 0.0109 0.0065 0.0086 0.0086 0.0065 0.0043 0.0022 0.0022 0.0022 用 matlab 求解优化模型，得到可能地点： 纬度 ''52.4'5442 ，经度 ''62.48'17112 ，内蒙古自治区锡林郭勒盟苏尼特右旗 六、模型的评价与推广 6.1 模型评价 1.优点：充分利用地理知识，挖掘几何关系，得到经纬度、日期、时间、太阳高度角 等与影长的关系；分析影长与各参数的变化规律时，采用控制变量法，考虑单一参数与 影长的关系；问题二、三、四用最小二乘法拟合，得到最优解；问题四，在视频图像上 直接建立坐标系和斜坐标系。 2.缺点：拟合算法存在误差，拟合精度不高；问题四中提取坐标点时，由于像素点密 集，产生较大误差；计算角度时与弧度之间的转化会使精度降低，引起误差。 6.2 模型推广 1.结合太阳影子定位技术对视频拍摄地点和日期进行确定，可用于根据周围物体影子 定位的定位系统。 2.视频拍摄地点和日期确定模型可以应用于视频数据分析领域，进行视频信息分析。 3.采用更高精确度的拟合方法，将进一步优化本文模型，适用于更高标准的实际应用。 七、参考文献 [1]http://www.baike.com/wiki/太阳高度角 （2015-09-11） [2]http://3y.uu456.com/bp-e4dd3dcc4afe04a1b071deb7-1.html (2015-09-11) [3]http://map.yanue.net/ (2015-09-13) [4]天津大学，基于视频中太阳影子轨迹的经纬度估计方法，2009.7.29 [5]赵书兰，MATLAB 编程与最优化设计应用，北京：电子工业出版社，2013 [6]章绍辉，数学建模，科学出版社，2010 八、附录 程序代码： Code1 画出直杆的太阳影子的长度随时间变化的变化关系 H=9:0.05:15; h=H-4*(120-116.3914)/60;%地方时 a=((360/24)*abs(h-12))/180*pi;%时角 b1=39+54/60+26/3600;%北京纬度 b2=b1/180*pi;%纬度弧度 n=295;%2015 年 10 月 22 日的积日 n0=79.6764+0.2422*(2015-1985)-((2015-1985)/4); t=n-n0; d=(2*pi*t)/365.2422;%日角 e=(0.3723+23.2567*sin(d)+0.1149*sin(2*d)-0.1712*sin(3*d)-0.758*cos(d)+0.365 6*cos(2*d)+0.0201*cos(3*d))/180*pi; %赤纬角 %e=asin(0.3979*cos(0.9856*(n-173))); c=asin(sin(b2)*sin(e)+cos(b2)*cos(e)*cos(a));%太阳高度角 L=3./tan(c);%杆的影子长度 plot(H,L); xlabel('北京时间 H/h'); ylabel('影子长度 L/m'); title('影子长度随时间变化规律'); Code2 作图分析太阳影子长度随纬度的时间变化关系 t=9; k=3; p=116.4/180*pi; b=[-1.6:pi/100:1.6]; l=sqrt(((k.*sin((sin(sin(b).*(0.3979.*cos(0.9856.*(295-173)))+cos(b).*cos(( (360./24).*abs((t-4.*(120-p)./60)-12))./... 180.*pi)).*sin(b)-(0.3979.*cos(0.9856.*(295-173))))/cos(sin(b).*(0.3979.*co s(0.9856.*(295-173)))+cos(b).*cos((360./24)... .*abs((t-4.*(120-p)./60)-12))./180.*pi)).*cos(b))).^2)+((k.*cos((sin(si n(b)*(0.3979*cos(0.9856*(295-173)))+cos(b)... .*cos(((360./24).*abs((t-4.*(120-p)./60)-12))./180.*pi)).*sin(b)-(0.397 9*cos(0.9856*(295-173))))/cos(sin(b)... .*(0.3979.*cos(0.9856.*(295-173)))+cos(b).*cos(((360/24).*abs((t-4.*(12 0-p)./60)-12))./180*pi)).*cos(b))).^2); plot(b,l); Code3 读视频并提取每一帧,处理图像得到影子端点坐标值 x=[]; y=[]; fileName ='D:Program FilesmatlabworkAppendix4.avi'; obj = VideoReader(fileName); numFrames = obj.NumberOfFrames;% 帧的总数 t=obj.Duration %视频的总时长（秒） v=obj.FrameRate % 视频帧速（帧/秒） frame1 = read(obj,1); imshow(frame1);%显示帧 [xa,ya]=ginput(1);%提取原点 for k = 1:3000: numFrames% 读取数据 frame = read(obj,k); frame1=im2bw(frame,0.768);%二值图像 imshow(frame1);%显示帧 imshow(frame);%显示帧 [x(k),y(k)]=ginput(1); imwrite(frame,strcat(num2str(k),'.bmp'),'bmp');% 保存帧 end Code4 处理附件 3，给出可能的地点和日期 Title "第三问附件 3"; Parameters p[-2.0,2.0],b[-1.6,1.6],k,n[0,366]; Variable t,l; Function l=sqrt(((k*sin((sin(sin(b)*(0.3979*cos(0.9856*(n-173)))+cos(b)*cos(((360/24) *abs((t-4*(120-p)/60)-12))/180*pi))*sin(b)-(0.3979*cos(0.9856*(n-173))))/co s(sin(b)*(0.3979*cos(0.9856*(n-173)))+cos(b)*cos((360/24)*abs((t-4*(120-p)/ 60)-12))/180*pi))*cos(b)))^2)+((k*cos((sin(sin(b)*(0.3979*cos(0.9856*(n-173) ))+cos(b)*cos(((360/24)*abs((t-4*(120-p)/60)-12))/180*pi))*sin(b)-(0.3979*c os(0.9856*(n-173))))/cos(sin(b)*(0.3979*cos(0.9856*(n-173)))+cos(b)*cos(((3 60/24)*abs((t-4*(120-p)/60)-12))/180*pi))*cos(b)))^2); Data; t l 13.15 3.533142184 13.2 3.546768029 13.25 3.561797643 13.3 3.578100715 13.35 3.595750783 13.4 3.61493428 13.45 3.635425983 13.5 3.657218272 13.55 3.680541115 13.6 3.705167836 13.65 3.731278025 13.7 3.758917911 13.75 3.788087888 13.8 3.818701015 13.85 3.850809619 13.9 3.88458522 13.95 3.919911828 14 3.956875992 14.05 3.99553479 14.1 4.035750835 14.15 4.077863059 Code5 计算可能的拍摄地点 Title "第四问 11"; Parameters p[-2.0,2.0],b[-1.6,1.6],k; Variable t,l; Function l=sqrt(((k*sin((sin(sin(b)*(0.3979*cos(0.9856*(194-173)))+cos(b)*cos(((360/ 24)*abs((t-4*(120-p)/60)-12))/180*pi))*sin(b)-(0.3979*cos(0.9856*(194-173))) )/cos(sin(b)*(0.3979*cos(0.9856*(194-173)))+cos(b)*cos((360/24)*abs((t-4*(1 20-p)/60)-12))/180*pi))*cos(b)))^2)+((k*cos((sin(sin(b)*(0.3979*cos(0.9856* (194-173)))+cos(b)*cos(((360/24)*abs((t-4*(120-p)/60)-12))/180*pi))*sin(b)- (0.3979*cos(0.9856*(194-173))))/cos(sin(b)*(0.3979*cos(0.9856*(194-173)))+c os(b)*cos(((360/24)*abs((t-4*(120-p)/60)-12))/180*pi))*cos(b)))^2); Data; t l 8.9 0.330858465 8.93 0.338844942 8.97 0.322074114 9 0.320477052 9.03 0.313886486 9.07 0.308904338 9.1 0.304733408 9.13 0.300737082 9.17 0.296578055 9.2 0.29003374 9.23 0.288434455 9.27 0.28363665 9.3 0.282188113 9.33 0.2771045 9.37 0.274840873 9.4 0.270043292 9.43 0.265107937 9.47 0.260988208 9.5 0.25688529 9.53 0.252885434 9.57 0.248885582 Code6 计算可能的拍摄时间和日期 Title "第四问 22"; Parameters p[-2.0,2.0],b[-1.6,1.6],k,n[1,366]; Variable t,l; Function l=sqrt(((k*sin((sin(sin(b)*(0.3979*cos(0.9856*(194-173)))+cos(b)*cos(((360/ 24)*abs((t-4*(120-p)/60)-12))/180*pi))*sin(b)-(0.3979*cos(0.9856*(n-173)))) /cos(sin(b)*(0.3979*cos(0.9856*(n-173)))+cos(b)*cos((360/24)*abs((t-4*(120- p)/60)-12))/180*pi))*cos(b)))^2)+((k*cos((sin(sin(b)*(0.3979*cos(0.9856*(n- 173)))+cos(b)*cos(((360/24)*abs((t-4*(120-p)/60)-12))/180*pi))*sin(b)-(0.39 79*cos(0.9856*(n-173))))/cos(sin(b)*(0.3979*cos(0.9856*(n-173)))+cos(b)*cos (((360/24)*abs((t-4*(120-p)/60)-12))/180*pi))*cos(b)))^2); Data; t l 8.9 0.330858465 8.93 0.338844942 8.97 0.322074114 9 0.320477052 9.03 0.313886486 9.07 0.308904338 9.1 0.304733408 9.13 0.300737082 9.17 0.296578055 9.2 0.29003374 9.23 0.288434455 9.27 0.28363665 9.3 0.282188113 9.33 0.2771045 9.37 0.274840873 9.4 0.270043292 9.43 0.265107937 9.47 0.260988208 9.5 0.25688529 9.53 0.252885434 9.57 0.248885582