范文一：太阳影子定位问题

太阳影子定位问题

摘要

如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，本文通过研究太阳影子定位技术，建立了基于太阳影子定位的数学模型。

第一问，首先，建立太阳高度角和影长的函数关系；接着，利用球面几何关系，推导出太阳高度角与太阳赤纬角、当地纬度和时角的函数关系；然后，推导赤纬角和时角的计算方法，参考有关资料用Bpirges赤纬角算法和Lamm时差算法对赤纬角和时角进行修正，提高模型精度；最后，将题目数据代入影长公式，借助matlab求得太阳影子长度的变化曲线，影长范围为?3.67m,7.64m?。

第二问，首先，根据影子长度和角度建立影子顶点坐标与当地经纬度的函数关系；接着，由于影长及坐标轴方向未知，引入相对杆长和方向角差分别消除杆长和坐标轴方向的影响；然后为了消除误差，使用最小二乘法拟合来进行计算；最后，用matlab编程求解当地经纬度，分别求得利用影长与角度拟合的两组结果，通过相互检验得到角度拟合出的结果更好，其解得直杆所处的地点坐标为?19.31N,110,18E?。

第三问，首先仍然根据影长与角度建立影子顶点坐标与当地经纬度、日期的函数关系，再利用最小二乘法进行曲线拟合，但是所得结果不能满足验证要求，其原因在于仅考虑影长或角度的关系而忽视了日期与经纬度之间的作用关系；因此，选择利用影子的顶点在地面上划过的弧长建立新的函数关系，该函数关系既体现影子长度又体现影子角度；最后，以该函数关系为基础用最小二乘法进行拟合，得到直杆所处的地点和日期，分别用影长和角度的函数关系对其解进行检验，结果可以相互验证，则得到附件二直杆所处的地点坐标为?38.99N,76.83E?，日期为2015年7月26日；附件三直杆所处的地点坐标为?29.07N,112.30E?，日期为2015年1月21日。 ，，

第四问，首先进行视频处理，每2分钟取视频一帧获得21张图像；然后，对图像进行裁剪、增加灰度对比度、锐化、边缘检测等处理，得到影子端点的图像坐标；接着，利用成像原理，根据图像坐标求得实际坐标；最后，利用前文方法拟合得到其地点和日期。直杆所处的地点坐标为?40.18N,122.91E?；当日期未知时，直杆所处的地点坐标为?39.86N,118.89E?，日期为2015年7月26日。 ，

关键词： 太阳高度角 太阳方位角 最小二乘法 相对弧长 图像识别

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

1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用你们建立的模型画出2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线。

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

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

4．附件4为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度为2米。请建立确定视频拍摄地点的数学模型，并应用你们的模型给出若干个可能的拍摄地点。

如果拍摄日期未知，你能否根据视频确定出拍摄地点与日期？

二、问题假设

(1) 地球为规则的球体

(2) 所有影长数据均在2015年采集

(3) 不考虑空气折射对影子长度与方向的影响

三、符号说明

问题一，本问是要求太阳影子长度的变化曲线，首先，考虑到将杆长与影子长度联系起来的是太阳高度角，所以如何求出准确的太阳高度角成为本文的关键；接着，尝试利用几何关系，推导太阳高度角，发现在当地纬度确定的情况下，只要求出太阳赤纬角和时角就可以；然后，我们考虑如何推导出太阳赤纬角与时角的求法，在查阅相关资料后，我们对这两个角度的求法进行了改进；最后，利用改进后的算法，考虑借助matlab求得太阳影子长度的变化曲线。

问题二，本问是要在建立以直杆底端为原点的坐标系，水平地面为xy平面的条件下，已知影子端点的坐标来求当地的经纬度。首先考虑根据影子长度建立影子顶点坐标与当地经纬度的函数关系，由于本问中杆长未知，考虑利用比商引入相对杆长消除杆长影响；接着，考虑根据影子方向建立影子顶点坐标与当地经纬度的函数关系，需要利用太阳方位角，由于本问中x、y轴方向与当地地理方向的关系并不已知，考虑引入x轴与当地正南方向的夹角为?，列出影子顶点坐标与太阳方位角函数关系，再利用做差引入影角变化消除?的影响；然后在求解两个函数关系时，为了更好地利用数据，消除误差，考虑利用最小二乘法拟合来进行计算；最后，考虑用matlab编程求解当地经纬度，并对结果进行分析。

问题三，本问与第二问类似，都是已知影子端点的坐标来求当地的经纬度，但区别是本问的拍摄日期也未知。首先，考虑仍然根据影长与影角建立影子顶点坐标与当地经纬度、日期的函数关系，再利用最小二乘法拟合求解；接着，考虑对影长与影角结果进行分析，发现解不能互相验证，考虑其原因在于单纯考虑影长或角度的关系时，日期与经纬度具有耦合作用；然后，考虑找一种既能体现影长又能体现方向的指标，通过这个指标来进行拟合，才能找到正确的位置和日期。考虑到影子的顶点在地面上划过的弧长既包含长度关系又包含角度关系，可以以弧长建立新的函数关系；最后，考虑以该函数关系为基础用最小二乘法进行拟合，来得到直杆所处的地点和日期，并且可以分别用影长和角度对其进行检验。

问题四，本问是要通过视频来求当地的经纬度与拍摄日期。首先，考虑进行视频处理，将视频录入matlab中，并将视频每2分钟取一帧，获得21张图像；然后，考虑进行图像处理，对图像进行裁剪、增加灰度对比度、锐化、边缘检测等处理，得到影子端点的图像坐标；接着，考虑利用成像原理，以图像坐标求得实际坐标；最后，考虑利用第三问中的弧长函数关系来得到直杆所处的地点和日期。

五、模型的建立与求解

5.1问题一模型的建立与求解

5.1.1引入概念

1.太阳直射点S

太阳直射点是地球和太阳中心连线与地球表面的交点。该点的太阳高度角为90°。太阳直射点随地球的自转和公转不断变化，经度坐标以天为单位周期变化，纬度坐标以年为单位周期变化。

2.赤纬角?

赤纬角又称太阳赤纬，它是地球赤道平面与太阳和地球中心的连线之间的夹角，即太阳直射点的纬度。赤纬角是由于地球公转造成的，以年为周期在南回归线与北回归线的范围内移动，即[?2326',2326']，促成了季节的变化。

3.时角?

时角是当地子午圈与过太阳直射点的子午圈之间的夹角，定义当地子午圈在东为正，在西为负，范围(?180°,180°]。时角与研究的地点和时间均相关。

5.1.2模型的建立 L

h

图1 太阳高度角地面示意图 图2 太阳高度角天体示意图 D

由图1易知l?Lcoth，求影长的关键在于求当时当地太阳高度角h，而h可根据图2进行求解。如图2，S为太阳直射点球，N为北极点，D为研究对象，则赤纬角、时角、研究对象纬度分别由?、?、?所示，太阳高度角即为球面弧DS的余角，球面三角形DSN中，边SN?90°??，DN=90°??，?SND=?，由球面三角形余弦公式可得：

cos(90°?h)=cos(90°??)cos(90°??)?sin(90°??)sin(90°??)sin(90°??)

化简可得：

sinh=sin?sin??cos?cos?cos? (1)

从（1）式中可看出，若当地纬度?已知，只需求得赤纬角?、时角?，即可求得当地太阳高度角h，进而求出影子长度l。

1.赤纬角?的计算

赤纬角是由于地球赤道平面与公转平面不重合产生的，它是一个近似以年为周期变化的量，其变化范围为?23?26?。赤纬角全年在南北回归线之间基本呈正弦曲线周期移动。每年夏至6月21日或22日赤纬达到最大值?23?26?，太阳位于地球北回归线正上空，是北半球日照时间最长的一天。随后赤纬角逐渐减少，至秋分9月21日或22日赤纬角等于0°时全球的昼夜时间均相等为。冬至12月21日或22日赤纬角减至最小值?23?26?，此时阳光斜射北半球，北半球昼短夜长。当春分3月21日或22日赤纬角又回到0°，如此周而复始形成四季，如下图：

[1]

图3 赤纬度变化图

图中箭头1表示上图观察方向，箭头2表示下图观察方向，S为过太阳中心平行于地球赤道平面的参考平面，PQ、MN分别是不同位置时地球中心距参考平面的距离，其中PQ表示夏至日时地球到参考平面的距离，其值PQ?r0sin?max，其中r0为日地距离，

?max为回归线纬度。 则任意时刻赤纬度??arcsinMN，MN?O'M?sin?max?r0?sin?公t?sin?max，其中r0

?公为地球的公转角速度，由于地球公转平均用时365.2天4，则

。综上，赤纬度计?公?36?0/36d5?.242?d0，.9t为地球离开春分点的时间（单位：天）

算公式为：

??arcsin(sin?公tsin?max) (2)

2.时角?的计算

时角是当地子午圈与过太阳直射点的子午圈之间的夹角，设其经度分别为?和?0，则有?????0。在一天即24小时内，太阳直射点绕地球转过一周，则其运动角速度等于地球自转角速度，?自?360?/24h?15?/h，又因为在北京时间12:00时，120°经线过

???太阳直射点，则在北京时间t时，?经线过太阳直射点，即120?（t-1）21??5

??0=12?0?（15t?）12。则时角的计算公式为：

?????120?+15（t?12） (3)

3.对计算公式的修正

事实上，由于地球公转轨道不是正圆，而是一个椭圆，日地距离r0及公转角速度?公都不是常数，地球日的时间也受公转的影响，赤纬角?及时角?的精确求解并不如此简单。历史上对赤纬角?及时角?的求解有多种方法，求解赤纬角?的方法包括Cooper法、

Whillier法、SpencerSpencer法、Stome法、Bpirges法，求解时角?的方法包括Lamm法、

法、Woolf法。

查阅相关文献，对以上方法进行简单介绍。

（1）赤纬角计算方法：

①Cooper法

?=23.45sin(2??284?n) 365

式中， n为该日期在一年中的日期序号，例如，当日期为1月1日时，n?1，当

日期为3月22日时，n?81。

②Spencer法

??0.006918?0.399912cos??0.070257sin??0.006758cos(2?)

?0.000907sin(2?)?0.002697cos(3?)?0.00148sin(3?)

式中，?为日角，??

③Stome法 2?(n?1)（n为该日期在一年中的日期序号）。 365

??sin?1?0.39795cos??

④Bpirges法 ??2?(n?173)??? ??365.242??

??0.3723?23.2567sin(wt)?0.1149sin(2wt)?0.1712sin(3wt)

?0.7580cos(wt)?0.3656cos(2wt)?0.0201cos(3wt) 式中，w?360，t?n?1?n0，365.2422

n0?78.801??0.2422(year?1969)??NT?0.25(year?1969)?

（2）时角计算方法：

当太阳中心连续两次通过地面静止观察者子午圈时称为一个太阳日。为了弥补由于地球绕日公转的速度不均匀以及赤道与黄道不在同一平面，导致一年中太阳日长短不一这一缺点，得到一个相对均匀的时间单位，人们引入一个假想的太阳，称为平太阳。平时生活所用的时间系统是基于平太阳时系统，所以在实时计算太阳时角时，关键是要将平时转换为视时。视时与平时的差称为时差eot ? equation of time?。所以，求解时角时应用时差修正当地时间在进行计算。与太阳赤纬相同，时差每时每刻都在变化，我们只能用一些不同精度的近似公式进行推算。

①Wloof法 [2]

eot?o.258cosx?7.416sinx?3.648cos(2x)?9.228sin(2x) 式中，x?

②Spencer法

360(n?1) 365.242

eot??0.000075?0.001868cos??0.032077sin??0.014615cos(2?)?0.04089(2?)?式中，??2?(n?1)

365

③Whillier法

eot?9.87sin(2B)?7.53cosB?1.5sinB(10) 式中，B?2?(n?81) 364

以上三种时差估算方法均以一年为周期，Lamm法选择时差的周期为四年，更贴近实际。

④Lamm法

2?kN2?kN??eot(N)??Akcos()?Bksin()?365.25365.25? k?0?5

式中，N为从每一个闰年开始为1至四年循环的最后一天，Ak和Bk的参数值见下表:

表1 Ak和Bk参数值

k

1

2 Ak/h Bk/h 2.0870e?4 9.2869e?3 ?5.2258e?2

?1.3077e?3

?2.1867e?3

1.5100e?4 0 ?1.2229e?1 ?1.5698e?3 ?5.1602e?3 ?2.9823e?3 ?2.3463e?4 3 4 5

以上是求解赤纬角和时角的典型算法，杜春旭等人在《一种高精度太阳位置算法》一文中对上述算法进行了分类比较，现引述其相关结论。 [2]

图4 四种赤纬角算法误差比较

上图显示近年每天世界时0时四种求解赤纬角算法的计算误差，可以看出，Bpirges法误差最小，最大计算误差为0.025。

图5 四种时差算法误差比较

上图显示近年每天世界时0时四种求解时角算法的计算误差，可以看出，Lamm法误差最小，最大计算误差为0.15。

综上所述，Bpirges赤纬角算法和Lamm时差算法精度远优于其他算法，本文将采用上述算法以提高模型的准确性。

5.1.3模型的求解

编写Bpirges赤纬角算法和Lamm时差算法的程序，通过matlab求解赤纬角?和时角?，进而求得影长。

下图为最终的到的2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线：

图6 太阳影子长度的变化曲线

从上图中可以看出，在北京时间9:00-15:00之间，位于天安门广场的3米高直杆在9:00时影长有最大值，为7.64m，随着时间推移，影长在9:00-12:18逐渐减小，到12:18影长达到最小值，为3.67m，随后，影长在12:18-15:00逐渐增大。

