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

摘要

通过太阳影子定位技术可以确定视频的拍摄地点和时间，为拍摄出更好的视频，掌握太阳影子的变化规律就变得尤为重要。本文主要综合运用了地理学、几何学、统计学、数学分析和高等代数等知识，并利用MATLAB,SPSS和mathematica等计算机软件，通过建立数学模型来研究影子长度的变化特征,进一步确定视频的拍摄地点和时间。

针对问题一，首先我们通过分析影子长度的影响因素得到与影子长度的关系（见表达式六）整理计算之后，就得到了影子长度的数学模型。

然后我们通过分析他们之间的关系，再利用MATLAB编程，得到了影子长度关于各个参数的变化规律（见图3到图7）。其次根据我们建立的模型，利用MATLAB编程画出了给定时间天安门广场3米高的直杆的太阳影子长度的变化曲线（见图8），然后在考虑折射率的情况下又画了一条变化曲线（见图9），最后进行了误差分析（见图10）。

针对问题二，我们采用了测试分析法（数据分析法和计算机仿真相结合），通过分析各个参量之间的关系，先以影长为目标做回归，用模型一的模型，通过SPSS进行拟合得到多组数据，再用MATLAB进行检验得到符合的两组经纬度。

然后我们又以太阳方位角为目标做回归,得到模型（见表达式12），其计算方法与影长做回归目标时一样。我们分步做了两次拟合，先用MATLAB拟合出经度，再做回归模型（见表达式14）最后得到经纬度和杆长。综上可知，肯定有一地点是在海南，还有一个地点可能在云南。

针对问题三，我们用问题二中的多项式回归，得到回归模型（见表达式17和20）

利用附件二得到的经纬度为和杆长，得到天数。利用附件三得到的经纬度为和杆长，得到天数

针对问题四，首先运用MATLAB软件，根据画面灰度，运用MATLAB软件，把视频转化成二值图，求得影子端点的像素坐标，然后根据相似原理，把像素坐标转化成水平面上的坐标（消去了视角的影响），进而求得影子的长度。用以上方法求得的数据，运用多次拟合的方法，得到该地的经纬度为，日期未知时，得到的经纬度与其相似。

【关键字】 影子长度 多项式拟合 太阳方位角 画面灰度

1、问题重述

如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。

1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用建立的模型画出已知时间天安门广场3米高的直杆的太阳影子长度的变化曲线。

2.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。将你们的模型应用于附件1的影子顶点坐标数据，给出若干个可能的地点。

3. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将你们的模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干个可能的地点与日期。

4．先利用软件提取视频中的数据，再根据数据改善模型，求出若干个可能的拍摄点。当拍摄日期未知时，确定出拍摄地点与日期。

2、问题分析

由题可知，本题具体分析如下：

问题一：本题要求建立影子长度变化的数学模型，这需要我们给出影子长度变化的影响因素。查阅文献了解到各个参量的定义及其表达式，然后联立即可得到影子长度变化的数学模型。分析影子长度关于各个参数的变化规律，首先我们要在保持其他参数不变的情况下，只改变一个参数，来研究影子长度的变化规律。对于具体问题的变化曲线，因为参数的值已经给出，带入模型，利用软件编程就可以画出它的变化曲线。

问题二：本题要求根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点，然后将模型应用于附件1的影子顶点坐标数据，给出若干个可能的地点。先用附件一中给的数据即顶点的与坐标，计算出影子的长度,然后用SPSS做回归拟合，得出的数据再用MATLAB进行检验。

问题三：问题三是问题二的拓展，建立数学模型确定直杆所处的地点和日期，比问题二多了一个未知量，我们可以采用问题二的模型和计算方法来解决本问题。

问题四：问题四是前两问的具体应用，只要求出视频中影子的长度就可以运用前面的模型求解。对于求取视频中的影长，可以用MATLAB软件编写程序，设定恰当的灰度阀值，把视频转化为二值图像。从图片右下角开始扫描杆子和方块（杆子底座），求得影子右端、杆子底部的坐标。由于是在三维空间中拍摄的，图片中物体的长宽比与实际的长宽比不同，可以根据杆子底盘的长宽比求得物体实际的长宽比。根据相似度原理，由杆子的实际长度、图片中的像素维度等，求得像素与实际长度的比。最后，用影子右端、杆子底部的坐标、物体长宽比、像素与实际长度的比，求出影子的长度。

