高教社杯全国大学生数学建模竞赛A题太阳影子定位

摘要通过太阳影子定位技术可以确定视频的拍摄地点和时间，为拍摄出更好的视频，掌握太阳影子的变化规律就变得尤为重要。本文主要综合运用了地理学、几何学、统计学、数学分析和高等代数等知识，并利用MATLAB,SPSS和mathematica等计算机软件，通过建立数学模型来研究影子长度的变化特征,进一步确定视频的拍摄地点和时间。针对问题一，首先我们通过分析影子长度的影响因素得到与影子长度的关系（见表达式六）整理计算之后，就得到了影子长度的数学模型。然后我们通过分析他们之间的关系，再利用MATLAB编程，得到了影子长度关于各个参数的变化规律（见图3到图7）。其次根据我们建立的模型，利用MATLAB编程画出了给定时间天安门广场3米高的直杆的太阳影子长度的变化曲线（见图8），然后在考虑折射率的情况下又画了一条变化曲线（见图9），最后进行了误差分析（见图10）。针对问题二，我们采用了测试分析法（数据分析法和计算机仿真相结合），通过分析各个参量之间的关系，先以影长为目标做回归，用模型一的模型，通过SPSS进行拟合得到多组数据，再用MATLAB进行检验得到符合的两组经纬度。然后我们又以太阳方位角为目标做回归,得到模型（见表达式12），其计算方法与影长做回归目标时一样。我们分步做了两次拟合，先用MATLAB拟合出经度，再做回归模型（见表达式14）最后得到经纬度和杆长。综上可知，肯定有一地点是在海南，还有一个地点可能在云南。针对问题三，我们用问题二中的多项式回归，得到回归模型（见表达式17和20）利用附件二得到的经纬度为和杆长，得到天数。利用附件三得到的经纬度为和杆长，得到天数针对问题四，首先运用MATLAB软件，根据画面灰度，运用MATLAB软件，把视频转化成二值图，求得影子端点的像素坐标，然后根据相似原理，把像素坐标转化成水平面上的坐标（消去了视角的影响），进而求得影子的长度。用以上方法求得的数据，运用多次拟合的方法，得到该地的经纬度为，日期未知时，得到的经纬度与其相似。【关键字】影子长度多项式拟合太阳方位角画面灰度问题重述如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用建立的模型画出已知时间天安门广场3米高的直杆的太阳影子长度的变化曲线。2.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。将你们的模型应用于附件1的影子顶点坐标数据，给出若干个可能的地点。3.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将你们的模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干个可能的地点与日期。4．先利用软件提取视频中的数据，再根据数据改善模型，求出若干个可能的拍摄点。当拍摄日期未知时，确定出拍摄地点与日期。问题分析由题可知，本题具体分析如下：问题一：本题要求建立影子长度变化的数学模型，这需要我们给出影子长度变化的影响因素。查阅文献了解到各个参量的定义及其表达式，然后联立即可得到影子长度变化的数学模型。分析影子长度关于各个参数的变化规律，首先我们要在保持其他参数不变的情况下，只改变一个参数，来研究影子长度的变化规律。对于具体问题的变化曲线，因为参数的值已经给出，带入模型，利用软件编程就可以画出它的变化曲线。问题二：本题要求根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点，然后将模型应用于附件1的影子顶点坐标数据，给出若干个可能的地点。先用附件一中给的数据即顶点的与坐标，计算出影子的长度,然后用SPSS做回归拟合，得出的数据再用MATLAB进行检验。问题三：问题三是问题二的拓展，建立数学模型确定直杆所处的地点和日期，比问题二多了一个未知量，我们可以采用问题二的模型和计算方法来解决本问题。问题四：问题四是前两问的具体应用，只要求出视频中影子的长度就可以运用前面的模型求解。对于求取视频中的影长，可以用MATLAB软件编写程序，设定恰当的灰度阀值，把视频转化为二值图像。从图片右下角开始扫描杆子和方块（杆子底座），求得影子右端、杆子底部的坐标。由于是在三维空间中拍摄的，图片中物体的长宽比与实际的长宽比不同，可以根据杆子底盘的长宽比求得物体实际的长宽比。根据相似度原理，由杆子的实际长度、图片中的像素维度等，求得像素与实际长度的比。最后，用影子右端、杆子底部的坐标、物体长宽比、像素与实际长度的比，求出影子的长度。得到上述数据之后，应用MATLAB进行多项式拟合，和应用SPSS软件进行非线性回归，两次拟合得到经纬度和日期。基本假设一年是365天；地球表面是球表面；地球的公转是正圆；大气层有折射无厚度；视频中杆的底盘是正方形，不考虑厚度；杆没有厚度。四、主要变量的符号说明为了便于描述问题，本文将问题中涉及的主要变量用下列符号来表示(如下表1所示)，有些变量将在文中用到时陆续说明。符号代表的含义符号代表的含义杆长影子的长度太阳时太阳赤纬角当地的纬度时角太阳的高度角折射后的太阳高度角日期序号折射率太阳方位角 原始像素的高度球坐标系切平面的坐标系表1五、模型的建立与求解5.1问题一建模和求解5.1.1.影子长度变化的数学模型（1）太阳时：太阳时是指以太阳日为标准来计算的时间，时间的计量以地球自转为依据，地球自转一周，计24太阳时，当太阳达到正南处为12:00。太阳时可以分为真太阳时和平太阳时，平太阳时假设地球绕太阳是标准的圆形，钟表所指的时称为平太阳时（简称为平时），我国采用东经120度经圈上的平太阳时作为全国的标准时间，即“北京时间”。（2）时角：时角是以正午12点为0度开始算，每一小时为15度，上午为负下午为正，即10点和14点分别为-30度和30度（因计算需要，把度数换算成了弧度）。因此，时角的计算公式为(弧度）(1)(2)其中为太阳时，是当地时间与北京时间的时间差（单位：小时）。赤纬角:赤纬角也称为太阳赤纬，即太阳直射纬度，其计算公式近似为(3)其中为日期序号，例如，1月1日为，10月22日为。太阳高度角：指某地太阳光线与该地作垂直于地心的地表切线的夹角。由图1可知，三点的坐标为由此可得：即太阳高度角的计算公式(4)其中为太阳高度角，为时角，为当时的太阳赤纬，为当地的纬度(北京的纬度为北纬)。图1由图2可知，即（5）图2（6）（6）将公式（1）（2）（3）带入（4），得出的值之后再代入（5）就能得到影子长度变化的数学模型（7）5.1.2影子长度关于各个参数的变化规律根据建立的数学模型，利用MATLAB编程（见附录），通过改变各个参数的值，观察影子长度的变化规律。首先，我们保持其他的变量不变，改变的值，由图3我们可以看到每天的最短影子长度在一年内呈周期性变化。图3（2）接下来我们保持其他的变量不变，改变经度的值，我们选取了东经和东经。由图4可以知道东经的度数越小，影子最短的时刻越往后。图4（3）其次我们在保持其他的变量不变，改变纬度的值，我们选取了北纬和北纬。由图5,图6可以知道北纬的度数越大，影子越长。图5图6（4）同样的，我们依然保持其他的变量不变，改变时间，由图7可以知道影子最短的时刻在中午附近，越远离正午，影子越长。图7（5）最后，我们依然保持其他的变量不变，改变杆长，由可知杆长与影子的长度成正比，即杆长越长，影子越长。5.1.3应用模型求解具体问题根据我们建立的模型，利用MATLAB编程画出2015年10月22日北京时间9:00-15:00之间天安门广场3米高的直杆的太阳影子长度的变化曲线如图8所示。图8然后我们又考虑了折射率，根据（8），得到的变化曲线如图9所示图9由图8和图9对比可以知道，是否考虑折射率差别小于千分之三，所以画了误差图10以便更好的观察误差。图10折射率还与大气层的厚度有关，但考虑到公式的复杂程度，这里不作考虑。5.2.问题二建模和求解5.2.1用一次拟合做回归5.2.1.1以长度为目标做回归根据问题一建立的数学模型（7）用附件一中给的数据即顶点的与坐标，计算出影子的长度,然后根据计算得到日期序号，然后利用SPSS以影长为因变量，不断给定杆长、纬度和一组初值做多次回归拟合(表达式见附录)，得到多组数据（部分数据在附录中给出），然后我们利用MATLAB进行检验。最后得到两组符合的拟合数据,经纬度表示就是经查阅地图可知这两个地点分别在海南岛上和云南省内。由于没有考虑太阳方位角的问题，分不清方向，所以我们将会在接下来做改进。5.2.1.2以太阳方位角为目标做回归太阳方位角是太阳在方位上的角度，它通常被定义为从北方沿着地平线顺时针量度的角。由图11可知太阳向量和法线的坐标为且（9）是正午的直射点，是杆的位置，是直射点，指向正东，指向正北，指向正上。在原坐标系上表示的是太阳的方向，切平面地表水平面（我们将其拿出来进行单独分析由于坐标系，过渡矩阵是在平面上的坐标。又因为图11由此可以得到然后结合各参量的表达式得到数学模型利用SPSS以为因变量不断给定杆长、纬度和一组初值做多次回归拟合,虽然我们考虑了太阳的方位角问题，但是画出的图形仍然存在误差。5.2.2用多次拟合做回归5.2.2.1确定该地的经度首先我们根据附录一中所给的数据用MATLAB画出散点图，由高度角、方位角和天文知识可以知道拟合的多项式为偶函数（关于正午时刻对称），为了得到更加精确的结果，我们采用了四次多项式模拟，然后用多项式拟合（程序见附录）多项式为（13），模拟出全天（从早上九点到下午五点）杆的影子长度图12然后利用mathematica求出曲线的根和顶点，得出最短影子长度出现的时间，进而确定该地的经度。5.2.2.2以影子长度为目标做回归根据上面求出的经度，就可以确定时角的值，然后利用SPSS以影子长度为因变量不断给定杆长和纬度初值做多次回归拟合，拟合的回归方程为（14）由以上公式得其中，可以把代入上式。由于我们版本的SPSS没有函数，所以我们将上面的回归方程化为（15）最后得到经纬度和杆长.综上可知，肯定有一地点是在海南，杆也可能在云南。5.3.问题三建模和求解5.3.1.1确定经度问题三是问题二的拓展，与第二问相比，该题的日期是未知量，即日期序号是未知量。所以这个问题我们也可用多次拟合得到结果。首先我们根据附录二中所给的数据用MATLAB画出散点图，我们依旧可以推得用偶次多项式进行拟合。然后用MATLAB求解得到经度。拟合的多项式为（16）模拟出全天（从早上九点到下午五点）杆的影子长度图13然后利用mathematica求出曲线的根和顶点，得出最短影子长度出现的时间,进而确定该地的经度为。5.3.1.2以影子长度为目标做回归根据上面求的经度，就可以确定时角的值，然后利用SPSS以影子长度为因变量不断给定参变量杆长和纬度初值做多次回归拟合，拟合的回归方程为（17)转化为SPSS用的方程为：（18）最后我们得到的经纬度为和杆长，得到天数即为5月20号，由对称性可知7月23号也是该题的答案，然后根据经纬度地图可以得到在新疆省内。5.3.2.1确定经度首先我们根据附录三中所给的数据用MATLAB画出散点图，用偶次多项式进行拟合。然后用MATLAB求解得到经度，拟合多项式为(19)模拟出全天（从早上九点到下午五点）杆的影子长度图14然后利用mathematica求出曲线的根和顶点，得出最短影子长度出现的时间进而确定该地的经度为5.3.2.2以影子长度为目标做回归根据上面求的经度，就可以确定时角的值，然后利用SPSS以影子长度为因变量不断给定参变量杆长和纬度初值做多次回归拟合，拟合的回归方程为(20)转化为SPSS用的方程为：(21)最后我们得到的经纬度为和杆长，得到天数即为11月3号，同样由对称性可知2月8号也是该问题答案。然后根据经纬度地图可以得到在湖北省内。5.4.问题四建模和求解5.4.1读取视频数据（1）读取图片运用MATLAB中的VideoReader函数可以知道视频共有40分钟41秒，每秒25帧，图片像素是1920×1080，影子变化不是很明显。为了分析结果明显，每隔2分钟取一次图片，以杆子和影子为图片中心，截取图片的像素为961×749，并设定恰当灰度阀值，将彩色图片转化为二值图像（程序见“将视频变为图片.m”），一共得到21张二值图像，其中的一张图片如下图：图15（2）求取影子端点像素坐标由图15可知，影子端点大致分布在图片的右下角，运用MATLAB软件编程（程序见“视频中影子端点的像素坐标.m”），从右下角开始扫描，寻找影子最右端，得到杆子的端点坐标如下图：图片序号1234567891011X’轴704704704706708708709710710712714Y’轴877869857849837827821809798790778图片序号12131415161718192021X’轴716714716717718718719719719719Y’轴773763755741734725717707690691表2由于杆子底座近似方形的，以杆子底部为其中心点，运用MATLAB软件编程（程序见“视频处理（杆子底部像素坐标）.m”），在图片中寻找方块，求其中心点，得出杆子底部的的像素坐标为。（3）求取坐标系下，影子端点的坐标以二值图片中杆子的底端为坐标原点，图片像素矩阵的纵列为轴，建立坐标系（单位：像素）,坐标图如下图：影子影子Y’轴X’轴O图16通过计算可以得出影子右端的坐标，如下表3：图片序号1234567891011轴-11-11-11-9-7-7-6-5-5-3-1轴808800788780768758752740729721709图片序号12131415161718192021轴1-112334444轴704694686672665656648638621622表3（4）求取XOY坐标系下，影子端点的坐标通过测量，得到原始图片测量高度为，原始图片中杆的测量高度为。根据相似原理，可得（22）其中原始图片中杆的像素高为，原始图片的像素高度。通过计算，可以得出像素。由题给的杆长，可得像素比，即图片中每337像素为1m长。以杆子的底端为坐标原点，正对镜头方向为X轴正方向，建立水平面XOY坐标系（单位：米），坐标图如下图：Y轴Y轴X轴O影子图17对于Y轴方向，物体在XOY坐标系中Y轴坐标、物体在坐标系中轴坐标和像素比的关系如下：由于图片来自三维空间，摄影机拍摄角度的不同，各物体在图片上的长度关系也会显示不同。通过测量，可以得出杆子底座长宽比=长/宽=4.29。对于X轴方向，物体在坐标系中X轴坐标、物体在坐标系中轴坐标的关系如下：（23）通过计算可以得出水平面XOY坐标系上影子右端的坐标，见表4：图片序号1234567轴-0.13989-0.13989-0.13989-0.11446-0.08902-0.08902-0.07630轴2.397632.373892.338282.314542.278932.249262.23145图片序号891011121314轴-0.06359-0.06359-0.03815-0.012720.01272-0.012720.01272轴2.195852.163202.139472.103862.089022.059352.03561图片序号15161718192021轴0.025430.038150.038150.050870.050870.050870.05087轴1.994071.973291.946591.922851.893181.842731.84570表4（5）求取影子长度通过影子右端的坐标可以计算出影子的长度，（24）影子长度（见附录--视频中影子长度）的散点图如下图：图185.4.2模型的建立与求解5.4.2.1经度的确立由以上方法得到时间，影子的长度数据，为了减少误差，我们根据散点图，排除两个异常点，并且为了数据的稳定性，我们对每两个相邻影子长度数据进行了取平均数处理。我们用上述处理过的数据用MATLAB画出散点图，用偶次多项式进行拟合。然后用MATLAB求解得到经度。拟合的多项式为（25）模拟出全天（从早上九点到下午五点）杆的影子长度图19然后利用mathematica求出曲线的根和顶点，得出最短影子长度出现的时间,进而确定该地的经度为5.4.2.2以影子长度为目标做回归根据上面求的经度，就可以确定时角的值，然后利用SPSS以影子长度为因变量不断给定纬度初值做多次回归拟合，拟合的回归方程为（26）由其中得到且的值代入上式。转化为SPSS用的回归方程为：（27）最后我们得到的经纬度为，然后根据经纬度地图可以得到在陕西省内。5.4.2.3当日期未知时的求解当日期未知时我们可以确定视频的拍摄地点和日期。我们依旧用多次拟合的方法来求解。首先用上述的数据和MATLAB拟合的方法得到经度，因为MATLAB拟合的方法不涉及到日期，所以经度不变，还是，时间也不变。根据上面得到经度，然后利用SPSS以影子长度为因变量不断给定纬度和日期序数初值做多次回归拟合，拟合的回归方程为（28）其中。转化为SPSS用的回归方程为：（29）最后我们得到的经纬度为，日期序号即6月22日。然后根据经纬度地图可以得到依旧在陕西省内。和第一次没有给日期的拟合结果相比较，地点基本相同，日期相差了22天。六、模型评价6.1模型的优点（1）在问题一的模型中，我们考虑了太阳的折射率，减小了影子长度的误差；（2）给出了在水平面上太阳的实际方向的变化函数；（3）给出了计算并描绘影子轨迹的MATLAB程序；（4）在问题二的求解中用了多种回归方法，互相验证；(5)在问题二和问题三当中，运用一次回归和二次回归，得到相似的结果。在问题四中，用MATLAB程序处理视频，使得误差较小。在问题四中，有无日期，可以得到相似的结果，说明了模型的有效性。6.2模型的不足：提取视频中的数据不够精确，忽略了底盘的高度6.3模型的改进:对于问题二的模型，我们没有考虑折射率的问题，接下来我们建立模型将折射率考虑在内。将影子长度分解到与轴上，得到（30）由表达式（n=1.000277）（7）可以得到（31）所以得到下面的数学模型（32）参考文献[1]伍光和，自然地理学，北京：高等教育出版社，2000。[2]吕林根许子道，解析几何，北京：高等教育出版社，2006.5。[3]茆诗松程依明濮晓龙，北京：高等教育出版社，2011.2。[4]华东师范大学数学系，数学分析上（下），北京：高等教育出版社，2010.7。[5]北京大学数学系几何与代数教研室前代数小组，高等代数，北京：高等教育出版社，2003.9。附录附录1：问题一全年最短影子变化曲线phi=.6965;pi=3.1415926;t0=12+0.2406;date=[295:365+294+1];L0=3;t=12;omega=(t-t0).*15/360*2*pi;theta=0.4092797.*sin(2*pi*(284+date)/365);w=asin(cos(phi).*cos(theta).*cos(omega)+sin(phi).*sin(theta));l=L0./tan(w);%%A=asin(cos(theta).*sin(omega)./cos(w));y1=cos(theta)*cos(omega)*sin(phi)-cos(phi)*sin(theta);x1=cos(theta)*sin(omega);xy=sqrt(x1.^2+y1.^2);x=x1./xy.*l;y=-y1./xy.*l;plot(date,l)xlabel('最短影子长度');ylabel('最短影子长度');title('全年最短影子变化曲线');gridon;附录2：问题一影子长度关于时间的变化曲线(没考虑折射率)phi=0.696512;pi=3.1415926;t0=12.240574;date=295;t=[9:0.1:15];omega=(t-t0)*15/360*2*pi;theta=0.4092797*sin(2*pi*(284+date)/365);w=asin(cos(phi)*cos(theta)*cos(omega)+sin(phi)*sin(theta));l=3./tan(w);plot(t,l)xlabel('时间');ylabel('影子长度');title('太阳影子长度的变化曲线');gridon;附录3：问题一考虑折射率后，影子长度关于时间的变化曲线