5.2问题二模型的建立与求解

5.2.1引入概念

1.太阳方位角F

太阳直射点和当地的大圆连线与当地南经线所成的球面角叫做太阳方位角，规定其

360?)，如图中F所示： 逆时针为正，范围[0，

图7 太阳方位角示意图

2.方位角的求解

同样在球面三角形DSN中运用球面三角形余弦公式可得：

cos(90°??)=cos(90°??)cos(90°?h)?sin(90°??)sin(90°?h)cos(180°?F)

化简可得：

cosF=sin?sinh?sin? (4) cos?cosh

5.2.2模型的建立

1.建立坐标系

根据题目要求，建立以直杆底端为原点的坐标系，水平地面为xy平面，但是xy轴的方向与地理方向的关系并不已知，直杆垂直于地面。

2.构建影子顶点坐标与当地经纬度的函数关系

（1）根据影长建立关系

由影子顶点坐标可计算得某时刻直杆影长：

li?根据太阳高度角亦可计算得某时刻直杆影长：

li?Lcothi

则

Lcothi?

在(5)式中xi、yi已知，hi可根据当地经纬度求得，而杆长L未知。所以引入相对杆长li=licothi消去杆长影响。最终建立的方程为：

?l1coth1

cothi?coth1

（2）根据影子方向建立关系 (6)

由于坐标系已建立，根据影子顶点坐标，还可以获得影子的方向信息，即影子顶点的坐标可与太阳方向角建立联系。又由于x、y轴方向与当地地理方向的关系并不已知，所以我们假设x轴与当地正南方向的夹角为?，逆时针为正，其范围(0,360]。

设某时刻影子与x轴夹角为?i，则

?i?arctan

又因为当地影子方向与太阳方向成180，则

Fi?(arctanyi xiyi???180)mod360 (7) xi

在（7）式中xi、yi已知，Fi可根据当地经纬度求得，而?未知。所以引入影角变化?Fi?Fi?F1??i??1消去杆长影响。最终建立的方程为：

Fi?F1?arctanyiy?arctan1 (8) xix1

可见我们可以根据某一时刻影子顶点的坐标，与当地经纬度建立两个方程，从而解得当地经纬度坐标。但是影子顶点坐标为测量值，各组数据均存在误差，而且只用一组数据求得的结果具有偶然性，没有充分利用所给数据。

3.数据拟合

为了更高效地利用数据，消除误差影响，提高计算精度，本文采用了最小二乘法对数据进行拟合。

下面简要介绍基本最小二乘法的统计学原理:

,,an)的依赖设物理量y与变量x1, x2,?, xl间有基于y?f( x1, x2,?,xl;a0,a1?

关系。

其中a0,a1,?,an是方程中需要确定的n?1个参数。最小二乘法就是通过m?m?n?1?个实验点(xi1,xi2,?,xil,yi) (i?1,2,?,m)确定出一组参数值(a0,a1,?,an),使由这组参数得出的函数值y?f(xi1,xi2,?,xi1;a0,a1,?,an)与实验值yi间的偏差平方和s(a0,a1,?,an)取得极小值。

最小二乘法一般采用多点测量来减小随机误差,使待求参数的个数小于方程式个数,即m?n?1。通过用最小二乘法进行统计处理,使未知数个数和方程个数相等,再进行求解得出a0,a1,?,an[3]。

可见，用最小二乘法拟合需要建立一个关系(x,y)?f(t,?,?)，本文之前建立的有关l、?F的关系均可以满足其要求。我们选择的做法是先利用其中一个关系进行拟合，并将拟合得到的结果代入另一个关系中检验，通过偏差平方和观察其吻合程度，选择吻合程度较高作为结果。

5.2.3模型的求解

使用matlab编程，分别用长度和方向关系根据最小二乘法拟合数据，即可求得对当地的经纬度，结果如下：

图8 拟合结果图

可见，用两种关系的拟合结果都不错，具有较高的吻合度。再将拟合结果代入另一个关系中进行检验，结果如下：

图9 检验结果图

可见，虽然利用影长关系所求解的拟合度更高，但其检验结果并不理想，而利用方向关系所求解在拟合和检验中都具有很高的吻合度。所以，该直杆所在地坐标为(19.31N,110.18E)。根据该坐标可以反求得到??103.4，即x轴方向为东偏北13.4

，

并不是正东方向。

下图为直杆所在地坐标在地图中的位置：

图10 第二问直杆地理位置

5.3问题三模型的建立与求解

在上一问中，本文先利用根据影长建立的关系进行拟合，并将拟合得到的结果代入根据影子方向建立的关系中检验，通过偏差平方和观察其吻合程度，选择吻合程度较高作为结果。利用影子长度与方向都可以有效地解决影子定位问题，两者相互验证，定位的准确度较高，因而，本问第一步先采用第二问的方法对问题进行求解。由于地球处于公转轨道的同一位置时，不同日期的影长和角度相同，故相同的情况会在各年重复出现，假设题目附件中的数据都是在2015年测的。

5.3.1使用第二问模型求解

在本问中当地的经纬度及日期未知，所以位置参数为经纬度与日期三个变量根据第二问的方法，依旧使用以下函数关系式：

（1）根据影长建立关系

（2）根据影子方向建立关系

?cothi cot1

Fi?F1?arctanyiy?arctan1 xix1

用matlab编程，以题目附件二中的数据为例，分别用影长和角度关系根据最小二乘法拟合数据，求得当地的经纬度及日期，结果如下：

图11 利用影长拟合结果

图12 利用角度拟合结果

由上图可以看出，直接使用第二问的方法所得的解不能互相验证，其原因在于单纯考虑影长或角度的关系时，日期与经纬度具有耦合作用，以影长因素为例，影长由太阳高度角决定，太阳高度角与当地纬度和赤纬度的差相关，当日期变化时，日期改变对影长的影响可以通过改变经纬度来抵消，如下图：

图13 日期与地点耦合关系图

这样就会产生一系列不同日期和地点的组合，并能使影长满足题中所给条件，而在这些日期与纬度的组合中，只有一点可以满足角度条件，但由于在拟合过程中单纯考虑影长条件，故无法找出这一正确组合。相应的单纯考虑角度因素，也会出现相应的问题，所以，解决此题的关键在于找到一种既能体现影长又能体现方向的指标，通过这个指标来进行拟合，才能找到正确的位置和日期。

通过以上结果可以看出，在增加一个未知参数后，单独利用影子长度与方向的方法的可靠性均降低了，这样的结果并不能令人信服。

5.3.2建立新的函数关系并求解

影子顶点与原点的连线反映影子的长度，连线与x轴的夹角反映影子的方向，而影子的顶点在地面上划过的弧长既包含长度关系又包含角度关系，下面根据弧长建立新的函数关系。 y

(xi,yi)?sixi?1

,yi?1

)

lili?1

ix

图14 弧长关系图

如图所示，由影子顶点坐标可以计算得第i段弧长：

?si?

再由余弦定理可以计算得第i

段弧长：

?si? 则

?xi?1?xi?

2?(yi?1?yi)2?li2?li?12?2cos(Fi?1?Fi)lili?1 (9) 对弧长进行归一化处理得到相对弧长，消去杆长的影响

?s?si?i?l1则，建立的新函数关系式为

?xi?1?xi?

2?cothi??cothi?1?cothicothi?1?(yi?1?yi)2????2cos(F?F) (10) ???i?1i2cothcothcoth?1??1?122

上式既能体现影长因素又能体现角度因素，以该式为基础用最小二乘法进行拟合，即可得到直杆所处的地点和日期。结果如下：

图15 附件二数据拟合结果

分别用影长和角度对其进行检验，得到如下结果：

图16 附件二数据检验结果

从图中可看出，该结果既满足影长条件，又满足角度条件，所以附件2中点的日期和坐标为：2015年7月26日，(38.99N，76.83E)。

5月1日，21同理，使用附件3中的数据可得到其日期和坐标为：201年

