2015全国数学建模竞赛A题一等奖

1 太阳影子定位太阳影子定位 摘要 视频拍摄地点和拍摄日期的确定在视频数据分析中占据着重要的位置， 本文 建立了影子长度与日期、时间、经纬度、直杆高度之间的函数关系模型，通过最 小二乘估计方法和小孔成像原理来确定视频拍摄的地点和日期。 针对问题一，根据太阳与地球之间的物理关系，建立了影子长度关于日期、 时间、经度、纬度、杆长 5 个变量的非线性函数模型一，并讨论了影子长度与这 5 个变量变化的关系，揭示了 5 个变量变化对影子长度的影响。把所建立模型应 用到问题一中的例子，给出了影子长度关于北京时间 9:00-15:00 的函数关系，并 描绘了影子长度与时间变化曲线图像。 针对问题二， 首先对太阳影子顶点坐标数据进行处理计算出每个时刻影子长 度；其次对模型一进行改进建立了影子长度关于日期、时间、经度、纬度、杆长 5 个变量的非线性函数模型二；然后根据数据（影子长度和时间）采用最小二乘 方法对模型二中未知参数 （经度、 纬度） 进行估计， 获得了若干个可能的地点 （经 度、纬度） ；最后讨论了问题二中附件一影子顶点坐标轴的方向，指明了坐标轴 确定方向的思路。 针对问题三， 首先对附件中太阳影子顶点坐标数据进行处理计算出每个时刻 影子长度，其次根据数据与模型二，类似采用最小二乘方法对模型中的未知参数 （经度、纬度、日期）进行估计，获得了若干个可能的地点（经度、纬度）和日 期。 针对问题四，引进“虚影子长度”变量，根据小孔成像原理和视频数据得到 “虚影子长度”变量各个时刻的大小，并建立了“虚影子长度”关于日期、时间、 经度、纬度、杆长、拍摄位置( , ) 7 个参数的非线性函数模型，并采用最小二 乘方法对模型中未知参数进行估计，得到杆的位置和若干个可能拍摄的地点。并 论证了，如果拍摄日期未知，能根据视频确定出拍摄地点与日期。 关键词：最小二乘估计法；非线性拟合；虚影子 2 一、问题重述一、问题重述 如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面， 太阳影子 定位技术就是通过分析视频中物体的太阳影子变化， 确定视频拍摄的地点和日期 的一种方法。 1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律， 并应用你们建立的模型画出 2015 年 10 月 22 日北京时间 9:00-15:00 之间天安门 广场（北纬 39 度 54 分 26 秒,东经 116 度 23 分 29 秒）3 米高的直杆的太阳影子 长度的变化曲线。 2.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确 定直杆所处的地点。将你们的模型应用于附件 1 的影子顶点坐标数据，给出若干 个可能的地点。 3. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型 确定直杆所处的地点和日期。 将你们的模型分别应用于附件 2 和附件 3 的影子顶 点坐标数据，给出若干个可能的地点与日期。 4附件 4 为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估 计出直杆的高度为 2 米。请建立确定视频拍摄地点的数学模型，并应用你们的模 型给出若干个可能的拍摄地点。 如果拍摄日期未知，你能否根据视频确定出拍摄地点与日期？ 二、问题分析二、问题分析 问题一，由于太阳高度角与直杆高度、影子长度存在三角关系，从这个突破 口，我们可以寻找出日期、时间、经纬度、直杆高度、影子长度之间的关系来建 立模型。 通过改变一个参数， 确定其它参数来分析影子长度与各参数之间的关系。 根据问题中给定的参数值就可以求解出直杆的太阳影子长度的变化曲线。 问题二，因为问题中所涉及的参数与问题一是一样的，通过问题一的模型， 利用已知参数，可以将对未知参数的求解问题转化成最小二乘估计问题，并确定 当满足条件时， 所得参数值即为所求。 关于坐标系的问题， 先假定坐标系的方向， 计算出方位角，根据方位角与附件 1 中数值的关系，旋转坐标系使得参数与附件 3 1 中数值一致，此时的坐标系方向即为附件 1 中坐标系的方向。 问题三，问题三是在问题二的基础上，待求参数增加了一个日期，所以可以 依照问题二的方法求解。 问题四，第一小问分析：首先，对视频的文件进行读取，采用 matlab 读取 视频 VideoReader 类，将读取出不同时刻直杆在摄像机中成像的直杆像长度 h ， 影子像长度 1 i l ，为了能够确定视频的拍摄位置，引入方位角，建立极坐标系。 其次， 因为摄像机的焦距的大小相对摄像机离杆的距离可忽略不计，所以在这里 为了计算简便，把摄像机凸透镜的成像原理转换成小孔成像原理，分别绘制直杆 的小孔成像图、影子的小孔成像图。接着，根据几何关系得出各变量的函数关系 表达式，最后，利用最小二乘优化方法得出视频拍摄地点的数学模型。第二小问 分析：问题“如果拍摄日期未知，能否根据视频确定出拍摄地点与日期？”经查 阅百度百科， 所谓的 “日期” 是指约定的日子和时间， 日期的标准格式为年月日。 对日期的定义比较模糊，有时日期也可以认为指年月日时分。那么针对本问题， 我们将分成两中情况。第一种：日期指月日（不考虑平年闰年） 、第二种：日期 指月日时分。 三、符号说明三、符号说明 N从元旦到计算日的总天数 h杆长，其单位是米 l影长，其单位是米 0 t北京时间，其单位是时，分 t地方时，其单位是时，分 d日期，其单位是年，月，日 0 北京所在的经度，其单位是度 直杆所在位置的经度，其单位是度 直杆所在位置的纬度，其单位是度 赤纬角，其单位是度 时角，其单位是度 4 高度角，其单位是度 方位角，其单位是度 p E标准时差，其单位是分 关于问题四中的符号说明在问题四中呈现。 四、模型假设四、模型假设 1、地球是一个规则的球体。 2、地球的自转速度是一个常数。 3、地球的公转速度也是一个常数。 4、不考虑闰年和平年。 5、东西经度统一为经度，南北纬度统一为纬度。 6、杆所在的地面是水平的。 7、直杆始终垂直于地面。 8、相机拍摄高度相对于地面为零。 五、问题一五、问题一 1.影长求解模型的建立 直杆的影子是由于太阳对直杆的照射产生的， 太阳光照射到直杆影子顶点的 入射方向和地平面之间的夹角实际上就是直杆高度和影子长度的关系， 即太阳高 度角与直杆高度、影子长度的关系tan( )lh的原理图，如图 1.1。其中，太阳 高度角又与赤纬角、时角等有关，通过对一系类关系的整理，建立了影长求解模 型。 5 h l 图 1.1 原理图 1.1 通过相关理论基础对赤纬角、时角、高度角进行解释 1.1.1 赤纬角 1 太阳赤纬角指地球赤道平面与太阳和地球中心的连线之间的夹角， 建筑设 计资料集中给出了赤纬角近似公式:例如：大寒日常 1 月 20 日，N=20，则 959120 。 )9500/1 (*)25.80sin(*45.23NN 例如：大寒日常 1 月 20 日，N=20，则959120 。 1.1.2 时角 2 单位时间地球自转的角度定义为时角， 规定正午时角为0， 上午时角为负值， 下午时角为正值。地球自转一周360，对应的时间为 24 小时，即每小时相应的 时角为15，每 4 分钟的时角为1。时角的计算公式是： 15 (t 12) 给定精确的地方时与标准时之间的转换关系： 00 4() p ttE 6 其中，系数 4 经度上每一度对应的是四分钟。 因为时差 p E的值对模型求解的影响不大， 为了模型求解的简便， 修正值 p E 可以忽略不计，公式可以简化为： 00 4()tt。 1.1.3 太阳高度角 太阳高度角是指太阳光的入射方向和地平面之间的夹角， 专业上讲太阳高度 角是指某地太阳光线与通过该地与地心相连的地表切线的夹角。 太阳高度角简称 高度角。太阳高度角的公式为： sin( )sin( )sin( )cos( )cos( )cos( ) 结合上面几个定义，整理出了以下几个关系式： 00 tan( ) sin( )sin( )sin( )cos( )cos( )cos( ) 23.45sin(80.25) (19500) 15(t 12) 4 () lh NN tt 由关系式，建立了影长关于参数日期、时间、直杆纬度、直杆经度、直杆高 度的模型一模型一： 2.各参数和影长关系的分析 将影长求解模型一简写为： 0 ( , , )lf hN t。其中参数有 5 个， 直杆高度h、 直杆纬度、日期d、时间 0 t、直杆经度。 为了要分析影长和各参数间的单独关系，我们可以给定其中 4 个参数的值， 来分析另外 1 个参数变化对影长的关系。不妨假设 5 个参数值（见表 1.1） ，从中 选择 4 个来确定分析。 表 1.15 个参数值 直杆高度h5 北京经度 0 116.3914 日期d3 月 10 日（N=70）直杆经度100.3710 时间 0 t9直杆纬度 39.9072 7 当其它 4 个参数给定时，分别作出影子与日期、时间、当地纬度、直杆长度的关 系图 1.2、1.3、1.4、1.5，如下图所示： 图 1.2图 1.3 图 1.4图 1.5 由以上 4 幅图可知，当其它参数给定时，影子长度与直杆高度成线性关系， 与日期、时间和纬度成非线性关系，并且与后三者存在周期性。 其中，影子长度与日期成年周期性，近日点最短，远日点最长；同样的，与 时间成日周期性，正午影子长度最短，并向两边逐渐变长；影子长度关于纬度的 变化规律也是成年周期性的，并且影子的长度随着纬度的增加成指数增长，在这 里因为给定的日期是 3 月 10 日，此时的太阳位置在南半球，而直杆位置在北半 球，太阳的照射范围是 90 度，所以当纬度达到 84 度左右时，就已经接近太阳的 照射尽头了， 这个时候的影子长度就接近最长， 当直杆纬度超过太阳照射范围时， 影子长度的值就为负值了。 8 3.问题求解 为了确定 2015 年 10 月 22 日北京时间 9:00-15:00 之间天安门广场（北纬 39 度 54 分 26 秒,东经 116 度 23 分 29 秒）3 米高的直杆的太阳影子长度的变化 曲线。在影长求解模型一中， 0 是北京经度，是直杆的经度，此时的直杆位置 在北京，即 0 。模型一可以简化为： / tan(arcsin(sin( )sin(23.45sin(80.25) (19500) cos( )cos(23.45sin(80.25) (19500)cos(15(t 12) lhNN NN 此时，日期为 10 月 22 日即295N ，39.9072 ，3h ，把这些值代入 后，l关于t就是一元函数。 因此可以计算不同时刻下影子的长度，见表 1.2。 根 据影长与时间变化数据拟合出直杆的太阳影子长度的变化曲线如图 1.6。 表 1.2 不同时刻下影子的长度 t 99.51010.51111.512 l6.68285.56584.81024.29373.95443.76063.6975 t 12.51313.51414.515 l3.76063.95444.29374.81025.56586.6828 图 1.6 太阳影子长度的变化曲线 9 六、问题二六、问题二 1.位置求解模型二的建立 问题二中，要求根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建 立数学模型确定直杆所处的地点。 根据太阳影子顶点坐标数据可以计算出每个时 刻对应的影长，如图 2.1 是附件 1 中影长与时间变化关系的散点图。 图 2.1 影长与时间的散点图 由影子长度与各参数的函数关系式 0 ( , , )lf hN t可知，影长与纬度、经 度、日期、时间等的关系为非线性。根据题目附件 1 数据知道108N （4 月 18 日） ，以及北京经度 0 =116.3914。影长函数 00 ( , , )lf hN t 仅由参数 0 , ,ht确定，而时间 0 t与影长l的数据可以根据题目附件 1 影子顶点坐标数据 获得（ 0, t l）见表 2.1。该问题就转化成函数模型形式 00 ( , , )lf hN t ，和已 知观测数据（ 0, t l） ，未知参数, ,h 的问题。这个问题的未知参数就可以采用最 小二乘估计来求解。 最小二乘估计 3 的目的就是估计参数, ,h 使得函数 0, f t在点 0i t1,2,in处的函数值与观测数据偏差的平方和达到最小。 即求满足如下条 10 件的函数 0 ,f t使得： 2 2 00 11 min, n iiii ii f tlf tl 其中 , , h 是待定的参数（即：纬度、经度、直杆的高度） ，而 就是 最小二乘法所确定的最佳参数，从而确定了直杆所处的地理位置（纬度和经度） 。 2.问题求解 由附件 1 将影子顶点坐标数据转换为影子的长度l，将北京时间以小时的形 式转换为数值 0 t ，得到表 2.1 如下。 表 2.1 不同时间对应的影长 时间 0 t影长l 时间 0 t影长l 14.71251.540231817 14.75131.579853316 14.81.21529695515.351.620144515 14.851.24905105215.41.661270613 14.91.2831953415.451.703290633 14.951.31799314915.51.74620591 151.35336404915.551.790050915 15.051.38938709115.61.835014272 15.11.42615285615.651.880875001 15.151.46339985315.71.927918447 15.21.501481622 并且直杆影子顶点坐标数据的测量日期为： 2015 年 4 月 18 日， 故有：N=108， 因此求得赤纬角 8042.10。 根据客观事实可以确定直杆在 618 时才有影子，故直杆位置与北京的经度 差约在0至120之间。根据直杆的影长给出直杆高度h的大约范围在 0m5m，在 用 MATLAB 求解时，需要给定待估参数的初始值，如表 2.2。 11 表 2.2 待估参数初始值 直杆的纬度直杆的经度 直杆的高度h 初始值 1 251103m 初始值 2 45704m 初始值 3 70265m 利用模型二求解出各参数值，见表 2.3。 表 2.3 拟合结果 直杆的纬度直杆的经度 直杆的高度h误差平方和 结果值 1 19.3237105.22852.0294m3.099408e 结果值 2 19.3237105.22852.0294m3.099408e 结果值 3 19.3237105.22852.0294m3.099408e 由表 2.3 中误差平方和可见通过拟合的方法求解直杆的所在的纬度、经度、 直杆高的准确度很好，直杆的纬度是19.3237，经度是105.2285，所以可能的地 点有四个。 3.坐标方向的确定方法 假设坐标系的 x0 轴方向为正东方向，y0 轴方向为正北方向。根据太阳方 位角的计算公式： cos/sincossin 由表 2.3 的数据，根据模型一求解出方位角的数值，判断附件 1 中， arctan( / )y x 与 90 是否相等？ 如果不相等，求解出arctan( / )y x 与180 的差值，当差值为正时，将假设 12 的坐标系按逆时针方向旋转等于差值的度数；当差值为负时，将假设的坐标系按 顺时针方向旋转等于差值的度数。旋转后的坐标系方向就是附件 1 中所给定的。 七、问题三七、问题三 通过对问题三的分析，可以知道问题三是在问题二的基础上，将日期从给定 的参数转变为求解的参数。即在模型二中，待定参数里增加了日期变量N，所 以依然可以采用模型二的方法来求解问题三。 将附件 2、3 中影子顶点坐标和北京时间转换成时间与影长的数据如表 3.1 所示。 表 3.1 时间与影长的关系 附件 2附件 3 时间影长时间影长 12.681.24725620513.153.533142184 12.731.2227945913.23.546768029 12.781.19892148613.253.561797643 12.83133.578100715 12.881353.595750783 12.931.1299174713.43.61493428 12.981.1078354813.453.635425983 13.031.08625420613.53.657218272 13.081.06508107213.553.680541115 13.131.04444626513.63.705167836 13.181.02426412613.653.731278025 13.231.00464031413.73.758917911 13.280.98549090813.753.788087888 13.330.96679049413.83.818701015 13.380.94858473513.853.850809619 13.430.93092788113.93.88458522 13.480.9137517513.953.919911828 13.530.897109051143.956875992 13.580.88097376214.053.99553479 13.630.86549225914.14.035750835 13.680.85050446814.154.077863059 给定参考值表 3.2。 13 表 3.2 待估参数初始值 直杆的纬度直杆的经度 直杆的高度h N 值 初始值 1 371103m 100 初始值 2 15603m 268 初始值 3 651325m 53 利用模型二求解出各参数，如表 3.3。 表 3.3 第三问的结果 附件 2 直杆的纬度直杆的经度 日期d直杆的高度h误差平方和 位置 1 39.892676.1852 5 月 23 日2.0008m1.682808e 位置 2 39.865276.1762 7 月 19 日1.9992m1.684008e 位置 3 39.892676.1852 5 月 23 日2.0008m1.682808e 附件 3 直杆的纬度直杆的经度 日期d直杆的高度h误差平方和 位置 1 32.9965106.6165 1 月 7 日3.0567m2.975408e 位置 2 32.8488106.6364 11 月 6 日3.0356m2.840108e 位置 3 31.2946107.0438 1 月 14 日2.6438m5.005207e 由表 3.3 中的误差平方和，可见通过拟合的方法求解直杆的所在的纬度、经 度、直杆高和日期有较高的准确性。 14 八、问题四八、问题四 首先，对视频的文件进行读取，采用 Matlab 读取视频 VideoReader 类，将读 取出不同时刻直杆在摄像机中成像的直杆像长度 h ，影子像长度 1 i l ，为了能够 确定视频的拍摄位置，引入方位角，建立极坐标系。其次，因为摄像机的焦距 的大小相对摄像机离杆的距离可忽略不计，所以在这里为了计算简便，把摄像机 凸透镜的成像原理 4 转换成小孔成像原理，分别绘制直杆的小孔成像图、影子的 小孔成像图。接着，根据几何关系得出各变量的函数关系表达式，最后，利用最 小二乘优化方法得出视频拍摄地点的数学模型。 1. 相关说明 （1）假设摄像机在水平地面拍摄。 （2）方位角:太阳光线在地面上的投影线与地面正南方向所夹的角， 180180,且cos/sincossin。考虑到对称性，本为仅讨论范围 在 1800。方位角在下文指的就是影子与 X 轴的夹角。 （3）极坐标系的建立：以正南方向作为 X 轴，以正东方向为 Y 轴，杆底作为原 点。 （4） 摄像机的拍摄坐标：先将摄像机视为一个质点O，则坐标为：,，表 示摄像机到杆底部的距离，表示摄像机位于杆南偏东的位置。 （5） 小孔成像 5 ：根据光的直线传播原理，用一个带有小孔的板遮挡在屏幕与 物之间，屏幕上就会形成物的倒像，把这样的现象叫小孔成像。如图 4.1。 图 4.1 小孔成像 光屏 小 孔 15 图 4.2 直杆影子成像的三维坐标图 （6） 将视频中的影子变化的动态图， 在三维坐标系的静态图大致如下， 这里以 为钝角为例，因为随着时间变化直杆像长度 h 始终不变，影子像长度 1 i l 随之变 化。于是将摄像机中的成像图分成直杆小孔成像平面图和影子小孔成像平面图。 这里绘制杆影子在摄像机中的成像静态三维图大致如下， 如图 4.2。 以下是图 4.2 中的符号说明： A 表示直杆； O 表示摄像机的拍摄位置，,； BC 表示真实影子， i lBC ； 是线段 BC 与线段 BD 的夹角； 方位角，BC 与 x 轴的夹角； 是 OB 与 x 轴的夹角； DB表示影子像， 1 i lDB； 光屏位置：过线段DB垂直于 OB 的平面； BD 表示虚影子，因为BOBD 根据光的直线传播原理，BC 所成的像可以认为是 同一时刻在同一速光线照射 BD 所成的像， 于是定义 BD 是 BC 的虚影子， 1 i lBD ； O 摄像机 A h i l BC 东（y） 南（x） D 光屏 B D 1 i l A h 16 2. 利用 Matlab 读取视频 VideoReader 类，提取 h ， 1 i l 对视频数据提取的主要步骤如下： Matlab 读取视频数据，帧的总数 n。 以 a 为步长，提取出 n/a 帧图像，将其表示为),.2 , 1(),( a n iif。 将)(if转换成灰度图)(i f 。 对图像进行腐蚀，去除图像的噪声。 用 PS 软件对杆顶、杆底、影子定点进行定位得出三点坐标，从而得出直杆像 的高度 h 、影子像的长度 1 i l 。 将视频每隔两分钟进行信息提取，共提取到 21 张图片，并通过 PS 软件读取到视 频中影子顶点、直杆底部和直杆顶部点的坐标数据如表 4.1、表 4.2. 表 4.1 影子像顶点坐标数据 时间XY时间XY 8.9059.1130.549.2755.5631.01 8.9358.9130.629.3055.3231.01 8.9758.4930.669.3354.7931.04 9.0058.0630.659.3754.4731.08 9.0357.7430.709.4054.1931.08 9.0757.5730.739.4353.6631.15 9.1057.1530.809.4735.4831.15 9.1356.6930.839.5053.2331.19 9.1756.3730.879.5352.9531.19 9.2056.1630.909.5752.6331.22 9.2355.8130.94 表 4.2 直杆底部和顶部坐标 直杆底部坐标直杆顶部的坐标 X30.6131.46 Y30.757.13 直杆底部和直杆顶部的横坐标不一致的原因：如图 4.3 直杆的截图，可知由 于杆的衔接位子错开，导致直杆的底部与顶部的 X 坐标有偏差。通过直杆顶部 Y 值和直杆底部 Y 值可以计算出视屏中直杆的高度为：23.62。 17 图 4.3 直杆截图 3. 小孔成像图 3.1 直杆（不考虑影子）的小孔成像平面图，如图 4.4。 图 4.4 杆的小孔成像 由于ABO与OBA相似，相似比为hh/，用 K 表示。 3.2 影子的小孔成像平面图 因为影子与正南方向的夹角（方位角）分别为锐角和钝角时，形成的几何图 形不同，我们猜测构造的函数表达式可能也有所不同。所示在构造函数表达式之 前先对方位角的大小是否影响函数表达式进行检验。 于是把影子的小孔成像图像 分成影子小孔成像平面图一、二。如图 4.5、4.7。 O A B A B h h 摄像机 （2 米） 18 当18090时： 、 图 4.5 影子小孔成像平面图（一） 几何关系： OBD与DBO相似，相似比=K 11jj kll 在OBDRt中： 1 tan j l O Oljtan 1 在OBC中，OBC,作OC边的垂线，如图 4.6 图 4.6OBC二维图（一） 通过三角函数关系得到： OOC180)(180 OOljsin)180sin( 整理得到： 东（y） A 南（x） B C O 摄像机 光屏 B D 1 j l D j l 1 j l B C O j l 19 )cos( )sin( tan j j l l O 当 900时： 图 4.7 影子小孔成像平面图（二） 几何关系： OBD与DBO相似，相似比为K 11ii kll 在OBDRt中： 1 tan i l O Olitan 1 在OBC中，OBC,作 OC 边的垂线，如图 4.8 图 4.8 OBC 二维图（二） 通过三角函数关系得到： OE OOlisin)sin( 东（y） B 南（x） CD 光屏 O i l 1 i l 摄像机 B D 1 i l A B CO i l E 20 整理得到与当 18090时一样的关系： )cos( )sin( tan i i l l O 由此可以认为不管为锐角还是钝角时，得到最终的函数表达式都是一致的。所 以影长的下标不需要区分ji,，统一用i表示。 4. 视频拍摄地点模型的建立 由于 )cos( )sin( tan i i l l O（1） 由方位角的表达式： cos/sincossin， 得到： )cos/sinarcsin(cos（2） 又因为影长与杆的高度关系式为： tan/hli （3） 将（2） （3）代入（1）式中，得到： )cos/sin(coscos(arcsin)tan/( )cos/sin(cossin(arcsin)tan/( tan h h O 再由：Olitan 1 得到视频拍摄地点的数学模型视频拍摄地点的数学模型： )cos/sin(coscos(arcsin)tan/( )cos/sin(cossin(arcsin)tan/( 1 h h li 其中： )cos()cos()cos()sin()sin()sin( )(4 )12(15 )9500/1 ()25.80sin(45.23 00 tt t NN 最后得到虚影子与各参数的函数表达式),( 0011 tNhlfl ii 。其中 1 i l 是通 21 过图 2.2 中相似关系得到，即： DB BA BD AB ,即 11ii l h l h h hl l i i 1 1 。 最后， 通过部分已知数据用曲线拟合的线性最小二乘法将模型中的待定系数求解 出来。 5.模型求解 5.1 第一小问模型求解 对于第一小问模型中未知参数是.,已知数据是., 001 tNhli用曲线拟合 的线性最小二乘法将模型中的待定系数求解出来。 最后给 3 个出参数的初始值得 到视频的 3 个可能拍摄地点。如下表 4.3、表 4.4 所示。 表 4.3 待估参数初始值 直杆的纬度直杆的经度 初始值 1 251103m2rad 初始值 2 45706m3rad 初始值 3 501209m4rad 表 4.4 拍摄地点的结果值 直杆的纬度直杆的经度 误差平方和 位置 1 29.6380107.59202.7874m3.2413rad3.18994e 位置 2 25.8220109.47161.4830m3.3049rad3.10274e 位置 3 80.3215105.44268.2649m2.9830rad3.10704e 由三组初始值得到的结果误差平方和都很小，又直杆的经度相对稳定，于是 可以认为得出的结果是较为准确的。对于拍摄地点的表述以初始值 1 为例：由第 22 一组初始值得到，直杆的位于地球纬度为29.6380，经度为107.5920，摄像机位 于直杆底部南偏东3.2413rad，相距2.7874m。 5.2 第二小问的模型建立、求解 对于第二小一问“如果拍摄日期未知，能否根据视频确定出拍摄地点与日 期？”经查阅百度百科，所谓的“日期”是指约定的日子和时间，日期的标准格 式为年月日。对日期的定义比较模糊，有时日期也可以认为指年月日时分。那么 针对本问题，我们将分成两中情况。第一种：日期指月日（不考虑平年闰年）、 第二种：日期指月日时分。 第一种情况：月日未知，利用视频拍摄地点的数学模型： )cos/sin(coscos(arcsin)tan/( )cos/sin(cossin(arcsin)tan/( 1 h h li 其中： )cos()cos()cos()sin()sin()sin( )(4 )12(15 )9500/1 ()25.80sin(45.23 00 tt t NN 得到虚影子与各参数的函数表达式),( 0011 tNhlfl ii 。 此时在该模型中未知参数是.,N已知数据是., 001 thli 任然通过曲线拟合的线性最小二乘法将模型中的待定系数求解出来。 最后得出的初始值和结果值分别如表 4.5、表 4.6 所示。 表 4.5 月日未知的参数初始值 直杆的纬度直杆的经度 N 初始值 1 251103m2rad90 初始值 2 45706m3rad200 初始值 3 501209m4rad160 23 表 4.6 月日未知的结果值 直杆的纬度直杆的经度 N误差平方和 位置 1 28.6494101.85675.5899m2.9336rad1353.3995e-04 位置 2 26.783999.57783.2937m2.9547rad1263.3437e-04 位置 3 70.9661120.888817.4648m3.2670rad1323.1092e-04 第二种情况：月日时分都未知： 当月日时分都未知时。我们需要对时分进行推算，当找出最适合的时间时接 着就可以用第一种情况进行求解。 那么通过视频中的信息推测当天的拍摄大致的 时间区间，最终得出拍摄位置的具体步骤如下： 绘制视频中的影子长度与时间的关系图，称其为图 a。 再绘制实际影子长度与时间的关系式，其中时间区间为 6:00-18:00，称 其为图 b。 我们由实际经验也可知道， 该曲线图近似为开口向上的二次函数图象。 由图 a 中曲线的变化规律可以拍摄的大致时间。即：曲线成上升趋势，可 以认为是在上午拍摄。同理，曲线成下降趋势，则认为是在下午拍摄。这里假设 为上午拍摄。 那么对拍摄的时间的可能取值范围缩短到了 6:00-12:00。 接着把该区间进行分段：6:00-7:00、6:01-7:01、 6:02-7:02.10:59-11:59、11:00-12:00。 将这些时间段代入视频拍摄地点的数学模型中， 比较拟合残差平方和， 选 取拟合残差平方和最小的模型。 比较拟合残差平方和， 选 取拟合残差平方和最小的模型。 此时残差平方和最小的模型残差平方和最小的模型所对应求出的日期、 时间、 地点就是视频拍摄 的日期、时间、地点。 综上所述，如果拍摄日期未知，我们能能根据视频确定出拍摄地点与日期。 九、结语九、结语 根据模型一，首先分析了影子长度与各个参数的关系，并画出了北京时间与 影子长度的关系图；其次对模型一进行改进并建立了模型二，采用最小二乘方法 24 估计出了若干个可能的地点，并且讨论了附件中影子顶点坐标轴的方向，由于时 间的关系，仅指明了坐标轴确定方向的思路；然后针对问题三，利用模型二，类 似采用最小二乘方法， 获得了若干个可能的地点和日期； 最后通过巧妙的引进 “虚 影子长度”变量，根据小孔成像原理和视频数据求出“虚影子长度”变量各个时 刻的大小以及建立的模型三， 采用最小二乘方法估计出了杆的位置和若干个可能 拍摄的地点，而且详细论证了，如果拍摄日期未知，我们能根据视频确定出拍摄 地点与日期，但是由于时间关系，没能完整的证明。 十、参考文献十、参考文献 1 杜春旭，王普，马重芳，吴玉庭，申少青，一种高精度太阳位置算法，新能 源及工艺， 卷起号：第二期，41-48，2010 年。 2 陈登鳌等， 建筑设计资料集 ，出版地：北京，中国建筑工业出版社，1994 年。 3 ccx19851205， 第七讲matlab实现非线性拟合， mQsVxr4aJXFVkMQNU5GyZUJahX2GC_t7FgP69Did9IJGvsJhUDz2ilSJggfbUQh qNxqVfDNZg_z YugyqjTu8eZ-YZyUty0yV5tL7OA2ccgm，2015 年 9 月 12 日。 4 百度百科，凸透镜成像规律， RHFu9_k9zUgxdPy5Faosl3CGkqRzFaeGxXuiWaHGK7DgSL9PtD9PHdTsVmNAZ OY5PN8SvTrxa l = ; for i = 1:length(t) omega = 15*(t(i) - 12); alpha = asin(sin(varphi*pi/180)*sin(delta*pi/180)+cos(varphi*pi/180)*cos(delt a*pi/180)*cos(omega*pi/180); temp = h/tan(alpha); l = l, temp; end plot(t, l) %t,l end 问题二程序： 根据附件 1 中影子顶点数据，求解若干个可能的地点 pro2 t0 = 14.714.7514.814.8514.914.951515.0515.1 15.1515.215.2515.315.3515.415.4515.515.55 15.615.6515.7; l = 1.149625826 1.1821989761.215296955 1.249051052 1.28319534 1.317993149 1.353364049 1.389387091 1.426152856 1.463399853 1.501481622 1.540231817 1.579853316 1.620144515 1.661270613 1.703290633 1.74620591 1.790050915 1.835014272 1.880875001 1.927918447; b0 = 70,26,5; delta = 10.8042; fun = inline(b(3)./tan(asin(sin(b(1)*pi/180)*sin(10.8042*pi/180)+cos(b(1)* pi/180)*cos(10.8042*pi/180)*cos(15*(t0 - 4*(116.3914-b(2)/60 -12)*pi/180),b,t0); b,r,j = nlinfit(t0, l, fun, b0) R = sum(r.2) 26 问题三程序： 根据附件 2 中的影子顶点数据，求解若干个可能得地点。 pro3_1 t0 = 12.68 12.7312.7812.8312.8812.9312.9813.0313.08 13.1313.1813.2313.2813.3313.3813.4313.4813.53 13.5813.6313.68; l = 1.2472562051.222794591.198921486 1.175428964 1.152439573 1.129917471.107835481.086254206 1.065081072 1.044446265 1.024264126 1.004640314 0.985490908 0.966790494 0.948584735 0.930927881 0.91375175 0.897109051 0.880973762 0.865492259 0.850504468; b0 = 37,110, 3, 100; fun = inline(b(3)./tan(asin(sin(b(1)*pi/180)*sin(23.45*sin(b(4) - 80.25)*(1 - b(4)/9500)*pi/180)*pi/180)+cos(b(1)*pi/180)*cos(23.45*sin(b(4) - 80.25)*(1 - b(4)/9500)*pi/180)*pi/180)*cos(15*(t0 - 4*(116.3914-b(2)/60 -12)*pi/180),b,t0); b,r,j = nlinfit(t0, l, fun, b0) R = sum(r.2) 根据附件 3 中的影子顶点数据，求解若干个可能得地点。 pro3_2 t0 = 13.15 13.213.2513.313.3513.413.4513.513.55 13.613.6513.713.7513.813.8513.913.951414.05 14.114.15; l = 3.5331421843.546768029 3.561797643 3.578100715 3.595750783 3.614934283.635425983 3.657218272 3.680541115 3.705167836 3.731278025 3.758917911 3.788087888 3.818701015 3.850809619 3.884585223.919911828 3.956875992 3.995534794.035750835 4.077863059; b0 = 65,132, 5, 53; fun = inline(b(3)./tan(asin(sin(b(1)*pi/180)*sin(23.45*sin(b(4) - 80.25)*(1 - b(4)/9500)*pi/180)*pi/180)+cos(b(1)*pi/180)*cos(23.45*sin(b(4) - 80.25)*(1 - b(4)/9500)*pi/180)*pi/180)*cos(15*(t0 - 4*(116.3914-b(2)/60 -12)*pi/180),b,t0); b,r,j = nlinfit(t0, l, fun, b0) R = sum(r.2) 27 问题四程序： 视频信息提取 pro4_1 fileName = appendix4.avi; obj = VideoReader(fileName); numFrames = obj.NumberOfFrames;% 帧的总数 for k = 1 : 3000 : numFrames% 每隔 2 分钟读取视频数据 frame = read(obj,k); imwrite(frame,strcat(num2str(k),.jpg),jpg);% 保存帧 end 确定视频的拍摄地点 pro4_2 t0 = 8.98.938.9799.039.079.1 9.139.179.2 9.23 9.279.3 9.339.379.4 9.439.479.5 9.539.57; l = 1.2066373271.198149813 1.180361781 1.162158431 1.14860483 1.141405903 1.123626041 1.104154221 1.090613019 1.081729056 1.066922788 1.056365566 1.046205243 1.023782345 1.010257492 0.998404278 0.976014837 0.968395333 0.957844157 0.945992067 0.93247313; b0 = 25, 110, 3, 2; % b(4)这里以弧度表示 fun = inline(b(3)./tan(asin(sin(b(1)*pi/180)*sin(21.8292*pi/180)+cos(b(1)* pi/180)*cos(21.8292*pi/180)*cos(15*(t0 - 4*(116.3914-b(2)/60 -12)*pi/180).*sin(b(4) - asin(cos(21.8292*pi/180)*sin(15*(t0 - 4*(116.3914-b(2)/60 -12)*pi/180)./cos(asin(sin(b(1)*pi/180)*sin(21.8292*pi/180)+cos(b(1)* pi/180)*cos(21.8292*pi/180)*cos(15*(t0 - 4*(116.3914-b(2)/60 -12)*pi/180)./(2./tan(asin(sin(b(1)*pi/180)*sin(21.8292*pi/180)+c os