phi=0.696512;pi=3.1415926;t0=12.240574;date=295;n=1.000277

t=[9:0.1:15];

omega=(t-t0)*15/360*2*pi;

theta=0.4092797*sin(2*pi*(284+date)/365);

w=asin(cos(phi)*cos(theta)*cos(omega)+sin(phi)*sin(theta));

l=3./tan(w);

plot(t,l)

xlabel('时间');

ylabel('影子长度');

title('影子长度变化曲线(考虑了折射率)');

grid

on;附录4：问题一影子长度变化误差图

phi=0.696512;pi=3.1415926;t0=12.240574;date=295;n=1.000277

t=[9:0.1:15];

omega=(t-t0)*15/360*2*pi;

theta=0.4092797*sin(2*pi*(284+date)/365);

w=asin(cos(phi)*cos(theta)*cos(omega)+sin(phi)*sin(theta));

w1=pi/2-asin(sin(pi/2-w)/n)

l=3./tan(w);

l1=3./tan(w1);

plot(t,l)

xlabel('时间');

ylabel('未考虑折射率和考虑了折射率的影子差值');

title('误差图');

grid

on;附录5：问题二以杆长的影子做回归SPSS用的表达式l*SQRT(1/((COS(N)*0.98321997*COS(0.2618*(t-12-t1))+SIN(N)*0.18242391)**2)-1)附录6：问题二以杆长的影子做回归SPSS得到的部分数据赋予初值拟合的结果Nt1lNt1l0.3-30.56.356-10.8402.2500.3-20.51.415-7.2550.0450.3-10.5-0.0731.1602.2570.300.5-0.0731.1602.2570.310.50.3350.7812.0500.320.50.3350.7812.0500.4-20.50.3350.7812.0500.5-20.51.415-7.2550.0450.6-20.50.3350.7812.2501-230.3350.7812.0500.4-130.3350.7812.0500.5-13-0.0731.1602.2571030.4191.4392.7290.713-0.0731.1602.2571130.4191.4392.7290.6-33-2.8070.7812.0500.9-33-2.8070.7812.0501.3-110.3350.7812.0501.2110.3350.7812.0501.21.520.3350.7812.0501.2121.7274.7450.04500.10.50.4191.4392.7290.60.220.3350.7812.0500.7-11-0.0731.1602.2570.8230.4191.4392.7290.813-0.0731.1602.2571130.4191.4392.729012-0.0731.1602.2571.40.121.415-7.2550.450附录7：问题二以太阳方位角做回归SPSS用的表达式0.98321997*SIN(0.2618*(t-12-t1))/SQRT(1-(COS(N)*0.98321997*COS(0.2618*(t-12-t1))+SIN(N)*0.18242391)**2)附录8：问题二多次拟合做回归求经度，用MATLAB模拟出全天（从早上九点到下午五点）杆的影子长度%%用多项式拟合，y0=[1.149625826 1.182198976 1.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];x0=[14.7 14.75 14.8 14.85 14.9 14.95 15 15.05 15.1 15.15 15.2 15.25 15.3 15.35 15.4 15.45 15.5 15.55 15.6 15.65 15.7];plot(y0,x0,'*');x=9:0.01:17;p1=polyfit(x0,y0,4);%p1为影子长度y关于时间x的二项式函数y1=polyval(p1,x);plot(x,y1)附录9：问题二的多次拟合做回归用SPSS以影子长度为目标做回归的表达式l*SQRT(1/((COS(N)*0.98321997*COS(0.2618*(t–12.71))+SIN(N)*0.18242391)**2)-1)附录10：问题三附录二中所给的数据用MATLAB求经度y0=[1.247256205 1.22279459 1.198921486 1.175428964 1.152439573 1.12991747 1.10783548 1.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];x0=[12.68333333 12.73333333 12.78333333 12.83333333 12.88333333 12.93333333 12.98333333 13.03333333 13.08333333 13.13333333 13.18333333 13.23333333 13.28333333 13.33333333 13.38333333 13.43333333 13.48333333 13.53333333 13.58333333 13.63333333 13.68333333];plot(y0,x0,'*');x=9:0.001:17;p1=polyfit(x0,y0,4);%p1为影子长度y关于时间x的二项式函数y1=polyval(p1,x);plot(x,y1)附录11：问题三附录二用多次拟合做回归，用SPSS以影子长度为目标做回归的表达式l*SQRT(1/((COS(N)*COS(0.4028*SIN(2*3.1415*(284+c)/365))*COS(0.2618*(t-14.7))+SIN(N)*SIN(0.4028*SIN(2*3.1415*(284+c)/365)))**2)-1)附录12：问题三附录三中所给的数据用MATLAB求经度%%用多项式拟合y0=[3.533142184 3.546768029 3.561797643 3.578100715 3.595750783 3.61493428 3.635425983 3.657218272 3.680541115 3.705167836 3.731278025 3.758917911 3.788087888 3.818701015 3.850809619 3.88458522 3.919911828 3.956875992 3.99553479 4.035750835 4.077863059];x0=[13.15 13.2 13.25 13.3 13.35 13.4 13.45 13.5 13.55 13.6 13.65 13.7 13.75 13.8 13.85 13.9 13.95 14 14.05 14.1 14.15];plot(y0,x0,'*');x=9:0.001:17;p1=polyfit(x0,y0,4);%p1为影子长度y关于时间x的二项式函数y1=polyval(p1,x);plot(x,y1)附录13：问题三附录三用多次拟合做回归，用SPSS以影子长度为目标做回归的表达式l*SQRT(1/((COS(N)*COS(0.4028*SIN(2*3.1415*(284+c)/365))*COS(0.2618*(t-12.65))+SIN(N)*SIN(0.4028*SIN(2*3.1415*(284+c)/365)))**2)-1)附录14：问题四将视频变为图片的MATLAB程序