得到上述数据之后，应用MATLAB进行多项式拟合，和应用SPSS软件进行非线性回归，两次拟合得到经纬度和日期。

3、基本假设

（1）一年是365天；

（2）地球表面是球表面；

（3）地球的公转是正圆；

（4）大气层有折射无厚度；

（5）视频中杆的底盘是正方形，不考虑厚度；

（6）杆没有厚度。

四、主要变量的符号说明

为了便于描述问题，本文将问题中涉及的主要变量用下列符号来表示(如下表1所示)，有些变量将在文中用到时陆续说明。

表1

五、模型的建立与求解

5.1问题一建模和求解

5.1.1.影子长度变化的数学模型

（1）太阳时：太阳时是指以太阳日为标准来计算的时间，时间的计量以地球自转为依据，地球自转一周，计24太阳时，当太阳达到正南处为12:00。太阳时可以分为真太阳时和平太阳时，平太阳时假设地球绕太阳是标准的圆形，钟表所指的时称为平太阳时（简称为平时），我国采用东经120度经圈上的平太阳时作为全国的标准时间，即“北京时间”。

（2）时角：时角是以正午12点为0度开始算，每一小时为15度，上午为负下午为正，即10点和14点分别为-30度和30度（因计算需要，把度数换算成了弧度）。因此，时角的计算公式为

(弧度） (1)

(2)

其中为太阳时，是当地时间与北京时间的时间差（单位：小时）。

（3）赤纬角:赤纬角也称为太阳赤纬，即太阳直射纬度，其计算公式近似为

(3)

其中为日期序号，例如，1月1日为，10月22日为。

（4）太阳高度角：指某地太阳光线与该地作垂直于地心的地表切线的夹角。由图1可知，三点的坐标为

由此可得：

即太阳高度角的计算公式

(4)

其中为太阳高度角，为时角，为当时的太阳赤纬，为当地的纬度(北京的纬度为北纬)。

图1

（5）由图2可知，即 （5）

图2

（6） （6）

将公式（1）（2）（3）带入（4），得出的值之后再代入（5）就能得到影子长度变化的数学模型

（7）

5.1.2影子长度关于各个参数的变化规律

根据建立的数学模型，利用MATLAB编程（见附录），通过改变各个参数的值，观察影子长度的变化规律。

（1） 首先，我们保持其他的变量不变，改变的值，由图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”），从右下角开始扫描，寻找影子最右端，得到杆子的端点坐标如下图：

表2

由于杆子底座近似方形的，以杆子底部为其中心点，运用MATLAB软件编程（程序见“视频处理（杆子底部像素坐标）.m”），在图片中寻找方块，求其中心点，得出杆子底部的的像素坐标为。

（3）求取坐标系下，影子端点的坐标

以二值图片中杆子的底端为坐标原点，图片像素矩阵的纵列为轴，建立坐标系（单位：像素）,坐标图如下图：

图16

通过计算可以得出影子右端的坐标，如下表3：

表3

（4）求取XOY坐标系下，影子端点的坐标

通过测量，得到原始图片测量高度为，原始图片中杆的测量高度为。根据相似原理，可得

（22）

其中原始图片中杆的像素高为，原始图片的像素高度。

通过计算，可以得出像素。由题给的杆长，可得像素比，即图片中每337像素为1m长。

以杆子的底端为坐标原点，正对镜头方向为X轴正方向，建立水平面XOY坐标系（单位：米），坐标图如下图：

图17

对于Y轴方向，物体在XOY坐标系中Y轴坐标、物体在坐标系中轴坐标和像素比的关系如下 ：

由于图片来自三维空间，摄影机拍摄角度的不同，各物体在图片上的长度关系也会显示不同。通过测量，可以得出杆子底座长宽比=长/宽=4.29。

对于X轴方向，物体在坐标系中X轴坐标、物体在坐标系中轴坐标的关系如下：

（23）

通过计算可以得出水平面XOY坐标系上影子右端的坐标，见表4：

表4

（5）求取影子长度

通过影子右端的坐标可以计算出影子的长度，

（24）

影子长度（见附录--视频中影子长度）的散点图如下图：

图18

5.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)在问题二和问题三当中，运用一次回归和二次回归，得到相似的结果。