(29.0N7，112.E30，其检验结果见附录。下图为附件二和附件三中直杆在地图中的位置：

图17 附件二直杆地理位置 图18 附件三直杆地理位置

5.4问题四模型的建立与求解

本问不同于前三问，没有直接给出直杆影子顶点的坐标，而是提供了视频，视频中具有大量的信息，包括天气、植被、过往的行人车辆，以及每一帧中都可以提取的时间和影子顶点的信息，结合已做的工作，提取有关影子顶点和时间的信息，并用上文方法求解。

5.4.1视频处理

这里以matlab为工作平台，将题目所给视频导入到matlab中。下图为导入后得到的结果：

图19 视频文件信息

从图中可知，该视频文件时长40分40秒，帧率25Hz，即每分钟有1500帧图像。由于附件1-3中都给出21组数据，所以以两分钟（3000帧）为间隔对视频进行截取，得到21幅视频截图，以此为基础提取影子顶点和时间信息。

图20 视频截图

5.4.2图像处理

首先通过观察的方法在图上点选直杆底部像素，读取其图像坐标为(868,880)，观察法受人为影响较大，具有一定误差，但误差在可接受范围内。

然后，提取影子顶点坐标，有以下几个步骤：

1.图像裁剪

通过对图像的观察，发现影子顶点只在x??867,1700?,y??780,930?范围内移动，如图所示

故以此范围对图像进行裁剪以降低计算量和计算复杂度。

2.灰度处理

用数字图像的直方图代替图像的灰度密度函数。灰度直方图是一个离散函数，表示数字图像每一灰度级与该灰度级出现频率的对应关系。为了使图像清晰，将图像的灰度范围拉开，让灰度频率较小的灰度级变大，让灰度直方图在较大的动态范围内趋于一致。用图像f(x,y)的直方图代替灰度的分布密度函数pf(f)，直方图均化后的图像g为[4] [5]：

g?T[f]??pf(u)du 0f

直方图均化后的图像g(x,y)在该处的灰度sk为

sk?T[rk]??l?0knl N

上式中，N为像素总数，L为灰度级的个数，nk为频数，rk为f(x,y)在像素(x,y)处的灰度。

3.图像锐化处理

为凸显影子边缘，将图像转换为黑白图像：

4.边缘检测

在图像中，边界表明一个区域的结束，另一个区域的开始。边缘检测利用物体差异来实现，包括灰度、颜色或者纹理特征。边缘检测实际上就是检测图像特征发生变化的位置。边缘检测有两个内容：抽出边缘点，剔除间断点，并将边缘连接成完成的线：

5.求影子端点图像坐标

图像处理至此，影子端点已十分清晰，只需计算白色图像最右端点的坐标即可。 5.4.3坐标转化

上一步中所得的坐标为图像的像素坐标，想要求解待求量，还需将图像的像素坐标转化成实际坐标，而坐标的转化要涉及到坐标系的转化。这里主要有：

1.图像坐标系：在图像上以摄像机光轴与图像平面的交点为原点的坐标系。该坐标系下有两种坐标，一种是以像素为单位的图像坐标，另一种是以毫米为单位的图像坐标，分别以(u,v)和(x,y)来表示，若设每一个像素在x轴与y轴方向上的物理尺寸为?x,?y，则(u,v)和(x,y)有如下关系：

u?x

?u0?x

yv??v0

?y

2.摄像机坐标系：以摄像机光心为原点的坐标系，坐标以(X',Y',Z')来表示，摄像坐标系原点与图像坐标系原点间的距离为焦距f。

3.实际坐标系：以实际地物为原点的坐标系，坐标以(X,Y,Z)来表示。(X',Y',Z')和

(X,Y,Z)的关系是：

?X'??X?

?'??Y?Y?????? ?Z'??Z???????1??1??

其中?是4?4旋转矩阵。

这里，本文采用线性成像模型，示意图如下：

Z’

图21 线性成像示意图

从图中可以看出如下的比例关系：

fX'x?'

Z fY'

y?'

Z

将三种坐标系下的坐标结合起来，就得到了图像的像素坐标与实际坐标之间的关系：

?u?

???Z'?v????1??

?X?

?Y??? ?Z????1?

?表示一个3?4的旋转矩阵。

在本问中已知杆长，就相当于有了空间上的已知点，利用这些已知点就可以求出旋转矩阵，从而根据影子端点的像素坐标，得到它的实际坐标。

用matlab计算，得到影子端点的实际坐标如下表所示：

表2 影子端点的实际坐标

根据上表，应用前文的方法解得当地坐标(40.18N,122.91E)，根据该坐标可以反求得到杆长L?1.72，与题中数据相近。其拟合结果和地理位置如下：

图22 第四问拟合结果

图23 第四问直杆地理位置

六、模型评价

6.1模型的优点

（1）查阅资料，利用了Bpirges赤纬角算法和Lamm时差算法对求解赤纬角和时角的方法进行了改进，得到了较为精确的结果，提高了模型的准确性。

（2）利用影子的顶点在地面上划过的弧长建立新的函数关系，对问题进行求解，能够同时反映长度关系和角度关系，提高了模型的可靠性。

（3）在处理图像时，通过图像裁剪、灰度处理等工作使图像分析预处理数据时更加简便、清晰。 6.2模型的缺点

模型没有考虑大气层对光的折射，导致分析影长与角度时可能会出现误差。

七、参考文献

[1]互动百科，赤纬角， http://www.baike.com/wiki/%E8%B5%A4%E7%BA%AC%E8%A7%92，2015.09.12。

[2]杜春旭，王普，马重芳，吴玉庭，申少青，一种高精度太阳位置算法，新能源工艺，2008（2）：41-48,2008。

[3]宋殿瑞，宋文臣，刘鹏振，最小二乘法探讨，青岛化工学院学报，1998（3）：296-301,1998。

[4]张毓晋，图像处理和分析，北京：清华大学出版社，1999.3:43-49。 [5]阮秋琦，数字图像处理学，北京：电子工业出版社，2001.1:180-208。

八、附录

1.附件三拟合及检验结果：

图24 附件三数据拟合结果

图25 附件三数据检验结果

2.相关程序 利用影长拟合：

l1=sqrt(X1.^2+Y1.^2); ll1=l1/l1(1); pl0=[phi lamda];

[PL1,var11]=lsqcurvefit(@ratio_l,pl0,th1,ll1); theta1=atan(Y1./X1);

dtheta1=(theta1-theta1(1))*180/pi;

[PL12,var12]=lsqcurvefit(@d_F,pl0,th1,dtheta1); ll1_check=ratio_l(PL1,th1);

plot(th1,ll1);

hold on

plot(th1,ll1_check,'*r');

title('利用影长拟合结果','fontsize',20); xlabel('时间','fontsize',20); ylabel('相对影长','fontsize',20); lege=legend('原始数据','拟合数据'); set(lege,'fontsize',20);

text(14.8,1.5,['偏差平方和',10,num2str(var11)],'fontsize',20);

text(15.1,1.1,['(',num2str(PL12(1),4),',',num2str(PL12(2),5),')'],... 'fontsize',20); 利用角度拟合：

theta1=atan(Y1./X1);

dtheta1=(theta1-theta1(1))*180/pi; pl0=[phi lamda];

[PL12,var12]=lsqcurvefit(@d_F,pl0,th1,dtheta1); dth1_check=d_F(PL12,th1);

F1_check=calc_F(T,PL12(1),PL12(2),th1); Delta=F1_check-180-theta1*180/pi; Delta_bar=mean(Delta); plot(th1,dtheta1); hold on

plot(th1,dth1_check,'*r');

title('利用角度拟合结果','fontsize',20); xlabel('时间','fontsize',20); ylabel('影角变化','fontsize',20); lege=legend('原始数据','拟合数据'); set(lege,'fontsize',20);

text(15.3,-3,['偏差平方和',10,num2str(var12)],'fontsize',20);

text(14.75,-7,['(',num2str(PL1(1),4),',',num2str(PL1(2),5),')'],... 'fontsize',20);

范文二：太阳影子定位问题(数学建模 )

太阳影子定位的计算模型

摘要

本文研究的内容是如何提取视频中物体的太阳影子的变化数据，并且依据该数据利用太阳影子定位技术来确定视频拍摄的地点及拍摄日期。

针对问题一，首先我们明确了太阳影子的长度l和日序N、北京时间tb、纬度φx、经度φy、真实杆高?这几个因素有关，然后我们查阅了相关资料，通过太阳高度角α、时角β、赤纬δ这几个量将上述几个参数联系起来，从而建立了有关影长变化的数学模型，根据该模型利用MATLAB 做出影长关于各参量的变化图，并分析了其变化规律，最后将问题一中所给条件代入模型中，作出了太阳影子长度变化曲线。

问题二实为问题一中模型的逆向思考。问题实质为，在已知影长l、日序N和北京时间tb的的情况下求可能的经度φx、纬度φy。为此我们引入了一个新的参量影长比Ci以消除附件所给数据中杆高?未知的影响，再利用影长比Ci的变化建立了直杆地点的空间匹配模型，结合附件1所给数据，采用最小二乘法拟合的方法求解得到可能的目标地点为：海南省东方市（108.7118°E，19.2360°N）。

问题三相较问题二又多了一个代求参量日序N，我们仍沿用上问的模型，根据附件2、3的数据用最小二乘法拟合的方法求解，得到附件2可能的测量地点和日期为：新疆维吾尔自治区喀什地区，5月26日（79.8255°E，40.0325°N，N=146）；附件3可能的测量地点和日期为：湖北省神农架林区，11月17日（110.3718°E，31.5198°N，N=321）。 针对问题四，我们首先将视频材料以3分钟为间隔得到14张静态图片，对这些图片进行灰度处理，然后采用Otus 最大类间误差法把杆子和影子从背景中分割出来，从而得到影子实际长度和灰度值坐标的转换关系，最终得到了8:54-9:33每隔3分钟的影子长度，再利用问题二模型得出视频的拍摄地点为：内蒙古鄂尔多斯市（109.45°E, 39.65°N，），假设拍摄日期也未知，则利用问题三的模型求解出的拍摄日序为203即7月22日，与真实时间7月13日误差较小，说明模型较为精确，得到的结果较为可靠。

关键词：太阳影子 最小二乘拟合 图像处理 灰度值坐标 Otus最大类间误差法

1. 问题重述

如何确定一个视频所处的拍摄地和拍摄日期是视频的数据分析的重要内容，而其中一种方法就是太阳影子定位技术，即通过分析视频里物体太阳影子的变化来确定视频拍摄地和日期。在本文中，我们将研究以下几个问题：

1. 由于太阳影子的长度是随着各种因素的变化而变化的，因此建立影子长度变化的数学模型，然后分析长度关于各个参数的变化规律，将 2015 年 10 月 22 日北京时间 9:00-15:00 之间天安门广场（北纬 39度 54 分 26 秒, 东经 116 度 23 分 29 秒）3 米高的直杆这些条件带入模型，画出影子长度的变化曲线。

2. 根据某个固定直杆在水平地面上的太阳影子的顶点坐标数据来建立数学模型，用以确定其所处的地点。再将附件1的影子顶点坐标数据代入该模型，求解出所有可能的地点。

3. 根据某个固定直杆在水平地面上的太阳影子的顶点坐标数据来建立数学模型，用以确定其所处的地点和日期。将附件2 和附件 3 的影子顶点坐标数据带入模型，一次性给出若干个可能的地点与日期。

4. 根据附件4里视频中太阳影子的变化，并且已经通过一定方法得知了杆子的高度，请建立一个能确定视频拍摄地点的数学模型，然后用该模型去确定可能的拍摄地点，如果拍摄日期也未知，考虑能否根据视频确定出拍摄地点和日期。

2. 问题分析

针对问题一，首先我们明确了太阳影子的长度l和日期N、北京时间tb、纬度φx、经度φy、真实杆高?这几个因素有关，然后我们查阅了相关资料，通过太阳高度角α、时角β、赤纬δ这几个量将上述几个参数联系起来，从而可以建立有关影长变化的数学模型，将问题一中所给条件代入模型即可作出太阳影子长度变化曲线。

针对问题二，附件1给出了2015年4月18日时，某个固定直杆在水平地面上的太阳影子的顶点坐标数据，但规定的 x轴方向和y轴方向未知。可将顶点坐标数据转化为影长l，从而经问题二转化为问题一的逆向求解，由于未知参量较多，先通过引入影长比Ci消除杆高?未知的影响，再采用最小二乘法拟合的方法【2】求解即可。

针对问题三，附件2、3给出的是某个固定直杆在水平地面上的太阳影子的顶点坐

标数据，且规定的 x轴方向和y轴方向未知，相较第二问多了一个待求参量日序N。沿用上述模型，算法复杂度仅有较小的常数级数增长，沿用第二问的算法理论上仍可行。

针对问题四，我们首先将视频材料以3分钟为间隔得到14张静态图片，对这些图片进行灰度处理，然后采用Otus 最大类间误差法把杆子和影子从背景中分割出来，从而得到影子实际长度和灰度值坐标的转换关系，最终得到了8:54-9:33每隔3分钟的影子长度，再利用问题二和问题三的模型即可求得拍摄地点和日期。

3. 模型假设与符号系统

3.1模型的假设

（1）假设地球是一个规则的球体；

（2）假设直杆在测量过程中相对地面是完全静止的；

（3）假设太阳光线为平行光线，且测量时天气晴朗；

（4）忽略大气层对太阳光的折射、高山阻挡、海拔高度等因素的影响；

（5）假设一天中太阳直射点的纬度不变。

3.2符号系统

符号

l

h

α

φx

φy

δ

β

b

N

tb

t 符号说明 太阳影子的长度 直杆高度 太阳高度角 纬度 经度 赤纬 时角 地球公转的相对角度 以每年 1 月 1 日为第一天开始计算的日序 北京时间 当地地方时 单位 米 米 度 度 度 度 / / 天 / /

4. 模型的建立与求解

4.1问题一的分析与求解

4.1.1模型的准备

我们先将本问题中需要用到的名词进行解释，具体如下：

（1） 地方时t ：以一个地方太阳升到最高的地方的时间为正午 12 时，将连续两个正

午12 时之间等分为 24 个小时，所成的时间系统，称为地方时。

φy?120°

t=tb+其中：tb为北京时间；φy为当地的地理经度。

（2）太阳高度角 【1】：对于地球上的某个地点来讲，太阳光线的入射方向和该地与地心相连的地表切面的夹角。

（3）时角β【2】：这是天文学中的名词，某个天体的时角被定义为该天体的赤经（RA ）与当地的恒星时（LST ）的差值，即恒星时减去赤经。所以，一个天体的时角表示该天体是否通过了当地的子午圈。它的数值则表示了该天体与当地子午圈的角距离，并借用时间的单位，以小时来计量（1HA = 15度）。如果以北京时间为标准时间，则有：

π β=15（t?12）() °

（4）赤纬δ：是地球赤道平面与太阳和地球中心的连线之间的夹角，由于地球的公转，它每时每刻都在变化着，其中在春分和秋分时刻等于零，而在夏至和冬至时刻有极值，分别为正负 23.442°，太阳赤纬角在公转运动中任何时刻的具体值都是严格已知的，利用查得的公式【2】可以得出结果

δ= 0.006918?0.399912cos b+0.070257sin b?0.006758cos 2b

+0.000907sin 2b?0.002697cos 3b+0.00148sin 3b

b=2π(N?1)/365

其中：b是地球公转的相对角度，N是以每年 1 月 1 日为第一天开始计算的日序。 上述定义中所提及的各角度汇总见下图：

图 1 太阳高度角、太阳赤纬示意图

4.1.2影子长度变化的数学模型

根据假设，视太阳光为平行光，竖立在A 地的直杆在太阳光CD 的照射下形成影子AD ，示意图如下：

图 2 直杆投影的几何模型

由图可知，直杆的长度?与其太阳影子l的长度之间有如下关系：

l=?cot α

联系以上各关系式可以得到如下的有关影子长度变化的数学模型：

l=?cot α sin α=sin φxsin δ+cos φxcos δcos β δ= 0.006918?0.399912cos b+0.070257sin b?0.006758cos 2b+ 0.000907sin 2b?0.002697cos 3b+0.00148sin 3b

b=2π(N?1)/365 π° β=15（t?12）() 120°?φy t=tb?{

其中，自变量为N、tb、φy、φx、?，因变量为l。

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

根据上述模型，我们分别取杆高?、北京时间tb、经度φy、纬度φx、日序N为唯一变量，代入求解得到以下结论：

（1）当北京时间tb、经度φy、纬度φx、日序N固定时，影子长度l随着真实杆高?的增大而增大。

（2）当真实杆高?、经度φy、、纬度φx日序N固定时，影子长度l随着北京时间tb的增加呈现先减少后增加的趋势。如下图：

图 3 影长随北京时间的变化规律

（3）当真实杆高?、北京时间tb、纬度φx、日序N固定时，影子长度l随着经度φy的变化出现两个峰值。如下图：

图 4 影长随经度的变化规律

（4）当真实杆高?、北京时间tb、经度φy、日序N固定时，影子长度l随着纬度φx的变化呈现出从赤道处往两极方向增长的趋势。并且出现一个峰值。如下图：

图 5 影长随纬度的变化规律

（5）当真实杆高?、北京时间tb、纬度φx、经度φy固定时，影子长度l随着日序N的增加呈现先减少后增加的趋势。如下图：

图 6 影长随日序的变化规律

4.1.4按题目条件求解

根据上述所建立的模型，将2015 年 10 月 22 日北京时间 9:00-15:00 之间天安门广场（北纬 39度 54 分 26 秒, 东经 116 度 23 分 29 秒）3 米高的直杆这些条件代入，即已知

?=3

φx=39.9072254°N

φy=116.3913982°E

N=294

tb=from 9 to 15

可以得到如下太阳影子长度变化的曲线。

图 7 题目所给时间地点的影长变化曲线

从图中，我们发现在北京时间 9:00-15:00 的过程中，太阳影子的长度先从7.23米单调递减直至3.61米，此时为北京时间12：12，此后影子长度单调递增直至5.95米。

4.2问题二的分析与求解

4.2.1问题分析

问题二中，附件1给出的是测量日期为2015年4月18日时，某个固定直杆在水平地面上的太阳影子的顶点坐标数据，且规定的 x轴方向和y轴方向未知。此问实为第一问的变形，即第一问是在已知杆长?、日序N、北京时间tb、经度φx、纬度φy的情况下求影长l，而本问是在已知日序N、北京时间tb、影长l的情况下求可能的经度φx、纬度φy。

因此本题需要考虑如下问题：

1. 如何建立模型，解决所求目标变量变化的问题；

2. 如何消除杆长?未知的影响；

3. 未知变量较一问有所增加，需考虑如何解决所求目标变量增加的问题

由于附件中并没有给固定直杆的长度?，因此不能继续依照直接公式推导求定值的方法，需引入新的参量影长比Ci，并采用最小二乘法拟合的方法【3】解决上述3个问题。

4.2.2模型的建立

（1）确定目标函数

已知

sin αi=sin φxsin δ+cos φxcos δcos βi{，i=1,2, …, n li=?cot αi

其中，αi表示在第i个给定条件（日序N、时刻t、经度φx、纬度φy一定）时的太阳高度角, βi表示在第i个给定条件（日序N、时刻t、经度φx、纬度φy一定）时的时角，li表示在第i个给定条件（日序N、时刻t、经度φx、纬度φy一定）时的影长。

可得该地相邻时间段的影长比

??cot αicot αiCi=||=||=|| i+1i+1i+1

而附件一所给数据中，影长比 li

Di=

所以，目标函数为 i+1i+1

nminξ=∑(Ci?Di) 2

i=1

（2）建立直杆地点的空间匹配优化模型

nminξ=∑(Ci?Di) 2 Ci=|Di=cot αii+1i=1i+1i+1sin α=sin φxsin δ+cos φxcos δcos β δ= 0.006918?0.399912cos b+0.070257sin b?0.006758cos 2b+ 0.000907sin 2b?0.002697cos 3b+0.00148sin 3b 2π(N?1) b= π° β=15（t?12）() 120°?φy {t=tb?4.2.3

模型的求解

为求解上述模型，我们采用最小二乘法，代入附件一的条件，运用M 进行拟合，最终迭代，得到近似解区域如下：

图 8 附件1中数据多次迭代得到的近似解区域

得到目标函数最小的情况下的经度及纬度如下：

表 1 附件1计算所得地点

经度 108.7118°E

纬度 19.2360°N

地点 海南地区

我们在Google 地图上查询了上述地点的实际位置：

图 9 附件1的拍摄地点

我们计算了所得数据情形下的影长比与附件1所给影长比的方差之和，其大小均小于10?5，验证了我们得到结果的准确性。

4.3问题三的分析与求解 4.3.1问题分析

问题三中，附件给出的只有某个固定直杆在水平地面上的太阳影子的顶点坐标数据，相较第二问多了一个未知参数，即此问是在已知北京时间tb、影长l的情况下求可能的经度φx、纬度φy、日序N。因此解决此问时，算法的复杂度较前一问仅有较小的常数级数增长，沿用第二问的算法理论上仍可行。

4.3.2模型的建立

由问题分析已知，本问较前一问仅多了一个未知参量日序N，模型、算法仍可沿用，此处不再赘述。

针对本问，由于太阳高度角α公式最终含有日期参量：

sin α=sin φxsin δ+cos φxcos δcos β

δ= 0.006918?0.399912cos b+0.070257sin b?0.006758cos 2b+

0.000907sin 2b?0.002697cos 3b+0.00148sin 3b

2π(N?1) b={

所以，在进行最小二乘法拟合时有3个需赋初值的参数：经度φx、纬度φy、日序N。 目标函数仍为相邻时间段的影长比的最小二乘拟合

n

minξ=∑(Ci?Di) 2

i=1

4.3.3模型的求解结果对比分析

依据上述理论，通过MATLAB 程序，代入附件2、3中的数据，多次迭代，得到近似解区域如下：

图 10 附件2中目标函数小于10?5解的分布图

图 11 附件3中目标函数小于10?5解的分布图

并收敛得到以下结果：

表 2 附件2中数据所得结果

经度 79.8255°E

纬度 40.0325°N

日期（日序） 地点

5月26日（146） 吐鲁番地区

表 2 附件3中数据所得结果

经度 109.45°E

纬度 39.65°N

日期（日序） 11月17日（321）

地区 湖北省

我们在地图上查询了上述地点的实际位置分别为新疆维吾尔自治区喀什地区附近，湖北省神农架林区附近：

图 12 附件2的拍摄地点

图 13 附件3的拍摄地点

4.4问题四的分析与求解 4.4.1问题分析

问题四中，直接以视频的方式给出了固定直杆的影长变化规律。此时需先将图片形

式的影长变化规律以坐标的形式进行转换，再将其转换为现实的坐标形式。这样就可将

问题转化为“已知日序N、北京时间tb、影长l，求可能的经度φx、纬度φy”，从而可利用问题二的模型，求出拍摄地点。

4.4.2数据预处理

影长需要以实际坐标中的形式进行展示，所以需要对视频资料进行预处理操作。具体步骤如下：

（1） 获取静态图片

视频展示的是一个动态过程，而确定视频中影子端点的坐标需要在静态画面下进行。为了与附件的时间间隔一致，我们选取了3分钟为时间间隔截取了视频中共14张图片, 并进行灰度处理，某一帧图片的灰度直方图如下：

图 14 某一帧图片的灰度直方图

（2） 目标与背景的分离

我们采用Otus 最大类间误差法把杆子和影子从背景中分割出来 ，其思想是：根据灰度特性，将图像分为目标和背景2部分，如果目标和背景之间的类间差越大，则说明构成图像的2部分的差别越大，因此类间方差最大的分割就意味着错分概率最小，计算以每个灰度值为阈值的分割的类间方差，那么类间方差最大的值即为阈值, 从而把杆子和影子从背景中分割出来。

（3） 图片的处理

每张图片的分辨率为1080×1920，我们利用MATLAB 读取了这些图片的灰度值，得到了灰度值矩阵，然后以影长在图片中灰度值的变化来刻画其真实长度的变化。具体如下:

a. 参考原点的选取

我们以直杆的接触地面的点为原点，提取其灰度值为（x 0, y 0） b. 参考坐标轴的选取

我们以参考原点灰度值的横坐标为x轴，纵坐标为y 轴 c. 实际影长的换算

由于题目中已经告诉我们实际的直杆长度，我们在灰度图中提取出直杆端点的灰度值（x1，y 1），影子端点的灰度值为（x2，y 2），则实际长度与灰度值坐标下的转换如下：

√（x1?x0）+（y1?y0）

√（x2?x0）+（y2?y0）

其中，?为真实杆长2米，l为真实影长。 最终我们从视频中提取得到如下数据：

北京时间 8:54 8:57 9:00 9:03 9:06 9:09 9:12 9:15 9:18 9:21 9:24 9:27 9:30 9:33

影子长度 2.144 2.099 2.056 2.011 1.969 1.926 1.886 1.853 1.800 1.763 1.726 1.685 1.645 1.621

4.2.3问题的求解

得到了影长的数据后，采用问题二和问题三的求解方法，此处不再赘述，最终求得结果如下：

经度 109.45°E

纬度 39.65°N

地区 内蒙古境内

我们在地图上查询了上述地点的实际位置为内蒙古鄂尔多斯市附近：

图 15 附件4拍摄地点

假设拍摄时间也未知，利用问题三的模型，可以求得拍摄日序为203，即7月22日，与真实的拍摄日期7月13日天的第193日较为接近，说明模型较为精确，得到的结果较为可靠。

5. 模型的优缺点及改进

5.1模型优缺点

优点：

1. 对于第二问，通过引入影长比Ci消除杆高?未知的影响，大大减少了时间复杂度； 2. 模型简单实用，易于计算机操作；

3. 模型计算结果准确性较大，定位较精确，可靠度较高。 缺点：

1. 对于第三问，由于引入了新参量日期，采用最小二乘法的时间复杂度较大； 2. 忽略了大气层对太阳光的折射、高山阻挡、海拔高度等因素的影响，使模型在不同情况下的适用性有所局限；

3. 本文所建立的模型对拍摄日期的求解无法精确到具体年份；

4. 直接用灰度值的变化的比例来刻画影子长度的变化规律，存在一定误差。

5.2模型的改进

1. 可以考虑大气层对太阳光的折射，从而增强模型的适用性。

2. 由于拍摄存在一定角度，可以再详细考虑图片中参考坐标系与世界坐标的转化。

6. 参考文献

[1] 费云霞 王春顺, 对太阳高度角的了解及其计算方法[J],中小企业管理与科技（上半月）,2008(01).

[2] 刘砚刚, 四季与赤纬角和时角[J],太阳能,1991(03). [3] 非线性最小二乘法

.

[4] 崔灿，张国华，一种基于单幅图像双消失点的摄像机标定方法，科学与技术工程，12 期.

[5]太阳高度角, 百度百科,

，9月1日.

[6]太阳赤纬, 百度百科,

，9月1日.

7. 附录

程序

%载入数据并保存

clear;clc; close all ;

%data=xlsread('附件1-3.xls',' 附件1'); %data=xlsread('附件1-3.xls',' 附件2'); data=xlsread(' 附件1-3.xls' , ' 附件3' ); save data data

%模型的参数影响

clear;clc; close all %%

%h 高度 h1=3

%N 日期 N1=4.18 31+28+31+18=108 %t_b 北京时间 t_b

%t_c 当地实际时间于北京时间的时差 %phi_x 为当地的地理经度 %%

%phi_y 为当地的地理经度 %l-h

phi_y=116;%经度 phi_x=39;%纬度 N=294;%日期

t_b=9;%北京时间

h=0:0.1:10;%杆高 figure;

l=shadow(phi_y,phi_x,N,t_b,h); %subplot(3,2,1) plot(h,l)

title(' 影子长度随着杆高的变化' ) xlabel(' 真实杆高h(m)') ylabel(' 影长/m') %% %l-t_b

phi_y=116;%经度 phi_x=39;%纬度 N=294;%日期 h=3;%杆高

t_b=7.5:0.1:17.5;%北京时间 l=shadow(phi_y,phi_x,N,t_b,h);

figure;

%df=diff(l); %find(df)

%plot([7.6:0.1:17.5],df) %subplot(3,2,2) plot(t_b,l)

title(' 影子长度随着北京时间的变化' ) xlabel(' 北京时间t_b') ylabel(' 影长/m') xlim([7.5 17.5]) %% %l-phi_y

phi_y=-180:180;%经度 phi_x=39;%纬度 N=294;%日期 h=3;%杆高

t_b=10;%北京时间

l=shadow(phi_y,phi_x,N,t_b,h); %subplot(3,2,3) plot(phi_y,l)

title(' 影子长度随着经度的变化' ) xlabel(' 经度\phi_y') ylabel(' 影长/m') %% %l-phi_x

phi_y=116;%经度 N=294;%日期 h=3;%杆高

phi_x=-90:90;%纬度 t_b=10;%北京时间

l=shadow(phi_y,phi_x,N,t_b,h); %subplot(3,2,4) figure; plot(phi_x,l)

title(' 影子长度随着纬度的变化' ) xlabel(' 纬度\phi_x') ylabel(' 影长/m') %% %l-N

phi_y=116;%经度 phi_x=39;%纬度 N=1:365;%日期 h=3;%杆高

t_b=10;%北京时间

l=shadow(phi_y,phi_x,N,t_b,h); %subplot(3,2,5) figure; plot(N,l)

title(' 影子长度随着日期的变化' ) xlabel(' 日期N' )

ylabel(' 影长/m')

function [l,alpha] = shadow (phi_y,phi_x,N,t_b,h) t_c=(120-phi_y)/15;%时差 st=t_b-t_c;%真太阳时 t=15*(st-12);%度数

b=2*pi*(N-1)/365;%弧度√

delta=(0.006918-0.399912*cos(b)+ 0.070257*sin(b)-... 0.00148*sin(3*b) );%弧度√

0.006758*cos(2*b)+0.000907*sin(2*b)-0.002697*cos(3*b)+ ...

alpha=asin(sin(phi_x*pi/180)*sin(delta)+cos(phi_x*pi/180)*cos(delta)*cos(t*pi/180)); l=h*abs(cot(alpha));%h为自变量 end

%问题一的程序 clear;clc; close all %p1

phi_y=116.3913;%经度 phi_x=39.9072;%纬度 N=294;

t_b=9:0.1:15;%北京时间 h=3;%杆高 figure;

l=shadow(phi_y,phi_x,N,t_b,h); plot(t_b,l)

title(' 影子长度随着北京时间的变化' ) xlabel(' 北京时间t_b') ylabel(' 影长/m') figure;

for N=1:40:360%日期 hold on end

plot([9.1:0.1:15],diff(shadow(phi_y,phi_x,N,t_b,h)))

lgd = num2str((1:40:360)','N=%d');

legend(lgd,'location' , 'Best' )

title(' 不同日期的影长变换率' )

%%问题二最优便利搜索

%%

%p2

%便利搜素

clear;clc;

close all

load data

D=sqrt(data(1:20,2).^2+data(1:20,3).^2)./sqrt(data(2:21,2).^2+data(2:21,3).^2); N=108;

xi=1e8;

tic

t_b=14.7:0.05:15.7;%北京时间

for phi_x=-90:1:90

for phi_y=-180:1:180

alpha=threat(phi_y,phi_x,N,t_b);

C=(cot(alpha(:,1:20))./cot(alpha(:,2:21)))';

if sum((C-D).^2)<>

phi_x_min=phi_x;

N_min=N;

end

end

end

for phi_x=phi_x_min-1:0.01:phi_x_min+1

alpha=threat(phi_y,phi_x,N,t_b); xi=sum((C-D).^2); phi_y_min=phi_y; for phi_y=phi_y_min-1:0.01:phi_y_min+1

C=(cot(alpha(:,1:20))./cot(alpha(:,2:21)))';

if sum((C-D).^2)<>

lat=phi_x;

riqi=N;

end

end xi=sum((C-D).^2); Log=phi_y;

end

toc

%%问题二的非线性拟合

clear;clc

close all

phi_x0=[];

phi_y0=[];

error=[];

for phi_x=-90:1:90 % 14-28

x=[phi_x phi_y];

if c<1e-5 for="" phi_y="-180:1:180" %="" 105="" -145="" c="">

error=[error c];

phi_x0=[phi_x0 phi_x];

end

end

end phi_y0=[phi_y0 phi_y];

scatter3(phi_x0,phi_y0,error,'filled' , 'MarkerEdgeColor' , 'k' , 'MarkerFaceColor' ,[1 0.8 0.5]); %scatter3(phi_x0,phi_y0,error,'filled','MarkerFaceColor',[1 0.8 0.5]);

xlabel('\phi_x')

zlabel('\xi')

%% ylabel('\phi_y')

%非线性最小二乘拟合

c0_x=phi_x0(k); k=find(error==min(error));

c0_y=phi_y0(k);

c0=[c0_x c0_y]; %c0=[50 100];%x初始点

%opt=optimset('tolx',1e-10,'tolfun',1e-10);

opt=optimset('Algorithm' , 'Levenberg-Marquardt' , ...

'Display' , 'iter' , 'tolx' ,1e-18, 'tolfun' ,1e-18, 'MaxFunEvals' ,4000, 'MaxIter' ,4000);

[ c, resnorm ,residual]=lsqnonlin('shadow_fit2',c0,[],[],opt);

function y= shadow_fit2(x,t_b)

%phi_x x(1)

load data

N=108; %phi_y x(2) D=sqrt(data(1:20,2).^2+data(1:20,3).^2)./sqrt(data(2:21,2).^2+data(2:21,3).^2); t_b=14.7:0.05:15.7;%北京时间

t_c=(120-x(2))/15;%时差

st=t_b-t_c;%真太阳时

t=15*(st-12);%度数

b=2*pi*(N-1)/365;%弧度√

delta=(0.006918-0.399912*cos(b)+ 0.070257*sin(b)-...

0.00148*sin(3*b) );%弧度√ 0.006758*cos(2*b)+0.000907*sin(2*b)-0.002697*cos(3*b)+ ...

alpha=asin(sin(x(1)*pi/180)*sin(delta)+cos(x(1)*pi/180)*cos(delta)*cos(t*pi/180)); C=cot(alpha(:,1:20))./cot(alpha(:,2:21));

C=C';

end y=sum((C-D).^2);

%%问题三的第一问

%最优搜索

%%

%p2

%便利搜素

clear;clc;

close all

load data

D=sqrt(data(1:20,2).^2+data(1:20,3).^2)./sqrt(data(2:21,2).^2+data(2:21,3).^2); %N=108;

xi=1e8;

tic t_b=14.7:0.05:15.7;%北京时间

for N=1:365

for phi_x=-90:1:90

for phi_y=-180:1:180

alpha=threat(phi_y,phi_x,N,t_b);

C=(cot(alpha(:,1:20))./cot(alpha(:,2:21)))';

if sum((C-D).^2)<>

phi_x_min=phi_x; xi=sum((C-D).^2);

phi_y_min=phi_y;

N_min=N;

end

end

end

end

for N=1:365

for phi_x=phi_x_min-1:0.01:phi_x_min+1

alpha=threat(phi_y,phi_x,N,t_b); for phi_y=phi_y_min-1:0.01:phi_y_min+1

C=(cot(alpha(:,1:20))./cot(alpha(:,2:21)))';

if sum((C-D).^2)<>

lat=phi_x;

riqi=N;

end

end

end

end

toc xi=sum((C-D).^2); Log=phi_y;

%最小二乘拟合

clear;clc

close all

phi_x0=[];

error=[];

riqi=[];

tic phi_y0=[];

for N=1:1:365

N

for phi_x=10:1:40 %

for phi_y=70:1:80 % 70-140

x=[phi_x phi_y N];

if c<1e-6 c="">

error=[error c];

phi_x0=[phi_x0 phi_x];

riqi=[riqi N];

end phi_y0=[phi_y0 phi_y];

end

end

end

scatter3(phi_x0,phi_y0,riqi,'filled' , 'MarkerEdgeColor' , 'k' , 'MarkerFaceColor' ,[1 0.8 0.5]); %scatter3(phi_x0,phi_y0,error,'filled','MarkerFaceColor',[1 0.8 0.5]);

xlabel('\phi_x')

zlabel('N' )

%% ylabel('\phi_y')

%非线性最小二乘拟合

c0_x=phi_x0(k);

riqi0=riqi(k); k=find(error==min(error)); c0_y=phi_y0(k);

%c0=[50 100];%x初始点

c0=[c0_x c0_y riqi0];

%c0=[-30 80 12];

%c0=[30 110 150];

%opt=optimset('tolx',1e-10,'tolfun',1e-10);

opt=optimset('Algorithm' , 'Levenberg-Marquardt' , ...

'Display' , 'iter' , 'tolx' ,1e-18, 'tolfun' ,1e-18, 'MaxFunEvals' ,1e4, 'MaxIter' ,1e4);

[ c, resnorm ,residual]=lsqnonlin('shadow_fit31',c0,[],[],opt);

c

load data

l=sqrt(data(:,2).^2+data(:,3).^2);

t_b=12.68:0.05:13.68;%北京时间

h = ganchang (c(1),c(2),c(3),t_b,l);

function y= shadow_fit31(x)

%phi_x x(1)

%phi_y x(2)

load data %N x(3)

D=sqrt(data(1:20,2).^2+data(1:20,3).^2)./sqrt(data(2:21,2).^2+data(2:21,3).^2); t_b=12.68:0.05:13.68;%北京时间

t_c=(120-x(2))/15;%时差

st=t_b-t_c;%真太阳时

t=15*(st-12);%度数

b=2*pi*(x(3)-1)/365;%弧度√

delta=(0.006918-0.399912*cos(b)+ 0.070257*sin(b)-...

0.00148*sin(3*b) );%弧度√ 0.006758*cos(2*b)+0.000907*sin(2*b)-0.002697*cos(3*b)+ ...

alpha=asin(sin(x(1)*pi/180)*sin(delta)+cos(x(1)*pi/180)*cos(delta)*cos(t*pi/180));

C=cot(alpha(:,1:20))./cot(alpha(:,2:21));

C=C';

end y=sum((C-D).^2);

%%问题三的第二问

clear;clc

close all

phi_x0=[];

phi_y0=[];

error=[];

riqi=[];

tic

for N=1:1:365

N

for phi_x=10:1:30 %

for phi_y=100:1:110 % 70-140

x=[phi_x phi_y N];

if c<1e-5 c="">

error=[error c];

phi_x0=[phi_x0 phi_x];

riqi=[riqi N];

end

end

end

end

phi_y0=[phi_y0 phi_y];

scatter3(phi_x0,phi_y0,riqi,'filled' , 'MarkerEdgeColor' , 'k' , 'MarkerFaceColor' ,[1 0.8 0.5]); %scatter3(phi_x0,phi_y0,error,'filled','MarkerFaceColor',[1 0.8 0.5]);

xlabel('\phi_x')

zlabel('N' )

%% ylabel('\phi_y')

%非线性最小二乘拟合

c0_x=phi_x0(k);

riqi0=riqi(k); k=find(error==min(error)); c0_y=phi_y0(k);

%c0=[50 100];%x初始点

c0=[c0_x c0_y riqi0];

%c0=[-30 80 12];

%c0=[30 110 150];

%opt=optimset('tolx',1e-10,'tolfun',1e-10);

opt=optimset('Algorithm' , 'Levenberg-Marquardt' , ...

'Display' , 'iter' , 'tolx' ,1e-18, 'tolfun' ,1e-18, 'MaxFunEvals' ,1e4, 'MaxIter' ,1e4);

[ c, resnorm ,residual]=lsqnonlin('shadow_fit32',c0,[],[],opt);

c

load data

l=sqrt(data(:,2).^2+data(:,3).^2);

t_b=13.15:0.05:14.15;%北京时间

h = ganchang (c(1),c(2),c(3),t_b,l);

function y= shadow_fit32(x)

%phi_x x(1)

%phi_y x(2)

load data %N x(3)

D=sqrt(data(1:20,2).^2+data(1:20,3).^2)./sqrt(data(2:21,2).^2+data(2:21,3).^2); t_b=13.15:0.05:14.15;%北京时间

t_c=(120-x(2))/15;%时差

st=t_b-t_c;%真太阳时

t=15*(st-12);%度数

b=2*pi*(x(3)-1)/365;%弧度√

delta=(0.006918-0.399912*cos(b)+ 0.070257*sin(b)-...

0.00148*sin(3*b) );%弧度√ 0.006758*cos(2*b)+0.000907*sin(2*b)-0.002697*cos(3*b)+ ...

alpha=asin(sin(x(1)*pi/180)*sin(delta)+cos(x(1)*pi/180)*cos(delta)*cos(t*pi/180)); C=cot(alpha(:,1:20))./cot(alpha(:,2:21));

C=C';

end y=sum((C-D).^2);

%计算杆长的程序

function h = ganchang (phi_y,phi_x,N,t_b,l)

t_c=(120-phi_y)/15;%时差

st=t_b-t_c;%真太阳时

t=15*(st-12);%度数

b=2*pi*(N-1)/365;%弧度√

delta=(0.006918-0.399912*cos(b)+ 0.070257*sin(b)-...

0.00148*sin(3*b) );%弧度√

0.006758*cos(2*b)+0.000907*sin(2*b)-0.002697*cos(3*b)+ ...

alpha=asin(sin(phi_x*pi/180)*sin(delta)+cos(phi_x*pi/180)*cos(delta)*cos(t*pi/180)); h=abs(l'.*tan(alpha));%h为自变量

end

%第四问

%处理视频获取影子长度的程序

clear;clc;

close all ;

obj = VideoReader('Appendix4.avi' ); %输入视频位置 numFrames = obj.NumberOfFrames;% 帧的总数 for k = 1:4500:63000% 帧25 1s

frame=rgb2gray(frame); frame = read(obj,k);%读取第几帧

% imshow(frame);%显示帧

end imwrite(frame,[num2str((k-1)/4500),'.jpg' ],'jpg' ); % 保存帧

范文三：太阳影子定位问题的分析与研究

太阳影子定位问题的分析与研究

摘要

通过分析物体的影子长度随时间变化的规律来确定视频的拍摄地点和日期，是视频数据分析的重要方面，在实际生活中有重要的应用价值。

针对问题一，首先查阅资料获取了任意给定日期的太阳直射点的纬度计算公式，然后讨论了地方时与北京时间的补差问题，在此基础上建立了影子长度变化的数学模型，分析了影子长度关于杆长、日期、地方时、杆子位置的经纬度以及年份的变化规律，给出了题目指定日期和时间段的天安门广场上的杆子的影子长度的变化曲线。

针对问题二，首先讨论了影子变化规律与杆子地理位置的关系，我们注意到，影子长度由长变短应该发生在正午之前，影子长度由短变长应该发生在正午之后，由此可以在经度上限定位置的可能范围，而影子按顺时针方向应该在太阳运行面的北侧，影子按逆时针方向应该在太阳运行面的南侧，由此可以在经度上限定位置的可能范围；然后，计算了附件1的影子长度和影子转动方向，为确定位置限定了范围；接着，建立并求解了多项式拟合模型，初步确定杆子位置经度在东经111.02度附近；在此基础上建立最优匹配搜索模型，求解得杆子的最可能位置是北纬19.2944度，东经108.7237度，位于海南省。

针对问题三，首先对附件2和附件3的数据进行了预处理，计算了影子长度，确定了影子均为顺时针转动，限定了可能位置的经纬度；增加日期为未知参数，建立最优匹配搜索模型，得到了一系列不同日期、经纬度相近的可能位置点，附件2的可能位置点分布在新疆自治区，附件3的可能位置点分布在湖北省。

针对问题四，首先采用Free Video to JPG Converter软件对视频进行处理，每2分钟截图一次，利用MATLAB的IMtool工具箱测绘出杆和影子长度，计算出截图时刻的影子长度，分别已知日期和未知日期的最优匹配搜索模型，求解得，已知日期的可能位置位于内蒙古，未知日期的匹配日期有2组，7月份的位置在内蒙古，1月份的位置在马来西亚。

关键词：太阳影子 经纬度 最优匹配搜索

1. 问题重述

确定视频的拍摄地点和日期是分析视频数据的重要内容，分析视频中物体的太阳影子变化是其中一种方法。附件1-3是不同地点和时间下的影子顶点坐标数据，附件4为一根直杆在太阳下的影子变化的视频。我们要分析解决以下问题。

（1）问题一：根据太阳影子长度的影响因素，建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律。应用所建立的模型画出今年10月22日北京时间9时至15时之间天安门广场上3米高的直杆的太阳影子长度的变化曲线。

（2）问题二：根据水平地面上某固定直杆的太阳影子顶点坐标数据，建立数

学模型来推断该直杆所处的地点。利用所建立的模型来推断在给定日期和时间段形成附件1中直杆顶点影子坐标的直杆的若干个可能的地点。

（3）问题三：根据水平地面上某固定直杆的太阳影子顶点坐标数据，建立数

学模型来推断该直杆所处的地点和日期。利用所建立的模型来推断在给定时间段形成附件2及附件3中直杆顶点影子坐标的直杆的若干个可能的地点和日期。

（4）问题四：建立通过已知日期和时间段的视频确定视频拍摄地点的数学模

型，应用模型给出附件4中视频若干个可能的拍摄地点。假如拍摄日期未知，尝试确定附件4中视频可能的拍摄地点与日期。

2. 问题分析

问题需要建立影子变化的数学模型，寻求影子关于各个参数的变化规律，根据影子的变化确定地点和日期。

问题一，要求建立影长变化模型，分析影长与各个参数的变化规律。需要通过日期和时刻分别确定太阳直射点的纬度和经度，根据立杆位置的经纬度确定太阳高度角，建立影长变化的数学模型，进而分析影长与该模型各个参量的变化规律。使用控制变量法，改变其中的某一个参量，研究这个因素对事物影响。

问题二，要求在已知日期的条件下根据顶点坐标确定立杆地位置。需要求出影子的长度与转角，通过研究影长关于地理位置的变化规律，确定参考范围，利用多项式拟合模型拟合影长与时间曲线，利用影长最小值的时间求出当地经度的不精确值。再基于最优匹配模型对参考范围进行搜索求解。

问题三，要求在未知日期条件下确定立杆地位置和日期。需要求出地理位置的参考范围，在问题二模型的基础上增加日期参量，忽略年份因素，再利用最优匹配模型对参考范围进行搜索求解。

问题四，要求通过视频信息确定位置，讨论日期对求解的影响。需要忽略视频拍摄角度误差。每间隔一定时间，获取影长信息。利用最优匹配模型对参考范围进行搜索求解。要求对搜索算法不断优化。

3. 模型假设与符号说明

3.1模型假设

1. 不考虑太阳光线受折射、海拔高度等因素的影响。

2．到达地球上的太阳光视为平行光线。

3．水平地面是地球球面上过该点的切面。

4．摄像机向直杆的垂线视为垂直于直杆与影长所在平面。

5．附件4的标识时间为北京时间。

3.2符号说明

符号名称

y-m-d

t

s

说明 日期 地方时 影长 经度差 符号名称 D h H 说明

公历天数 杆长 太阳高度角 第i月的总天数

4. 模型的建立与求解

4.1问题一的模型建立与求解

4.1.1 太阳直射点纬度的确定

太阳直射点是地球表面太阳入射角度为90°的地点即日地连心线与地球表面的交点，太阳直射点所在经线的地方时为12时。通过对日期的计算，可以确定太阳直射点的纬度，便于接下来影长变化模型的建立与求解。

太阳直射点纬度 的计算公式 如下：

式中 称日角， ? ，

式中D为理论天数，即日期在年内的顺序号， ( ) ( ) ，其中Y为年份， ( )为向下取整。

4.1.2 对时间差的修正

地球是自西向东自转，东边时间总比西边早。同一时刻，地球上不同经度的时间不同，称为地方时。区时是某地所在时区的中央经线的地方时。区时经度每 15 度相差一小时，地方时经度每1 度相差 4 分钟。

对于附件一，北京时间是北京所在的东八区的区时，并不是北京（东经116.4°）地方时，它们满足：

北京时间-北京地方时=14.4分钟

4.1.3 影长变化的数学模型的建立

如图4-1，设地球半径为R， A地的经线交赤

道于D点，地方时为 ,某日太阳直射点在B地， B

地的经线交赤道于C点。则A地的纬度 的大小为

( ），太阳直射点B的纬度 的

大小为 ( )。

由地理时区知识 可得AB两点的经度差 为

( ) ( )

、 分别是太阳直射点B和立杆地点A的经

度。 ???? 与A地水平地面的法向量????? 夹角，大小等于 ，考虑实际，对A地 是阳光? 白昼时刻有：| | 。

O为原点，以OD及其延长线为X轴，地轴为Z轴建

， ， )， ， ， )，

图4-2是AOB平面截地球所得圆面，因为视太阳光线平行，显然：

因此可得

????? ， ????? ? ?

如图4-3，设在阳光照射下，A地有直杆AE，在水平地面上的影子为AF，阳光与AE的夹角为 ，则显然满足

( )

天文学中规定 ，H为太阳高度角，太阳

高度角H与影长s成负相关关系，当太阳高度角越大，影

长越小。

经过以上推导，建立影长变化的数学模型如下：

( ) ( )

( ) ( ) {

式中，s为影长,h为杆高，H为太阳高度角，太阳直射点坐标为（ ， ），直杆所在地坐标为（ ， ）， 为太阳直射点与直杆所在地的经度差，y为年份， 称日角，D为理论天数， ( )为向下取整。

4.1.4 影子长度变化规律的分析

通过建立影长变化的数学模型，我们可以分析杆长h、日期D、地方时t、当地经度 、当地纬度与影长关系。

1）杆长h与影长s的关系

这里， 可视为定值，一定为正值， 是自变量为 的正比例函数。随着 增加， 线性增加。

2）日期D与影长s 的关系

我们按照季节将全年划分成四个阶段，计算出该季节在当年的公历天数起止，进行分析，如表4-1所示：

表4-1四季时间划分依据表

夏天 6.25-9.23 季节 起止日期 公历天数起止

秋天

9.24-12.22 第267天-第356天

冬天 12.22-次年3.20 第356天-次年第79天

（1） 对南北回归线外向两极区域分析

基于影子变化的数学模型，利用matlab绘制出北京时间9:00时，天安门广场（北纬39度54分26秒,东经116度23分29秒），3米高的杆在太阳影子长度的曲线，横轴表示公历理论天数，纵轴表示影子长度。虚线依次代表春分、夏至、秋分、冬至日。

见图4-4。

由图联系地理知识可知，天安门广场位于北回归线以北。在冬天和春天，太阳直射点从南回归线向北回归线移动，影子长度随着日期的增加而逐渐减短，夏至日时，太阳直射点到达北回归线，影子长度最短；在夏天和秋天，太阳直射点从北回归线向南回归线移动，影子长度随着日期的增加而逐渐增长，冬至日时，太阳直射点到达南回归线，影子长度最长。

（2） 对南北回归线以内向赤道区域分析

我们以坐标为北纬3度,东经116度23分29秒地点为例，分析3米高的直杆，北京时间9:00时在该地点太阳影子长度的曲线，横轴表示公历理论天数，纵轴表示影子长度。虚线依次代表春分日、夏至日、秋分日、冬至日，见图4-5。 图4-5

该地点位于赤道与北回归线之间。在冬天，太阳直射点从南回归线向赤道移动，影长随着日期的增加而逐渐减短；在春天，太阳直射点从赤道经过该地点向北回归线移动，未经过该地点时，影长随着日期的增加而逐渐减短，经过该地点后，影长随着日期的增加而逐渐增长；在夏天，太阳直射点从北回归线经过该地点向赤道移动，未经过该地点时，影长随着日期的增加而逐渐减短，经过该地点后，影长随着日期的增加而逐渐增长。在秋天，太阳直射点从赤道向南回归线移动，影长随着日期的增加而逐渐增长。

3）地方时t与影长s的关系

基于影子变化的数学模型，绘制出2015年10月22日北京时间07：30-17：30时天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆在太阳影子长度的曲线，横轴表示时间，纵轴表示影长，见图4-6。

图4-6

由图可知，在07：30-17：30，影子长度随着时间先减小，后增大，在两端斜率趋势明显增大。10月22日时，太阳直射南半球，全球太阳从东偏南升起，西偏南落下，影长先减小后增大，在当地时为12点时，取得最小值。太阳未升起之前和已落下之后，没有阳光，无法产生影子，影长无确定值，趋向无穷。

4）当地经度 与影长s的关系

地球是自西向东自转，东边的时间也比西边的早。跨过180°经线时，规定日期应减少一天。据此，我们画出3米高的直杆处于天安门广场所在纬圈（北纬39度54分26秒）不同位置点（经度取值为东经90°向东至西经120°），在北京时间2015年10月22日09：00的曲线。

为了便于讨论我们分别画出当地时为2015年10月22日和2015年10月21日的曲线，坐标轴表示经度（‘+’表示东经，‘-’表示西经），纵轴表示影长，见图4-7

图4-7

由图可知，影长随着经度取值向东移动，先减短后增大。地球是自西向东自转，东边的时间比西边的早。在当地时为正午12点的地点(大约在东经165°)影长取得最小值。当地时接近日出时刻与日落时刻，影长增长趋势明显增大。图线与影长与时间图线近似。

5）当地纬度 与影长s的关系

基于影子变化的数学模型，绘制出2015年10月22日北京时间09：00时天安门广场所在经线（东经116度23分29秒）上南纬40°-北纬40°的地点3米高的直杆在太阳影子长度的曲线，如图：横坐标表示纬度，纵坐标表示影长。

图4-8

由图可知，在南纬40°-北纬40°范围内，影长随着地点北移先减小后增大。在南纬15°取得最小值。北半球秋季，太阳直射点离开赤道向南回归线移动。太阳直射点所在地的影长最短，同一经线上，越远离太阳直射点所在位置，在此时刻影长越大。

6）年份 与影长s的关系

我们选取一年中某一天为例，忽略一年内的影长s变化，分析年份Y与影长s的关系。基于影子变化的数学模型，画出1985年-2020年10月22日北京时间09：00时天安门广场（北纬39度54分26秒,东经116度23分29秒）上3米高的直杆在太阳影子长度的散点图，横坐标表示年份，纵坐标代表影长。

图4-9

由图看出，年份不同，影长也所变化。历法中年度天数与地球实际公转周期具有偏差，实际时间与理论时间有差，图示现象符合客观实际。

4.1.5 影长变化的数学模型的求解

结合问题一数据计算，在这里，日期是2015年10月22日，则D为295，天安门广场坐标为北纬39度54分26秒,东经116度23分29秒，直杆高为3m，北京时间（东八区，取东经120度）为9:00-15:00。利用matlab编写程序得到图：

图4-10

我们对9:00-15:00以20分钟为频度给出影长，给出各时刻的影长：

时间 11:00 12:00 13:00 14:00

影长 4.072983 3.688905 3.820433 4.509125

表4-1 时间影长关系表 时间 影长 11:20 12:20 13:20 14:20

3.883981 3.67657 3.981103 4.899798

时间 11:40 12:40 13:40 14:40

影长

3.757348 3.719983 4.207632 5.401753

可知：在9:00-15:00内，直杆影子长度先减小后增加，在12小时15分影长最小值取得3.674 米。

4.2问题二的模型建立与求解

4.2.1 影长转动角度与地理位置的关系 地球自西向东自转，在地球上看，太阳一天内东升西落，影子指向由西向东。但是太阳直射点的位置影响影子的转向，如图

影子长度由长变短应该发生在正午之前，影子长度由短变长应该发生在正午之后，由此可以在经度上限定位置的可能范围，而影子按顺时针方向应该在太阳运行面的北侧，影子按逆时针方向应该在太阳运行面的南侧，由此可以在经度上限定位置的可能范围。

E

图4-11

当地时为12时的太阳高度称为正午太阳高度。由天文常识可知，一天中正午时太阳高度最大，影长s最小。通过拟合曲线，确定影长s最小点的北京时间，可以得到立杆地点的经度。

利用附件1的 21组影子顶点坐标(x,y)，得到： 影长s：

√

以y轴为基线的角度 ：

√

结果如下表所示

表4-2 附件一影长参数表

14:48 14:51 14:54 14:57 15:00 15:03 15:06 15:09 15:12 15:15 15:18 15:21 15:24 15:27 15:30 15:33 15:36 15:39 15:42

1.1038 1.1383 1.1732 1.2087 1.2448 1.2815 1.3189 1.3568 1.3955 1.4349 1.4751 1.516 1.5577 1.6003 1.6438 1.6882 1.7337 1.7801 1.8277 0.5085 0.5142 0.5198 0.5255 0.5311 0.5368 0.5426 0.5483 0.5541 0.5598 0.5657 0.5715 0.5774 0.5833 0.5892 0.5952 0.6013 0.6074 0.6135 65.26536 65.69006 66.1037 66.5023 66.89424 67.27196 67.63757 67.99574 68.34385 68.68767 69.01821 69.34462 69.66153 69.9735 70.2804 70.57914 70.8719 71.1595 71.44476 1.215297 1.249051 1.283195 1.317993 1.353364 1.389387 1.426153 1.4634 1.501482 1.540232 1.579853 1.620145 1.661271 1.703291 1.746206 1.790051 1.835014 1.880875 1.927918

（1）纬度范围的估计

附件二给出该日日期为2015年4月18日，太阳直射点位于北回归线和赤道之间。

基于影子旋转方向与地理位置关系的讨论，通过分析附件一数据旋转角度与时间的关系，发现旋转角随着时间逐渐较小，即该点位于北半球。据此我们可知该点纬度范围在 。

（2）经度范围的估计

基于影子旋转方向与地理位置关系的讨论，通过分析附件一影长与时间的关系，发现影长随着时间逐渐增大（表4-2），发现当地时已过正午12点，该地在该时刻当地时为正午12点的地点以东，据此我们将经度范围选定在 。（‘+’表示东经，‘-’表示西经）

4.2.3 多项式拟合模型的建立

多项式拟合模型 是解决曲线拟合最常用的方法，我们采用该模型对影长与北京时间的函数进行拟合。

整体上考虑近似函数 ( )同所给数据点( )( )误差 ( )

( )的大小,常常利用误差平方和∑ 则多项式数据拟合的具 进行衡量，

体做法是：假定给定数据( )( )，求 ( ) ， 为所有次数不超过

( )的多项式构成的函数类， 现求 ( ) ∑ ，使得

∑(∑ )

为 的多元函数，因此上述问题即为求 ( )的极值问题。由多元函数求极值的必要条件，得

∑(∑ )

即

∑(∑ ) ∑

可以证明，上式的系数矩阵是一个对称正定矩阵，故存在唯一解。从上式中解出 ( )，从而可得多项式

( ) ∑

4.2.4 多项式拟合模型的求解

利用Excel软件,绘制21组北京时间与影长的数据的散点图。

图4-12

由图结合第一问的影长变化模型，图线可看作抛物线，开口向上，影长在当地时为12点处，取得最小值。故最高次m取2。

用matlab中的polyfit函数对影长和时间进行拟合，得到拟合函数

图4-13 图4-14

求出，在北京时间12.5984点时，当地时为12点。由 (当地时 北京时间）可得，当地经度大概在东经111.0236°。

4.2.5 最优匹配搜索模型的建立

（1）问题一参量的简化

基于问题一影子长度变化规律的分析，年份对影子长度的影响较小，忽略年份因素。我们观察到，太阳直射点纬度与时间的关系成正弦图像，我们将计算公式近似为

( )

() ( )

（2）最佳匹配目标函数的建立

基于影长变化数学模型， 我们建立高度、经度、纬度范围的最佳匹配模型：

∑( ? )

?

( )

( )

( ){

( , )

（3）最优匹配搜索模型的建立

最优匹配搜索模型是使用搜索方法取得最佳匹配模型目标函数的最优值。 考虑求解无约束最优化问题 ( ) ， 的下降迭代法。设 ( )连续可微，并利用它的一阶导数表达式。在每个迭代点 ( ),利用它的梯度向量 ( )( ( )),计算一个满足

( ) ( )

的下降方向 ( )，再沿方向 ( )进行搜索确定一个适当的步长 ( )，得新的迭代点

( ) ( ) ( ) ( )

为确保迭代序列{ ( )}， 能收敛于的一个局部最优解 ,步长 ( )的选取应使 ( ( ))充分小于 ( ( ))。当 ( )满足预先指定的精度要求时，迭代终止。并把 ( )接受为 的一个近似。

其主要工作量集中于线性搜索过程中确定 ( )所必须的函数值与导数值的计算。就理想情形而言，步长 ( )应取为函数 ( )在方向 ( )上的极小值，即

( ) ( ( ) ( ))

argminf(x)是f(x)达到最小值时的x的取值。 相应于求解一个非线性方程：

( ( ) ( )) ( )

4.2.6 最优匹配搜索模型的求解

基于多项式拟合模型的求解，得出当地经度大概在东经111°，且该地此时必须具有光照条件。据此，我们将经度范围选定在 。

我们在高度范围 ，经度范围 ，纬度范围 内基于最优匹配搜索模型编程进行搜索。

基于附件一长度精度为 ，目标函数的阈值应取在 ，允许系统误差导致的最低位偏差，利用matlab 调用fminsearch函数，得到符合条件的唯一解。

是

杆长 /m 利用经纬度地图，发现该点位于我国海南省，符合客观情况。

纬度/° 经度/° K

图4-16

4.3问题三的模型建立与求解

4.3.1 数据预处理 （1）对附件二数据的分析

表4-3 附件二影长角度参数表

北京时间 12:47 12:50 12:53 12:56 12:59 13:02 13:05 13:08 13:11 13:14 13:17 13:20 13:23 13:26 x坐标(米) -1.1813 -1.1546 -1.1281 -1.1018 -1.0756 -1.0496 -1.0237 -0.998 -0.9724 -0.947 -0.9217 -0.8965 -0.8714 -0.8464 y坐标(米) 0.2048 0.2203 0.2356 0.2505 0.2653 0.2798 0.294 0.308 0.3218 0.3354 0.3488 0.3619 0.3748 0.3876

角度

-80.1645 -79.1977 -78.2035 -77.1912 -76.1444 -75.0733 -73.9763 -72.8489 -71.6889 -70.4973 -69.2718 -68.017 -66.7269 -65.3951 影长

1.198921 1.175429 1.15244 1.129917 1.107835 1.086254 1.065081 1.044446 1.024264 1.00464 0.985491 0.96679 0.948585 0.930928

13:35 13:38 13:41

-0.7719 -0.7473 -0.7227 0.4246 0.4366 0.4484 -61.1861 -59.705 -58.1824 0.880974 0.865492 0.850504

基于影子旋转方向与地理位置关系的讨论，通过分析附件二影长与时间的关系，发现影长随着时间逐渐增大，如表…，发现当地时未过正午12点，该地在该时刻当地时为正午12点的地点以西，据此我们将经度范围选定在 。（‘+’表示东经，‘-’表示西经）

基于影子旋转方向与地理位置关系的讨论，通过分析附件二数据旋转角度与时间的关系，发现旋转角随着时间逐渐较小，即该点位于北半球。据此我们将纬度范围选定在 ( ) ,这里 ( )是影子比对日的太阳直射点的纬度。

（3） 对附件三数据的分析

表4-4 附件三影长角度参数表

北京时间 13:15 13:18 13:21 13:24 13:27 13:30 13:33 13:36 13:39 13:42 13:45 13:48 13:51 13:54 13:57 14:00 14:03 14:06

x坐标(米) 1.2791 1.3373 1.396 1.4552 1.5148 1.575 1.6357 1.697 1.7589 1.8215 1.8848 1.9488 2.0136 2.0792 2.1457 2.2131 2.2815 2.3508 y坐标(米) 3.3242 3.3188 3.3137 3.3091 3.3048 3.3007 3.2971 3.2937 3.2907 3.2881 3.2859 3.284 3.2824 3.2813 3.2805 3.2801 3.2801 3.2804 角度

21.04594 21.94681 22.84478 23.73783 24.625 25.50915 26.38613 27.2587 28.12479 28.98502 29.83868 30.68585 31.52718 32.36044 33.18778 34.00774 34.82086 35.62614

影长

3.561798 3.578101 3.595751 3.614934 3.635426 3.657218 3.680541 3.705168 3.731278 3.758918 3.788088 3.818701 3.85081 3.884585 3.919912 3.956876 3.995535 4.035751

基于影子旋转方向与地理位置关系的讨论，通过分析附件三影长与时间的关系，发现影长随着时间逐渐增大，发现当地时已过正午12点，该地在该时刻当地时为正午12点的地点以东，据此我们将经度范围选定在 。（‘+’表示东经，‘-’表示西经）

基于影子旋转方向与地理位置关系的讨论，通过分析附件二数据旋转角度与时间的关系，发现旋转角随着时间逐渐较小，即该点位于北半球。据此我们将纬度范围选定在 ( ) ,这里 ( )是影子比对日的太阳直射点的纬度。

4.3.2 模型的参数修正

问题三是在问题二的基础上缺少了日期参量，无法直接得到太阳直射点的纬度。 我们增加考虑日期参数，利用无约束最优化线性搜索模型进行求解。日期范围全年，视该年为平年，取全年365天，即 。

我们建立高度、经度、纬度、日期范围的无约束最优化线性搜索模型：

∑( ? )

?

( )

( )

( )

∑ {

( ， ( ) , , , )

4.3.3 模型的求解 （1）对附件二的求解

我们在高度范围 ，经度范围 ，纬度范围 ( ) 内基于无约束最优化线性搜索模型编程进行搜索。

基于附件二长度精度为 ，目标函数的阈值应取在 ，允许系统误差导致的最低位偏差，利用matlab 调用fminsearch函数，得到符合条件的解，如下表所示：

表4-5 问题三附件二求得经纬表

杆高/m 2.049543 2.055079 2.060361 2.065397 2.070173 2.074692 2.078951 2.080809 2.076689 2.072305 2.067664 2.062765 2.057622 2.052227 2.046588 2.040719 2.034607 2.028267 纬度/° 40.76755 40.86568 40.95922 41.0481 41.1323 41.21179 41.28652 41.31919 41.24685 41.1698 41.08807 41.0017 40.91071 40.81515 40.71507 40.61048 40.50146 40.38803 经度/° 80.02031 80.05044 80.07888 80.10584 80.13117 80.155 80.17732 80.18696 80.16548 80.14245 80.1179 80.09176 80.06418 80.03494 80.00413 79.97188 79.93788 79.9023 月份 6 6 6 6 6 6 6 7 7 7 7 7 7 7 7 7 7 7 日数 1 2 3 4 5 6 7 6 7 8 9 10 11 12 13 14 15 16

K

2.89E-08 3.17E-08 3.47E-08 3.77E-08 4.08E-08 4.40E-08 4.71E-08 4.85E-08 4.54E-08 4.23E-08 3.92E-08 3.61E-08 3.31E-08 3.02E-08 2.75E-08 2.50E-08 2.27E-08 2.07E-08 （2）对附件三的求解

我们在高度范围 ，经度范围 ，纬度范围 ( ) 内基于无约束最优化线性搜索模型编程进行搜索。

基于附件二长度精度为 ，目标函数的阈值应取在 ，允许系统误差导致的最低位偏差，利用matlab 调用fminsearch函数，得到符合条件的解，如下表所示：

表4-6 问题三附件三求得经纬表

杆高/m 2.989941 3.008068 3.006719 2.988776 2.9711 2.9537

纬度/° 32.53331 32.65778 32.64847 32.5255 32.40562 32.28886 经度/° 110.288 110.2709 110.2722 110.2891 110.3059 110.3224 月份 2 2 11 11 11 11 日数 3 4 8 9 10 11 K

3.48E-08 3.07E-08 3.10E-08 3.52E-08 4.13E-08 4.93E-08

利用经纬度地图，发现附件二和附件三地点分别集中位于我国新疆维吾尔族自治区和湖北省，符合客观情况。

图4-17 图4-18

4.4问题四的模型建立与求解

4.4.1 数据的准备

附件4为一根直杆在太阳下的影子变化的视频，时长40分钟。已知杆长h为2米。我们采用Free Video to JPG Converter软件对视频进行处理，每两分钟截图一次，获得21张图片，利用matlab中imtool工具箱测绘出视频中杆长与影长长度。由于时间较短，测量影子角度变化对影长误差，认为摄像机向直杆的垂线垂直于直杆与影子的平面。

实际 米已知，即可得出比例，从而求出影长。

实际 测量

结果如下表：

表4-7 问题四测量得影长

实际 测量

时间 8:57:06 8:59:06 9:01:06 9:03:06

测量 683.69 683.69 683.69 683.69

测量 787.64 776.66 767.12 758.09

实际

2.3041 2.272 2.2441 2.2176

9:09:06 9:11:06 9:13:06 9:15:06 9:17:06 9:19:06 9:21:06 9:23:06 9:25:06 9:27:06 9:29:06 9:31:06 9:33:06

4.4.2 问题的求解

683.69 683.69 683.69 683.69 683.69 683.69 683.69 683.69 683.69 683.69 683.69 683.69 683.69 729.60 720.04 711.05 703.55 694.04 683.53 674.02 664.52 658.01 645.51 642.26 637.72 626.34 2.1343 2.1063 2.08 2.0581 2.0303 1.9995 1.9717 1.9439 1.9249 1.8883 1.8788 1.8655 1.8322

利用获得的21组数据，基于无约束最优化线性搜索模型进行搜索：

∑( ? )

?

( )

( )

( ){ ( ， ( ) )

由于测量误差的存在，我们将K的阈值取到 ，得到一组数据:

图4-18

使用经纬度地图查询，发现该地点位于内蒙古自治区。

4.4.3 日期未知情况下的求解

假如日期未知，即缺少了日期参量。增加考虑日期参数，利用无约束最优化线性搜索模型进行求解。日期范围全年，视该年为平年，取全年365天，即 。

∑( ? )

?

( )

( )

( )

∑

{

( ， ( ) , , ) 由于测量误差的存在，我们将K的阈值取到 ，得到140组数据（见支撑材料answer442），将经纬度作图4-19，可发现部分点位集聚，我们在各自范围内取目标函数K最小的参数，得到两组数据，见图4-20和图4-21。

图4-19

图4-20

图4-21

经度 111.6°

纬度 41.54°

时间范围 7月

地点范围

中国内蒙古

5. 模型的评价与推广

5.1模型优点

1.建立了影长变化的数学模型，较为准确地确定了各参数与影长的定量关系，合理地对杆高，经纬度，日期，时间，年份与影长关系进行定性分析。

2.采用最小二乘法拟合模型，利用多项式拟合方法对影长与时间变化曲线进行了拟合，得到了经度的不精确值，使得搜索范围得到了优化。

3.建立了优化匹配模型，利用无约束最优化线性搜索方法对地理范围进行了搜索，算法优化度较高，同时求得的解较为精确，达到 。

5.2模型缺点

1.在对日期的求解过程中忽略了年份对太阳直射点纬度的变化影响。 2.没有考虑真太阳时与平太阳时的误差。

3.没有考虑摄像机拍摄角度导致的影长测量误差。

5.3模型改进

1.优化匹配模型增加年份参量，考虑年份变化对影长的影响 2.修正太阳自转对时角的影响。

3.通过坐标系转换，消除视角内实际影长与视频影长的误差

5.4模型推广

本文利用影长信息求解立杆点位置与日期，是定位技术中一种比较方便的方法，对易于测量的数据求得了较为精确的结果。应用前景较为广阔。本文只是粗浅地涉及了有关太阳影子定位技术的一点边缘，今后对于这种定位计算还应当有更加简洁使用的方法。

6. 参考文献

[1]百度百科，太阳直射点，http://baike.baidu.com/view/93613.htm

[2]汪和平，影子与季节_纬度的关系，中学数学，2008.09

[3]邵雷，最小二乘法的基本原理和多项式拟合，数学参考，2006.02

[4]赵家奎，无约束最优化线性搜索方法的研究，武汉工业大学学报，1995.12

7. 附录

7.1问题一探究日期与影长关系的程序

clc clear

spheroid=referenceEllipsoid('earth');

[xa,ya,za] = enu2ecef(0,0,0,39.9072254,116.3913982,0,spheroid); j=1; n0=79.6764+0.2422*(2015-1985)-fix((2015-1985)/4);

for i=1:1:365

d=i-n0;

wd=0.3723+23.2567*sin(2*pi*d/365.2422)+0.1149*sin(2*2*pi*d/365.2422)-0.1712*sin(3*2*pi*d/365.2422)-0.758*cos(2*pi*d/365.2422)+0.356*cos(2*2*pi*d/365.2422)+0.0201*cos(3*2*pi*d/365.2422); jd=120+(12-9)*15;

[xb,yb,zb] = enu2ecef(0,0,0,wd,jd,0,spheroid);

ra=xa^2+ya^2+za^2; rb=xb^2+yb^2+zb^2;

rab=(xa-xb)^2+(ya-yb)^2+(za-zb)^2; p=acos((ra+rb-rab)/2/sqrt(ra)/sqrt(rb));

h=3;

l(j)=tan(p)*h; t(j)=i; j=j+1; end plot(t,l) hold on

plot([80,80],[0,16],'--') plot([175,175],[0,16],'--') plot([266,266],[0,16],'--') plot([356,356],[0,16],'--') axis([ 0 365 0 16])

7.2问题一影长变化图线的绘制程序

clc clear

n0=79.6764+0.2422*(2015-1985)-fix((2015-1985)/4); d=295-n0;

wd=0.3723+23.2567*sin(2*pi*d/365.2422)+0.1149*sin(2*2*pi*d/365.2422)-0.1712*sin(3*2*pi*d/365.2422)-0.758*cos(2*pi*d/365.2422)+0.356*cos(2*2*pi*d/365.2422)+0.0201*cos(3*2*pi*d/365.2422); spheroid=referenceEllipsoid('earth');

[xa,ya,za] = enu2ecef(0,0,0,39.9072254,116.3913982,0,spheroid); j=1;

for i=9:0.01:15 jd=120+(12-i)*15;

[xb,yb,zb] = enu2ecef(0,0,0,wd,jd,0,spheroid); ra=xa^2+ya^2+za^2; rb=xb^2+yb^2+zb^2;

rab=(xa-xb)^2+(ya-yb)^2+(za-zb)^2; p=acos((ra+rb-rab)/2/sqrt(ra)/sqrt(rb)); h=3;

l(j)=tan(p)*h; t(j)=i; j=j+1; end plot(t,l) grid on

7.3问题二多项式拟合程序 clc clear

point=xlsread('附件1-3','附件1','B4:C24'); x=point(:,1); y=point(:,2); z=sqrt(x.^2+y.^2); t=(14.7:0.05:15.7)'; xishu=polyfit(t,z,2); syms x

df=diff(0.148902881082785*x^2-3.75188316236736*x+24.1275007161760); f=inline('0.148902881082785*x^2-3.75188316236736*x+24.1275007161760'); x=solve(df); y=f(x); tt=double(x); jd=120+(12-tt)*15;

7.4最优匹配搜索模型的求解程序

function f=problem2fun(x)

global year month day shadowlength currenttime residual=0; for i=1:21

residual=residual+(shadow(year,month,day,currenttime(i),x(1),x(2),x(3))-shadowlength(i))^2; if(residual>0.01) residual=10000; break; end end

f=residual; end

%计算影长

function [shadow] =shadow(year,month,day,time,h,latitude,longitude) n0=79.6764+0.2422*(year-1985)-fix((year-1985)/4); T=ahargana(year,month,day)-n0; theta=2*180*T/365.2422;

delta=0.3723+23.2567*sind(theta)+0.1149*sind(theta*2)-0.1712*sind(theta*3)-0.758*cosd(theta)+0.3656*cosd(theta*2)+0.0201*cosd(theta*3); t=120+(12-time)*15-longitude;

H=asind(sind(latitude)*sind(delta)+cosd(latitude)*cosd(delta)*cosd(t)); shadow=h/tand(H); end clc clear

data_problem2_1=xlsread('???t1-3','???t1','B4:C24'); global year month day shadowlength currenttime

shadowlength=sqrt(data_problem2_1(:,1).^2+data_problem2_1(:,2).^2); year=2015;month=4;day=18; currenttime(1)=14+42/60; for i=2:21

currenttime(i)=currenttime(i-1)+0.05; end i=1; for h=1:1:3

for lati=-90:1:90 for longi=90:1:130

[x,fval,exitflag]=fminsearch(@problem2fun,[h,lati,longi]); if exitflag==1

answer2(i,:)=[x(1),x(2),x(3),fval]; i=i+1; end end

end end

for i=1:length(answer2) if answer2(i,4)

answer22(i,:)=answer2(i,:); end end

answer22=unique(answer22,'rows');

7．5计算影长程序

function [shadow] =shadow(year,month,day,time,h,latitude,longitude) n0=79.6764+0.2422*(year-1985)-fix((year-1985)/4); T=ahargana(year,month,day)-n0; theta=2*180*T/365.2422;

delta=0.3723+23.2567*sind(theta)+0.1149*sind(theta*2)-0.1712*sind(theta*3)-0.758*cosd(theta)+0.3656*cosd(theta*2)+0.0201*cosd(theta*3); t=120+(12-time)*15-longitude;

H=asind(sind(latitude)*sind(delta)+cosd(latitude)*cosd(delta)*cosd(t)); shadow=h/tand(H); end

7.6计算积日程序

%?????yè?￡¨ahargana￡?

function [ahargana] = ahargana(year,month,day) if mod(year,4)==0 if month==1 ahargana=day; elseif month==2 ahargana=31+day; elseif month==3 ahargana=60+day; elseif month==4 ahargana=91+day; elseif month==5 ahargana=121+day; elseif month==6 ahargana=152+day; elseif month==7 ahargana=182+day; elseif month==8 ahargana=213+day; elseif month==9 ahargana=244+day; elseif month==10 ahargana=274+day;

elseif month==11 ahargana=305+day; else

ahargana=334+day; end

else

if month==1

ahargana=day; elseif month==2 ahargana=31+day; elseif month==3 ahargana=59+day; elseif month==4 ahargana=90+day; elseif month==5

ahargana=120+day; elseif month==6

ahargana=151+day; elseif month==7

ahargana=181+day; elseif month==8

ahargana=212+day; elseif month==9

ahargana=243+day; elseif month==10 ahargana=273+day; elseif month==11 ahargana=304+day; else

ahargana=334+day; end

end

end

30

范文四：定位问题

定位问题方法

随着it 行业的发展，软件测试的普及，再加上测试行业入行门槛要求相对较低，在各大小城市中，软件测试工程师如雨后春笋般地涌现；在这个竞争激励的年代，想要在测试行业有所突破，需要掌握扎实的理论知识外，还需要掌握不断发展和改进的测试技能，以下为我在测试工作中对软件测试定位问题方法做的一些总结，欢迎大家多批评指正。

在软件测试过程中，发现问题，可能不是什么难事，但要精准的定位到问题出在哪里，需要用到一定的技能和经验积累；

而研发人员，往往也都希望测试人员对问题能做到精准的描述，以便能快速定位到问题所在模块，进行快速分析和解决；

一、抓包分析

1、在使用浏览器测试过程中，如何判断是前端的问题，还是后端的问题，单单从报的错误，或者无响应，或者返回空白页面上，都无法定位到问题出在哪里，这个时候，就需要借助一些工具来定位问题，首先从发起的http 请求开始定位；

2、基于浏览器的http 请求，许多浏览器都有自带的插件，可以抓取到http 的请求和返回，例如chrome ，firefox 都有自带的抓包工具，开启之后，都可以看到你手工操作工作中抓取到的http 请求和返回报文；如下图为chrome 浏览器按F12抓取到进行搜索的请求和返回。

3、除了自带的插件可以抓包外，还有许多第三方的抓包工具，例如fiddler ，以下简单介绍一下fiddler 抓取的报文；

打开fiddler ，按F12开启抓包；

对点击收取邮件进行操作，进行报文抓取

抓取到的http 请求和返回数据如下：

4、这就是一个简单的获取邮件列表时抓取到的http 请求和返回数据，其他的请求进行类似的操作即可抓取到，如何判断是http 请求的报文有问题，还是响应结果有问题，就需要根据接口设计规范文档中的描述来作为依据，进行判断；

5、fiddler 可以模拟发起请求，通过Composer 页签，把获取到的请求分别填写到请求url ，请求body ，点击Execute 即可（需要保证当前的会话有效）

二、看日志分析

排除掉前端的问题之后，需要逐层分析，从后台服务来进行分析定位，当然，定位分析问题的前提条件是，需要对业务流程非常清楚，清楚各模块之间的交互；假如后台服务的平台为linux ，则需要对linux 的

一些常用命令进行掌握，采用tail -f 、 more , grep , 等命令来跟踪日志，最终定位到具体模块的错误信息，或者接口请求、返回结果

三、查询数据库分析

在测试过程中，除了通过前端的体现来对用例的执行结果进行判断外，还需要对数据库的结果来进行检查，例如对于一个订单的提交，前端返回成功，如果单凭前端返回成功来判定该用例执行成功的话，很有可能在后台数据处理出现错误，该用例实际的执行结果是失败的； 这种情况下，就需要对数据库表中的值进行查看验证。

四、排除法定位

我们发现问题时，往往直接提交bug ，但有些问题是在特定账号或者特定环境下才会出现，而当研发对问题进行修改时，发现问题不能重现而打回，为了能对问题进行精准定位，则需要用排除法来定位问题的重现条件。

范文五：会展定位问题

会展应有鲜明的主题,没有主题的会展是不能够吸引观众的,招展也就无法开展。

根据会展内容的不同,可以将会展市场分为:

① 以某种高科技产业或优势产业为依托举办的专业性科技博览会、交易会;

② 将某些产业与内外贸易相结合而开展的产品交易会、展销会;

③ 以宣传本地人文资源如文化、艺术、体育等为宗旨举办的博览会、展示会;

④ 以重要的城市为中心举办的综合性的国际会议及大型的博览会、展销会。会展公司在确定会展主题时,应充分 发挥本地的资源优势,使会展呈现出鲜明的主题,提高会展的竞争力。

明确的主题是会展明确细分市场, 也就是确立会展的受众目标 —— 参展商和观众的关键, 是会展突出个性特点的标 志。要确定成功的主题,会展公司应该从以下几个方面人手:

(1) 选择优势领域

所谓优势领域,就是选择在巾国较为发达或颇具发展前景的领域或行业。展会的主题定位于这些优势领域,才能 充分显示展会对于相关行业发展的带动作用和展会显著的营销效应。目前,我国在机床、信息、汽车和环保技术等领 域具有特定的优势和广阔的市场发展前景。

例如,北京地区举办的大型国际展览会中,机床、通信、纺机、印刷、冶金、汽车等专业性展览会已经进入世界 先进会展的行列,其中被全球会展权威机构 —— 国际博览联盟 (UIF)认可的品牌展览会就有 6个。同时,根据欧洲信息 技术监理会发布的全球信息通信技术市场报告, “ 早在 1998年, 中国的 IT 市场就已经超过韩国, 到 2003年可能要超过 日本成为亚太地区最大的市场。中国电信业也已成为世界最具发展潜力的市场,并且至少在 21世纪的第一个 l0年里, 这种地位将不会被动摇。 ” 因此,当德国汉诺威会展公司准备进军中国市场时,就将展会的主题定位于这些领域。 (2) 突出区域特色

一个地区会展品牌与其本身的经济与产业发展特点是密切相关的。会展公司要创建名牌展会,就应结合自身条件, 利用区域优势,塑造有区域特色的展会。

例如, 以时尚产品著称于世的法国巴黎, 正是因时装、化妆品等成功展会而使其享有 “ 会展之都 ” 的美称;而被誉为 “ 购物天堂 ” 的香港也是以珠宝、 皮草、 玩具等会展著称。 深圳近几年高新技术产业异军突起,目前高新技术产业产值已 占到国民生产总值的 60%。深圳会展主题的选择就紧紧围绕深圳市今后的产业发展方向,突出高新技术城市和新兴城 市的鲜明特点,打造独具特色的城市会展品牌,增强了区域特色和吸引力,带动了相关产业的发展。深圳高交会的成 功走的正是这条道路。

(3)体现专业性质

众所周知,专业性已成为世界范围内会展发展的一大趋势,从 l992年开始,世界综合性的大型博览会已由专业化 的博览会代替。

例如,汉诺威工博会虽是综合展,但却是由若干个专业展组成的,如机器人展、自动化立体仓库展、铸件展、低 压电器展、灯具展、仪器仪表展、液压气动元件展等。这些专业展一般两年办一次,这样尽管工博会年年办,但下面 的各个专业展主题却不重复,而且每个专业展的规模和水平均居世界一流,成为各个行业的名牌展会。专业展会与博 览会相比,针对性更强,展览项目不易重复,能够更加深入地促进行业贸易的发展,充分体现展会专业合作与交流的 渠道作用。

(4)拓展发展空间

目前国内大大小小、形形色色的会展主题已不计其数,要在会展主题的选择上实现创新,已经很难找到空白的会 展主题。因此,会展公司可以与原有的展会合作,通过展会主题的收购与兼并来拓展市场,增强展会的竞争力。已被 国际展览联盟认可,与德国杜塞尔多夫冶金铸造展和美国克里夫兰钢铁展并列为世界三大冶金及热加工展览会之一的 中围国际冶金工业展就是一个成功的典范。它的迅速崛起在很大程度上就得益于主办者通过展会主题的收购与兼并, 塑造出国际化品牌的会展。德国杜塞尔多夫展览公司进入中国市场,也准备采取与中国同行合作办展的方式来拓展市 场。因此,积极拓展原有展会的空间,不仅可以有效扩大展会规模,有利于形成品牌展会,同时也是解决重复办展问 题的有效途径

转载请注明出处范文大全网 » 太阳影子定位问题