[filename,pathname,fileindex]=uigetfile('*.avi','选择视频文件','video.avi','Multiselect','on');

if

ischar(filename)

%只有选择了文件才进行以下计算

video=VideoReader([pathname

filename]);

%%General

Settings:

%Duration

时间

%Name

视频名称

%Path

路径

%Tag

=

VideoReader

%Type

=

mmreader

%UserData

=

[]

%%Video

Settings:

%BitsPerPixel=24

%FrameRate

视频采集速率

%Height

高度

%NumberOfFrames

总帧数

%VideoFormat

图像模式

%Width

=

宽度

LEN=video.NumberOfFrames;

%获得视频长度

dir=strcat(pathname,strrep(filename,'.avi',''),'\pic');

mkdir(dir);

fn=strrep(filename,'.avi','');

%thresh

=

graythresh(200);

for

k=1:3000:LEN-1

%若read到len，常会报错如下???

MATLAB:read:readTimedOut，read到len-1就好了

frame=im2bw(read(video,k),0.8);

imwrite(frame(166:914,800:1760),strcat(dir,'\',fn,'-avi-',int2str(k),'.bmp'),'bmp');%把每帧图像存入硬盘

end

elseif

iscell(filename)

navi=length(filename);

for

n=1:navi

video=VideoReader([pathname

filename{n}]);

%%General

Settings:

%Duration

时间

%Name

视频名称

%Path

路径

%Tag

=

%Type

=

mmreader

%UserData

=

[]

%%Video

Settings:

%BitsPerPixel=24

%FrameRate

视频采集速率

%Height

高度

%NumberOfFrames

总帧数

%VideoFormat

图像模式

%Width

=

宽度

LEN=video.NumberOfFrames;

%获得视频长度

dir=strcat(pathname,strrep(filename{n},'.avi',''),'\pic');

mkdir(dir);

fn=strrep(filename{n},'.avi','');

for

k=1:3000:LEN-1

%若read到len，常会报错如下???

MATLAB:read:readTimedOut，独到len-1就好了

frame=im2bw(read(video,k),0.8);

imwrite(frame(166:914,800:1760),strcat(dir,'\',fn,'-avi-',int2str(k),'.bmp'),'bmp');%把每帧图像存入硬盘

end

disp(strcat(num2str(n),'/',num2str(navi),'

:

"',filename{n},'"

Finished!',datestr(now,13)));

end

else

return

end

msgbox('所有帧提取完毕,已写入磁盘！','提示');附录15：问题四求视频中影子端点的像素坐标的MATLAB程序%%在运行完“视频变为图片.m”后，将图片的名称依次变为1,2,3...21.bmp，方便读取%%把图片的存储路径改为d:\Appendix4\pic%%读取视频转化后图片中影子端点的像素坐标%%需要多次手动改变图片的名称序1,2,3...21.bmpimax=749;jmax=961;i0=0;j0=0a=imread('d:\Appendix4\pic\21.bmp');%imshow(a)%通过显示图片，可知影子端点处于右下角p=0;q=0;forj=1:1500;%从图片右下角开始读fori=30:70;ifa(imax-i,jmax-j)+a(imax-i,jmax-j-1)++a(imax-i,jmax-j-2)==0%在二值图中，1表示“黑色”p=i;q=j;break;endendifp>0break;endendi0=imax-pj0=jmax-qfori=1:20forj=1:20a(i0+i,j0+j)=0;endendimshow(a)附录16：问题四用MATLAB做视频处理（杆子底部像素坐标）

%%在运行完“视频变为图片.m”后，将图片的名称依次变为1,2,3...21.bmp，方便读取

%%把图片的存储路径改为d:\Appendix4\pic

%%需要多次手动改变图片的名称序1,2,3...21.bmp

imax=749;jmax=961;isum=0;jsum=0;ii=0;jj=0;

bimax=100;bjmax=140;cc(1:100)=0;dd(1:140)=0;

a=imread('d:\Appendix4\pic\21.bmp');