（6）在问题四中，用MATLAB程序处理视频，使得误差较小。

（7）在问题四中，有无日期，可以得到相似的结果，说明了模型的有效性。

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('全年最短影子变化曲线');

grid on;

附录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('太阳影子长度的变化曲线');

grid on;

附录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得到的部分数据

附录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.bmp

imax=749;jmax=961;i0=0;j0=0

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

%imshow(a) %通过显示图片，可知影子端点处于右下角

p=0;q=0;

for j=1:1500; %从图片右下角开始读

for i=30:70;

if a(imax-i,jmax-j)+a(imax-i,jmax-j-1)++a(imax-i,jmax-j-2)==0 %在二值图中，1表示“黑色”

p=i;q=j;

break;

end

end

if p>0 break;

end

end

i0=imax-p

j0=jmax-q

for i=1:20

for j=1:20

a(i0+i,j0+j)=0;

end

end

imshow(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');

   for i=1:bimax;

       c(i)=0;

       for j=1:bjmax;

           b(i,j)=a(imax-bimax+i,j);

           c(i)=c(i)+b(i,j);

       end

       if c(i)

       end

   end

    for j=1:bjmax;

        d(j)=0;

       for i=1:bimax;

           d(j)=d(j)+b(i,j);

       end

       if d(j)

       end

   end

   for i=1:bimax

       if cc(i)==1 isum=isum+i;ii=ii+1;

       end

   end

   imid=floor(isum/ii)

    for j=1:bjmax

       if dd(j)==1 jsum=jsum+j;jj=jj+1;

       end

   end

   jmid=floor(jsum/jj)

       for i=1:5;

        for j=1:5;

    b(imid+i,jmid+j)=1;

        end

       end

    imshow(b)

    

    imax-bimax+imid

    jmid

附录17问题四用MATLAB程序求影子长度

%%x,y为影子端点像素坐标

%%a、b为杆子底部坐标

%%dist为影子的实际距离

%%xsbi为实际中1米在图片中的像素

%%shbi为由于拍摄视角造成三维空间中x、y的比例

x=[704 704 704 706 708 708 709 710 710 712 714 716 714 716 717 718 718 719 719 719 719];

y=[877 869 857 849 837 827 821 809 798 790 778 773 763 755 741 734 725 717 707 690 691];

a=715;b=69;

xsbi=337;

shbi=10.5/2.45;

xx=shbi*x/xsbi;xy=y/xsbi;

xa=shbi*a/xsbi;xb=b/xsbi;

cx=abs(xx-xa);cy=xy-xb;

for i=1:21

dist(i)=sqrt(cx(i).^2+cy(i).^2);

end

plot(dist,'*');

xlabel('图片序号');

ylabel('影子长度（单位：米）');

title('视频中影子长度')

x0=x-a;y0=y-b;%图片中坐标

x1=xx-xa;y1=xy-xb;%现实坐标

附录18：问题四求经度的MATLAB程序

%%用二项式拟合，

y0=[2.389854505 2.360232568 2.329913992 2.299019022 2.265844435 2.241888651 2.214762181 2.180452622 2.151972552 2.12185101 2.096477741 2.074222964 2.04751724 2.014937759 1.983945017 1.960311962 1.935241396 1.908689893

];

x0=[8.916666667 8.95 8.983333333 9.016666667 9.05 9.083333333 9.116666667 9.15 9.183333333 9.216666667 9.25 9.283333333 9.316666667 9.35 9.383333333 9.416666667 9.45 9.483333333

];

plot(y0,x0,'*');

x=9:0.01:17;

p1=polyfit(x0,y0,2);%p1为影子长度y关于时间x的二项式函数

y1=polyval(p1,x);

plot(x,y1)

本文来源：https://www.2haoxitong.net/k/doc/94d09789250c844769eae009581b6bd97e19bc41.html