2015年A题一等奖论文

太阳影子定位问题

摘要

目前，如何确定视频的拍摄地点和拍摄日期是计算机视觉的热点研究问题，是视频数据分析的重要方面，有重要的研究意义。本文通过建立数学模型，给出了通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的方法。

对于问题一，建立空间三维直角坐标系和球面坐标系对直杆投影和地球进行数学抽象，引入地方时、北京时间、太阳赤纬、杆长、太阳高度角等五个参数，建立了太阳光下物体影子的长度变化综合模型。求解过程中，利用问题所给的数据，得到太阳赤纬等变量，将太阳赤纬等参量代入模型，求得了北京地区的9:00至15:00的影子长度变化曲线，当12：09时，影子长度最短；并分析出影长随这些参数的变化规律，利用控制变量法思想，总结了五个参数与影子长度的关系。最后进行模型检验，将该模型运用于东京、西藏两地，得到了这两座城市的影长变化规律曲线，发现变化规律符合实际两地实际情况。

对于问题二，为了消除不同直角坐标系带来的影响，将实际坐标转换为二次曲线的极坐标，建立了极坐标下基于多层优化搜索算法的空间匹配优化模型。求解时，先将未知点的直角坐标系的点转换为极坐标，然后设计了多层优化搜索算法，通过多次不同精度的搜索，最后得出实际观测点的经纬度为东经E 115︒北纬N 25︒。同时对模型进行验证，实地测量了现居住地的某个时间段的值，通过模型二来求解出现居住地的经纬度，分析了误差产生的原因：大气层的折射和拟合误差。

对于问题三，将极坐标转换后的基本模型转换为优化模型，建立了基于遗传算法的时空匹配优化模型。将目标函数作为个体的适应度函数，将经度纬度及日期作为待求解变量，用遗传算法进行求解，得到可能的经度纬度及其日期：北纬20度，东经114度，5月21日；北纬20度，东经114度，7月24日；东经94.5度，北纬33.8度，6月19日。最后，将遗传算法与多层优化搜索算法进行对比分析，得出遗传算法的求解效率和求解精度均优于多层次搜索算法。

对于问题四，首先将视频材料以1min 为间隔进行采样得到41帧（静态图片），将这些静止图片先利用matlab 进行处理，后进行阀值归一化处理，得到这些帧的灰度值矩阵。在图片上建立参考模型，获得影子端点的参考位置。利用投影系统和模型二，建立了基于图形处理的视频拍摄地点搜索模型。利用模型二中多层搜索算法，求得满足精度的最优地点。最优的地点是：东经119, 北纬48.7，在内蒙古的呼伦贝尔市。同时假设日期是未知量，将模型四与模型三相结合，得到了可能的地点和时间，并分析了可能出现误差的原因，最后回答了当视频日期未知，也可以确定其位置和日期。 最后，给出了模型的优缺点和改进方案。

关键词：极坐标化，多层优化搜索算法，遗传算法，图像处理，MATLAB

1. 问题重述

1.1问题背景

随着现代科技的发展，日常生活中摄像机的应用越来越普遍。无论是个人家庭还是组织单位，都通过摄像机来录制各种视频以分享信息，例如实时视频监控、记录自然景观、观测气象信息等。而通过视频来确定拍摄地点的地理位置信息是目前计算机视觉领域的热点研究问题之一。一个视频的地理位置能够提供当地气候、平均温度、平均降雨量、植物索引、地表概况、海拔高度和人口密度等大量背景信息[1]。因此从视频中确定地理位置是一项有很大潜力应用空间的技术。

1.2问题描述

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

试建立数学模型讨论下列问题：

1. 建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用所建立的模型画出2015年10月22日北京时间9:00-15:00之间天安门广场3米高的直杆的太阳影子长度的变化曲线。

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

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

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

2. 问题分析

2.1问题一分析

问题一要求分析投影长度随各参数的变化规律，建立影子长度变化的数学模型。首先对直杆建立空间三维坐标系，将地球简化成规则球体建立球面坐标系。在这两个坐标系中，通过几何证明，运用向量知识可分析出影响影子长度的各种参数，得出地球上某日白天某时刻影子顶端在地平面上的具体位置，由此可以给出影子长度的变化规律。

2.2问题二分析

问题二要求根据某固定直杆在水平地面上的太阳影子顶点坐标数据及日期数据，建立数学模型确定直杆所处的地点。与第一问有相似之处，但分析附件所给数据，发现附

件中只给出x 、y 坐标值，而并没有给出xy 轴的准确方向，所以考虑将直角坐标转换成极坐标，来消除由于不同坐标系选取所造成的影响。 2.3问题三分析

问题三与问题二有相似处，区别是第三问附件没有提供日期，需要根据直杆影子端点坐标确定直杆所在地点的经纬度和日期。具体的日期可以由太阳直射点纬度来确定，而根据问题二中的模型，xy 坐标与太阳直射点纬度有关。如果继续用第二问的模型来求解，需要不断改变太阳直射点纬度来拟合极坐标方程，这样做算法复杂度会很大。所以考虑对问题二模型进行修改，不采用拟合，而直接建立与待求点经纬度以及日期有关的目标函数，通过约束经纬度范围来缩小待求点的可行域，从而简化算法复杂度。

2.4问题四分析

问题四中，直接以视频的方式给出了固定杆长的距离变化规律。将图片形式的影长变化规律以坐标的形式进行转换，转换为现实的坐标形式。这样就可以利用问题二的模型，整合现有的算法，求出拍摄地点。

3. 模型假设与符号系统

3.1模型的假设

（1）假设地球为一个规则的球体。

（2）由于日地距离远大于地球半径，所以假设太阳光线为平行光。 （3）假设地球上某地的水平地面是地球球面上过该地的切面。 （4）假设不考虑太阳光线穿过大气层时所发生的折射。 （5）假设一天中太阳直射点的纬度不变。

（6）假设不考虑太阳的视面角、高山阻挡、海拔高度等因素的影响。 （7）假设不考虑阴天没有阳光的情况。

3.2符号系统

问题一符号系统

符号 意义 单位

直杆所在地纬度值 度 α

β 太阳直射点的纬度 度

A 、B 两地经度差 度 θ

ϕ 太阳光线与直杆的夹角 度

直杆长度 米 h

直杆影长 米 L

t 地方时 时 t 0 北京时间 时

E 直杆所在地的经度 问题二、三符号系统

意义

直杆所在地纬度值 太阳直射点的纬度 附件1中第一组坐标的y 值

极径 直杆长度 极角 问题四符号系统 意义 固定杆长度

实际长度与灰度值坐标下的转换比例

投影系统

度

符号 α β

单位 度 度 米 米 度

y 1 ρ h θ

符号 L k

P

单位 米

4. 问题一的建模与求解

4.1问题分析

在问题一中，为了描述直杆影子长度变化的动态过程，首先以直杆为z 轴，建立空间三维坐标对直杆影子的变化进行数学抽象。再将地球作为规则球体建立球面坐标系，利用空间解析几何与平面解析几何的知识，对两个坐标系中的相关向量与角度进行分析，分析出影响影子长度的参数，得到影子端点在坐标系中的位置表达式。由此可以求出影子长度随各个参数的变化规律。

建模流程图如下所示：

图4.1问题一建模流程图

4.2模型准备

为了建模的方便，先给出一些地理名词的解释和一些数据的预处理方法。

4.2.1名词解释[6]

地方时：以一个地方太阳升到最高的地方时间为正午12时，将连续两个正午12时之间等分为24个小时，所成的时间系统。它是观测者所在的子午线的时间。

北京时间：是中国采用北京时区的区时作为标所在的东八准时间。北京时间并不是

北京（东经116.4°）地方的时间，而是东经HF 120°地方的地方时间。

太阳赤纬：是地球赤道平面与太阳和地球中心的连线之间的夹角。

太阳直射点：地球表面太阳光射入角度(即太阳高度角 )为90度的地点，它是地心与日心连线和地球球面的交点。

太阳高度角：对于地球上的某个地点，太阳高度角是指太阳光的入射方向和地平面之间的夹角；专业上讲是指某地太阳光线与通过该地与地心相连的地表切线的夹角。

4.2.2数据预处理

（1）经纬度转换

在问题一中，天安门广场的坐标是用经纬度（度分秒）的形式给出的。为了下面建模求解的方便，将其统一转换成以“度”为单位。

换算方法为：分位数除以60，秒位数除以3600。 所以，天安门广场的纬度可以转换为：

39︒54'26''=39+54÷60+26÷3600=39.907︒ 经度可以转换为：

116︒23'29''=116+23÷60+29÷3600=116.391︒

（2）北京时间与地方时的转换[9]

问题中所给出的时刻为北京时间，而北京时间指的是东经120°地方的地方时，并不是问题中地点的地方时。所以先要将所给的北京时间转换成相应的地方时。

转换规则为：东经度120度地区，每增加1度，加4分钟。 所以有转化公式：

⎧t 0+(E-120)*4,E>120t =⎨

t -(120-E)*4, E

其中，E 表示直杆所在地点的经度，t 0是北京时间，t 是直杆所在地方的地方时。

用此公式对问题一中的北京时间进行操作，得到直杆所在地的地方时，如下表所示：

要研究影子的变化，需要建立空间三维坐标对直杆影子的变化进行数学抽象。通过对直杆和地球分别建立了两个空间直角坐标系，用空间解析几何和向量知识，可以确定两个坐标系上各点之间的位置和角度关系。

4.3.1建立直杆处空间三维坐标系

根据假设，视太阳光线为平行光，以直杆所在地点的正东方向为x 轴，以正北方向为y 轴，以直杆直立即垂直于x 0y 平面的方向为z 轴，建立空间直角坐标系，得到直杆在x 0y 平面的投影与光线的位置关系，如下图所示：

图4.2直杆空间三维坐标系

其中，AE 是与过A 处的经线相切的方向向东的单位向量；AK 是A 处地平面内方向向北的单位向量。AH 是A 处垂直于xOy 平面的直杆，AF 是该直杆在xOy 平面内的投影，HF 是当天太阳光线的照射方向，照射方向与直杆所成角度∠FHA =ϕ。

4.3.2建立直杆在地球上的宏观空间球面坐标系

根据假设，可视地球为规则球体O 1，过直杆底端A 处的经线与赤道交于D 点，B

点为某日的太阳直射点，过B 点的经线与赤道交于C 点。

以O 为原点，以OD 所在直线为x 轴，以地轴ON 所在直线为z 轴建立空间直角坐标系

O -xyz ，如图4.3所示：

图4.3 直杆在地球上的空间三维坐标系

4.3.3确定各点之间的位置和角度关系

（1）α为直杆所处位置的纬度数，并且-90 ≤α≤90 。 若A 地在北半球，则α>0，若A 地在北半球，则α

（2）亦即上面提到的赤纬，并且-23 26' ≤β≤23 26' 。 β为太阳直射点B 地的纬度，

（3）θ为A 地与B 地的经度差，t 是地方时。 对于某日A 地白昼t 时刻：θ=(12-t ) ⨯15 (0≤t ≤24) 。

（4）∠AOB =∠AHF ，证明过程如下：

由假设可知，太阳光线是一簇簇的平行线，所以HF //BO ，如图4.4，圆O ’是过A,B 两地的大圆，于是∠AOB =∠AHF ，证毕。

由以上分析可得：

设地球半径为R ，∠AOD =α, ∠BOC =β, ∠DOC =θ, ∠AOB =ϕ1，则有：

⎧AE =(0,1,0),AK =(-sin α,0,cos α) ⎪

⎨A (R cos α,0, R sin α), B (R cos βcos θ, R cos βsin θ, R sin β)

⎪

⎩cos ϕ=cos =cos αcos βcos θ+sin αsin β

图4.4 过A 、B 两地，以地球中心为圆心的圆

4.3.4确定日影坐标及其长度

（1）确定影子端点的横坐标

AH h

=如图4.4，在Rt ∆AHF 中，HF =，其中h 为直杆长度。 cos ϕcos ϕ

δ设HF 与AE 所成角为，则cos δ=cos =cos =-cos βsin θ

cos δ

h ，即F 点在平面如图4.2，对HF 在AE 上的正投影AJ ，有AJ =HF cos δ=

cos ϕ

A -xy 上投影端点的横坐标：

-cos βsin θ

x =h （4.1）

cos αcos βcos θ+sin αsin β

（2）确定影子端点的纵坐标

AK γ设HF 与成角为，cos γ=cos =cos =sin αcos βcos θ

cos γ

h ，即点-cos αsin β，如图1，对HF 在AK 上的正投影AG ，有AG =HF cos γ=

cos ϕ

F 在A -xy 上投影端点的纵坐标：

sin αcos βcos θ-cos αsin βy =h （4.2）

cos αcos βcos θ+sin αsin β

（3）确定日影坐标的长度

已知直杆投影端点的横纵坐标，并且直杆底端即为坐标原点，所以可以得到直杆影长：

L =4.3）

4.3.5影子长度变化的综合模型

根据上面的分析，太阳光下物体影子的长度变化综合模型为：

⎧L =⎪

-cos βsin θ⎪

x =h ⎪⎪cos αcos βcos θ+sin αsin β

（4.4） ⎨

⎪y =sin αcos βcos θ-cos αsin βh ⎪cos αcos βcos θ+sin αsin β⎪ ⎪⎩θ=(12-t ) ⨯15(0≤t ≤24)

4.4模型的求解

4.4.1赤纬角β求解

太阳赤纬角β在每一年的任何时刻的值都是可求的，其计算公式为[7]： ⎧β=0.3732+23.2567sin θ1+0.1149sin 2θ1-0.1712sin 3θ1+0.758cos θ1⎪+0.3656cos 2θ+0.0201cos 3θ

11⎪

⎪2πT

（4.5） θ=⎨1

365.22⎪

⎪T =N -N 0⎪

⎩N 0=79.6764+0.2422*(year -1985) -floor((year-1985) /4)

其中式中θ1为日角，即θ1=2πt /365.2422；N 为积日，即日期在年内的顺序号，如平年12月31日为365，闰年的12月31日是366。year 为计算时刻所在的年份，floor 为向下取整函数。

根据问题一中的2015年10月22日，可知积日N=295，year=2015，所以可以求出：

T =215.0576，θ1=3.6996(弧度）

根据上述所求结果，得到：

β=-10.8636（度）

4.4.2直杆所在地与太阳直射点之间的纬度差θ的求解

纬度差计算公式有：

θ=(12-t )*15(度)

其中t 为直杆所在地的地方时。

将4.2.2中由北京时间转换出的地方时t 代入以上公式，可以得到不同时刻，直杆所在地点与太阳直射点的纬度差θ的变化值，如下所示：

4.4.3影长变化的求解结果

由于在很短时间内，影子不会出现大的变化，所以可以认为1分钟内，影子长度是近似不变的。将这段时间分为361个时间段，每一分钟是一个小时刻，将这个时刻的影长作为这一分钟内的影子长度。

将上面计算出来的β和θ代入影子端点的坐标和影子长度表达式，得到每一分钟，平面直角坐标系内影子端点的坐标变化值和影子长度变化值。由于数据较多，这里只给出每隔30分钟的数据样点，结果如下表所示：

由上表，可以作出天安门广场3米高的直杆在太阳下影子长度的变化曲线，如下所示：

图4.5直杆影子长度随时间的变化曲线

结论：

（1） 直杆的影长从9时开始，先减小，减小至北京时间12:09时，影长达到最

短，为3.673731米，之后开始增大。

（2） 10月22日北京正处于秋末，太阳直射点在赤道和南回归线之间，此时正

午时分直杆的影长比其本身更长。

（3） 北京时间12:00的影长为3.679733米，比12:09时稍长，这也进一步说

明北京时间并不是指示北京的地方时。

4.5分析影子长度和各参数之间的变化规律

问题中要求分析影子随各参数的变化情况，首先，根据4.3.5中的模型，可以看出影长L 和β、θ、α有关。而赤纬β是关于日期的函数，θ是关于地方时t 的函数，t 又是关于经度的函数。

所以综上可知，影响影子长度的参数有：直杆所在地的经纬度、地方时、当前的日期。

以影子长度与纬度的变化关系为例，研究直杆同一时刻同一经线上不同纬度地点的影长变化，将β，θ均视为定值，设：

k 1=cos βsin θ，k 2=cos βcos θ，k 3=sin β

则影子端点的坐标为：

-k 1⎧

x =h ⎪cos αk 2+sin αk 3⎪

（4.6）⎨

⎪y =sin αk 2-cos αk 3h ⎪cos αk 2+sin αk 3⎩

为了简明地表达二者之间的关系，取时刻为当地时间12点，即θ=0，取日期为问题所给10月22日时的太阳直射点赤纬，即β=-10.86 ，则：

k 1=0，k 2=0.982，k 3=-

0.188

所以影子端点的坐标为：

⎧x =0⎪

15.95⎨

y =⎪0.982-0.188tan α ⎩

此时影子的长度为：

15.95

0.982-0.188tan α

由此可以作出影子长度随纬度变化的变化趋势，如下所示：

L =

图4.6东经120度上影子随纬度的变化规律图

结论：

（1） 在东经120度上，直杆影子长度随着纬度的增加而逐渐增加，在纬度近似

为N 76︒时，影长开始陡增，在北纬79︒达到一个远大于正常情况的极大值，越过此极大值之后，影长又开始陡减，在纬度近似为N 82︒时减少速率逐渐平缓。

（2） 如下图所示：当太阳直射点纬度不是0度时即不直射赤道时，影子最长点

会出现在小于北纬90︒的某个纬度处，并且此时的影长接近无穷长，这就是图中在北纬79︒出现一个极高峰值的原因。

图4.7日照光线示意图

其他因素以此为例，进行同样的分析，就可得到各因素与影长的变化关系，正午影长随日期的变化如下图所示：

图4.8 天安门广场午时影长随日期变化规律图

结论：

（1） 午时天安门广场的影长随日期的变化规律为：一年中从第一天开始随着日期

的变化，影子长度先减小，达到一个最小值，再增大。

（2） 2015年天安门广场午时影长最短的一天是一年中的第173天。

4.5模型检验

将问题一模型运用到2015年的10月22日的其他城市。在这里，取西藏（东经91.11，北纬29.97）和东京（东经138.6，北纬35.5）为检验的对象。

由上面的模型，计算出在西藏和东京，一根3m 长的直杆在太阳下得到的影子长度随北京时间的变化曲线分别是：

图4.9 西藏在北京时间9:00至15:00的影子长度变化曲线

图4.10 东京在北京时间9:00至15:00的影子长度变化曲线

结论：

（1） 西藏的地方时比东八区（东经120度）区时晚2小时左右，所以西藏的正

午时间为北京时间14:00左右，模型规律与实际的影长曲线规律是相符的。 （2） 东京的地方时比东八区（东经120度）时间早1小时左右，所以东京的正

午时间为北京时间11:00左右，模型规律与实际的影长曲线规律是相符的。

5. 问题二的建模与求解

5.1问题分析

问题二中，附件给出的仅仅是直杆所处地平面上未知x 轴方向和y 轴方向的坐标值。

为了消除观测者在观测时任意选定坐标轴造成的影响，

对题目所给的坐标数据进行平移

处理，再将直角坐标转换成极坐标，给出影子端点轨迹的极坐标方程。再根据附件中的数据，对含有参数的极坐标方程进行拟合，得出相关参数值。对于每一个确定经纬度和日期的观测点，代入极坐标方程可以得到相应函数值。以该函数值最接近0为目标，建立基于多层优化搜索算法的空间匹配优化模型。

建模流程图如下所示：

图5.1问题二建模流程图

5.2模型的准备

由于附件中并没有给出直杆的原长，所以需要先对直杆的长度进行估算，下面直杆长度的计算需要用到正午太阳高度角的概念。正午太阳高度角H 指的是一天中最大的太阳高度角。计算公式如下所示：

H =90︒-|β-α|

其中，β为太阳直射点的纬度，α为直杆所在地的纬度。

5.3模型的建立

通过第一问求得的影子端点坐标，得到影子端点的直角坐标系下的轨迹方程，再建立极坐标系，将处理过后的xy 坐标转换成极坐标，给出极坐标下的轨迹方程。

5.3.1确定影子端点的轨迹方程

由模型一可知影子端点的横纵坐标表达式为：

-cos βsin θ⎧

x =h ⎪

cos αcos βcos θ+sin αsin β⎪

⎨

sin αcos βcos θ-cos αsin β⎪y =h ⎪cos αcos βcos θ+sin αsin β⎩

移项代入化简得：

y sin αsin β+h cos αsin β⎧

cos βcos θ=⎪h sin α-y cos α⎪

（5.1） ⎨

x sin β⎪-cos βsin θ=

⎪h sin α-y cos α⎩

两边平方消去θ，可以得日影端点F 在直杆底端所在平面A -xy 上的轨迹方程为：

x 2sin 2β-y 2cos(α+β)cos(α-β) +hy sin 2α-h 2sin(α+β)sin(α-β) =0（5.2）

5.3.2将直角坐标转换为极坐标

为了消除不同地点观测者选取坐标轴方向的随机性对坐标产生的影响，将直角坐标转换成具有统一极点和极轴的极坐标系。

（1） 对原先直角坐标进行预处理

对于附件1中21组xy 坐标的数据，保持他们横坐标不变，纵坐标都减去附件1中第一组坐标的纵坐标值，即：

⎧x ' =x

（5.3）⎨

y ' =y -y ⎩1这样的处理相当于对投影端点轨迹曲线作了平移，将第一组数据平移到了x 轴上，

如下图所示：

图5.2对附件1中坐标进行预处理

（2）对影子端点建立极坐标系

以直杆所处位置底端为极点，以极点到第一组影子端点连线方向为极轴，建立极坐标系，如下图所示：

图5.3对影子端点建立的极坐标系

（3）将模型一中求得的轨迹方程转换成极坐标方程

应用上述规则，预处理之后的坐标为：x '=x , y '=y -y 1，其中x , y 代表原先的21组数据值。

为了将轨迹方程转换成极坐标方程，令x '=ρsin θ, y '=ρcos θ，所以有：

x =ρcos θ, y =ρsin θ+y 1，（5.4）

代入式（4.8），求得影子端点轨迹的极坐标方程f (ρ, θ) =0：

ρ2cos 2θsin 2β-ρ2sin 2θcos(α+β)cos(α-β) +ρsin θ(sin2α⨯h -2y 1cos(α+β)cos(α-β)) +

hy 1sin 2α-h 2sin(α+β)sin(α-β) -y 12cos(α+β)cos(α-β) =0

为了便于观察，将极坐标中变量的系数用常数符号表示，得到了确定日期确定经纬度以及确定杆长的固定直杆的影子端点随时间变化轨迹的极坐标方程，如下：

k 1ρ2cos 2θ-k 2ρ2sin 2θ+k 3ρsin θ+k 4=0（5.5）

其中：

⎧k 1=sin 2β⎪

⎪k 2=cos(α+β)cos(α-β)

（5.6）⎨

k =sin 2α⨯h -2y cos(α+β)cos(α-β) 1⎪3

⎪k =hy sin 2α-h 2sin(α+β)sin(α-β) -y 2cos(α+β)cos(α-β) ⎩411

5.3.3对已给数据进行极坐标二次曲线拟合

由以上推导，待求位置固定直杆在所在地平面上的投影的极坐标方程形式已知。所以可以用极坐标系下的二次曲线模型对已知xy 坐标的21组数据进行拟合：

A ρ2cos 2θ-B ρ2sin 2θ+C ρsin θ+D +ε=0

其中A , B , C , D 为回归系数，ε为随机误差，服从均值为0的正态分布，ρ, θ为21组平面直角坐标经过坐标转换之后得到的极坐标值，拟合精度为10-7

拟合出的极坐标方程为：

W (ρ, θ) =A 0ρ2cos 2θ-B 0ρ2sin 2θ+C 0ρsin θ+D 0=0( 5.7)

5.3.4建立直杆地点的空间匹配优化模型

（1）确定目标函数：

n

其中ρi , θi 为每一搜索点处的极坐标值。

min F i =

∑(ρ, θ)

i

i

i =1

n

(i =1, 2,..., n ) (5.8)

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

n

⎧

(ρi , θi ) ∑⎪

⎪min F i =i =1(i =1, 2,..., n )

n ⎪

⎪W (ρ, θ) =A ρ2cos 2θ-B ρ2sin 2θ+C ρsin θ+D =0

0000

⎪⎪

⎨x =ρcos θ, y =ρsin θ+y 1

(5.9)⎪-cos βsin θ

⎪x =h

cos αcos βcos θ+sin αsin β⎪⎪sin αcos βcos θ-cos αsin β⎪y =h

cos αcos βcos θ+sin αsin β⎪⎩

5.4模型的求解

为了求解上述模型，设计如下多层优化搜索算法： Step1，先对搜索点杆长的近似估计，得到合理的杆长；

Step2，再确定搜索点处影子端点极坐标，作为判断所搜点是否可行的依据； Step3，由陆地经纬度范围确定搜索点的论域；

Step4，以5度为步长进行顶层搜索，搜索出下层搜索的搜索点论域； Step5，以更小步长进行下层步长进行搜索，再次更新搜索点论域；

Step6，是否满足精度要求，不满足精度要求，再次进行step5，否则，当前值是最优值。

5.4.1搜索点杆长的近似估计

由于附件中没有给出具体杆长，所以在每一个搜索点，可以根据5.2中所提到的正午太阳高度角给出杆长的近似估计。

以东经90度北纬30的搜索点P 17,25为例，杆长计算步骤如下： （1）根据搜索点经度将所给数据北京时间信息转换成当地时间，t 0=14.7是附录所给北京时间，可以算出当地的地方时为：

t k =t 0-(120-Lo (i )) /15=14.7-(120-90)/15=12.7

（2）对所给横纵坐标(x (k ) , y (k ) )(k =1,2,...,21) 用最小二乘法进行拟合，求得横纵坐标随时间推移的线性关系：

y =ax +b

拟合结果如下：

y =0.147x +0.3475

（3）确定当地时间为正午12点时的影子长度

可将影子横纵坐标均视为随着时间的推移线性变换，所以得到当地时间正午时刻影子端点的横纵坐标：

(1)(21)(1)(1)

⎧⎪x 0=x -(x -x )(t k -12)

(5.10) ⎨

y =ax +b ⎪0⎩0

其中t

(1)k

为附录数据的第一组的当地时间。

l =直杆正午时刻的影长l ：

以P 17,25点为例，正午时刻的影长为：

⎧x 0=1.0365-(1.8277-1.0365)(12.7-12) =0.4827⎨

⎩

y 0=0.147⨯0.4827+0.3475=0.4185

l =0.6389

（4）根据正午时刻太阳高度角近似估计杆长

h i , j =tan(90-β-La (j ))

以P 17,25点为例，估算出的杆长为：

h 17,25=tan(90-10.3686-30) =1.1763

5.4.2确定搜索点处影子端点极坐标

首先根据式[8]得到影子端点的横纵坐标值：

-cos βsin θk ⎧

x (k ) =h i , j ⎪i , j

cos αcos βcos θ+sin αsin β⎪k

(5.11)⎨

sin La (j )cos βcos θ-cos La (j )sin βk ⎪y (k ) =h i , j i , j

⎪cos La (j )cos βcos θk +sin La (j )sin β⎩

其中La (j ) 为搜索点纬度，β为该日赤纬角在本题中为定值，θk 为观测点与当地时刻t k 太阳直射点的经度差：θk =(12-t k ) ⨯15 (0≤t k ≤24; k =1,2,...,21)

再依据坐标转换规则得到影子端点的极坐标值：

⎧ρ(k ) =i , j ⎪⎪

(5.12)⎨y i , j (k ) -y i , j (1)

⎪θi , j (k ) =arctan

x i , j (1)⎪⎩

以P 17,25点为例，第五组的极坐标值计算公式如下：

⎧ρ(5)==1.173

i , j ⎪⎪⎨y i , j (5)-y i , j (1) θ(5)=arctan =0.024⎪i , j

x i , j (5)⎪⎩

以此类推得到了P 17,25点的21组极坐标值，如下表所示：

表5.1

P 点的21组不同时刻极坐标值

5.4.3搜索点论域的确定

求解以上模型的关键是搜索方法的简化，经过查找资料，得到了七大洲大致的陆地经纬度范围：

5.4.4求解结果

在每一大洲的经纬度范围内，取五度为一步长，确定了每一个大洲搜索点论域，进行初步搜索

初步计算结果如下所示：

经过顶层搜索，确定了待求位置的取值区间。将取值区间更细划分，取1度为步长，进行下层搜索来进一步细化搜索过程。

下层搜索计算结果：

北纬N 25︒，东经E 115︒（具体地点为江西赣州）

5.5模型验证

为了验证该模型的正确性，我们进行了实地测量。

取40厘米长的直杆，于14:30至15:30在现居住地（E113，N30）进行了影长的坐标采样，得到了相关的数据，数据见附录。

利用模型二对此实测数据进行求解，得出的结果为：E115，N25

模型验证的结论：

（1）和实际的地点存在出入，但是误差相对较小。

（2）误差来源：未考虑太阳折射的误差、拟合曲线的误差，实地测量的误差。 （3）改进算法可使误差减小。

6. 问题三的建模与求解

6.1问题分析

问题三与问题二的主要区别在于是否已知日期，第三问需要根据所给的xy 坐标，求出地点和日期。日期可以通过太阳直射点的纬度求得，但在第二问中每次拟合时，都要不断改变太阳直射点纬度，非常繁琐。所以考虑修改模型，不采用拟合，而是直接构造与坐标点有关的、分辨度更高的目标函数。

6.2模型的建立

在模型二基础上进行修改，根据问题二的模型可知：

-cos βsin θ⎧

x =h ⎪cos αcos βcos θ+sin αsin β⎪

x =ρcos θ, y =ρsin θ+y 1⎨

sin αcos βcos θ-cos αsin β⎪y =h ⎪cos αcos βcos θ+sin αsin β，⎩

由于直杆所在地与太阳直射点的经度差用θ表示，会与极角混淆，所以这里用θ1表示，有：

-cos βsin θ1⎧

⎪x =cos αcos βcos θ+sin αsin βh ⎪1

（6.1）⎨

sin αcos βcos θ-cos αsin β1⎪y =h ⎪cos αcos βcos θ+sin αsin β⎩1

可以知道，在搜索过程中，每一个确定时刻已知经纬度的地点都可以通过上述式子求出其极坐标值。

将这些极坐标数据进行标准化处理，与附件中给出的确定时刻的影子端点真实坐标值进行比较，构造分辨度更高的目标函数。

(1)极坐标数据标准化

对一个经纬度确定的搜索点相应时刻的21组极坐标值ρk 进行线性标准化处理：

ρk -ρmin ⎧ρ' =⎪k

ρmax -ρmin ⎪

，k =1, 2,...,21（6.2） ⎨

⎪θ' =θk -θmin

k ⎪θmax -θmin ⎩

其中ρmin 和θmin 是附件2或3中所有数据极坐标中的最小值；ρmax 和θmax 是附件2或3中所有数据极坐标的最大值。

(2)构建目标函数

根据上面分析，问题三的目标函数为：

全国大学生数学建模竞赛一等奖论文

21

其中ρ0k , θ0k 为附件2或3中真实数据的极坐标值。

min d k =

k =1

21

（6.3）

建立出的时空匹配优化模型如下所示：

21

⎧⎪

⎪min d k =k =1

21⎪

⎪ρk -ρmin ρ' =⎪k

ρmax -ρmin

⎪⎪

⎪θ' =θk -θmin

（6.4） ⎨k θ-θ

max min

⎪

⎪x =ρcos θ, y =ρsin θ+y 1⎪

-cos βsin θ1

⎪x =h ⎪cos αcos βcos θ1+sin αsin β⎪

⎪y =sin αcos βcos θ1-cos αsin βh ⎪cos αcos βcos θ1+sin αsin β⎩

6.3模型的求解

由于第三问中新增加了一个决策变量日期，若再使用搜索算法，会加大搜索难度，加深算法复杂度，搜索时间会非常长，所以选择采用遗传算法求解第三问。

遗传算法实现的流程图如下：

图6.1 遗传算法流程图

全国大学生数学建模竞赛一等奖论文

6.3.1确定个体的适应度函数

参考第二问求解过程中各坐标值的确定过程，根据式（6.3）中目标函数与待求解变量的定量关系，求得目标函数，作为遗传算法中的个体适应度函数，即

F =d k =

k =21

21

（6.6）

6.3.2使用遗传算法对问题进行求解

利用遗传算法对附件2和附件3中数据进行求解，通过MATLAB 编程得到一些可能的时空坐标，如下表所示：

6.4

两种算法的对比分析

基于第二问的搜索方法，对其中参数之间的耦合关系进行修正，依然能用搜索算法对问题进行求解，将两种方法的求解过程进行对比，结果如下： （1）搜索时间对比

得到5-10组可能数据所需时间对比如下表所示：

结论：遗传算法求解效率更高

（2）搜索精度对比

两种算法五组目标函数的值结果对比如下所示：

全国大学生数学建模竞赛一等奖论文

结论：遗传算法拟合精度更高

7. 问题四的建模与求解

7.1问题分析

问题四中，直接以视频的方式给出了固定杆长的距离变化规律。此时需要将图片形式的影长变化规律以坐标的形式进行转换，同时转换为现实的坐标形式。这样就可以利用问题二的模型，求出拍摄地点。

7.2数据的预处理

影长变化规律需要以实际坐标形式进行展示，所以需要对视频资料进行预处理操作。

7.2.1将视频格式的数据转换为图片格式的数据

视频展示的是一个动态过程，而确定视频中影子端点的坐标需要在静态画面中进行。根据我国视频采用的PAL 制[1]，每分钟的画面帧数为25帧，影子在一分钟内的变化情况是十分微小的，所以，选取1分钟作为时间间隔来截取视频的静态图片。

以1分钟为时间间隔截取视频画面，得到的静态图片有41张，每张图片的分辨率为1920*1080。利用matlab 软件，读取这些图片的灰度值，以灰度值的变化来刻画杆子的影长的变化规律。

7.2.2确定影子端点在图片中的坐标

（1）将得到的静态图片灰度值矩阵预处理

为了消除静态图片中其他灰度值对刻画影长变化规律的影响，首先将静态图片的灰度值进行归一化，将图片的灰度值均与同一变量进行比较，在这里归一化的方法是阀值[3]归一化方法：

⎧0, h i -h 1

⎩255, h i -h 1>H （7.1）

H =10

其中H i 是归一化的灰度值，h i 是第i 张静态图片的灰度值，h 1是第一张图片的灰度

值，H 为事先所选的合适阀值，在这里选取灰度值为10为阀值。这种归一化方法消除了静态点的灰度值微小变化产生的影响。同时放大了属于影子长度变化灰度值的变化规律。

（2）确定端点在图片中的坐标位置

为了方便的刻画端点处的变化规律，将计算机处理灰度值的范围限制在端点处，在这里通过二维网络搜索方式[3]，搜索到应该将图片的灰度值坐标限制在（800~1000,800~1000）处。

全国大学生数学建模竞赛一等奖论文

7.3模型的建立

在7.2中求出了端点在图片中的坐标，为了找出视频可能的拍摄地点，需要知道端点在实际情况下的坐标。然后利用问题二的模型建立基于图形处理的确定视频拍摄地点模型。

全国大学生数学建模竞赛一等奖论文

7.3.1实际坐标与图片参考坐标的关系

实际坐标与图片中参考坐标系的投影关系如下[2]：

(x , y ,1) T =P ⋅(X,Y,Z,1) （7.4） P =[p 1, p 2, p 3, p 4]（7.5）

其中(x , y ) T 是所求点在参考坐标系中的坐标，[p 1, p 2, p 3, p 4]是投影矩阵，(X,Y, Z) 是实际坐标。

其中[p 1, p 2, p 3]是消隐点矢量[2]。p 4是照相机中心在参考坐标系的投影，此时为：

p 4=[1.935,3.440,0]T

7.3.2消隐矢量的求解

图7.2消隐点示意图

消隐点如图所示，需要求出每个消隐点所在直线，在这里选取固定杆底座作为参考物，求得本平面的消隐点矢量。选择消隐点的直线上的两个点，即三个消隐点有六个点，其灰度值坐标公式：

e ={(-62，01，) ,(43,01，), （43，01，) ，(43，，211) ，(43，01，),(59,4,1）}

求消隐点矢量的方法是：利用上面同一消隐点所在直线的端点，求出每个端点坐标叉乘e 1⨯e 2的值分别为：{a i , b i , c i }。

然后构建M 矩阵：

⎧a *a a *ba *c⎫⎪⎪

M =∑⎨b*ab*bb*c⎬（7.6）

i =1⎪⎪⎩a *cb*cc *c ⎭

求出M 对应的最小特征矢量即为消隐点矢量。

在求解上述的投影系统方程式之后，利用上述过程就可以得到实际的坐标。

n

7.3.3建立基于图形处理的视频拍摄地点搜索模型

模型的二的求解方法是通用的，在已知实际坐标的情况下，利用模型二对地点进行搜索，建立了基于图形处理的确定视频拍摄地点模型。

全国大学生数学建模竞赛一等奖论文

⎧⎪⎪T (x , y ,1) =P ⋅(X,Y, Z,1) ⎪⎪

⎨P =[p 1, p 2, p 3, p 4]

（7.7）⎪n

⎪(ρi , θi ) ∑⎪min F =i =1

(i =1, 2,..., n ) i ⎪n ⎩

7.4模型的求解

Step1求解实际空间坐标

利用上述所求的参考坐标、投影矢量P 和实际坐标和参考坐标之间的投影公式，可求解至各个点的实际坐标。

利用归一化处理静态图片上的静态点，得到影子端点的实际坐标如下表所示：

Step2求解视频拍摄地点坐标

利用模型二的多层次搜索方法进行求解，此处不再赘述。求解得出的结果为： 北纬N 48.7︒，东经E 119︒（具体地点为内蒙古呼伦贝尔市）

7.5对附加问题的分析

假设当时间是未知的，将第四问模型与第三问模型相结合，来求解第四问。首先将问题四转换为优化模型，确定适应度函数。

利用问题三的求解过程，求出以下的4个可能点和可能时间：

结论：

（1） 与已知的数据相比较，是可以求出的大致的时间范围； （2） 所求的数据与真实数据存在一定误差。

全国大学生数学建模竞赛一等奖论文

（3） 日期均相近，说明遗传算法在求解过程中上陷入了局部最优解。

问题四中提出的问题，“若视频拍摄日期未知，能否根据视频确定出拍摄地点与日期？”可以从充分性和必要性两个方面作出回答。

充分性：根据问题四处理数据的过程，在任何时间内均可以使用此种方法进行数据的预处理，经过时间的预处理，得到实际坐标，将影子变化规律转换为实际坐标的变化规律。此时，将坐标带入问题三中进行求解即可。

必要性：根据图像的影子变化规律在极坐标归一中，规律是存在的，总有一个时刻某个地点的影子是与之相同。同时我们设计的算法中，对模型的求解分为初步搜索和细化搜索，有比较强的适应性。

8. 模型的优缺点及改进

8.1模型优缺点 优点：

1、 设计的多层优化搜索算法可实现性高，能够满足用户的要求 2、 利用极坐标坐标变换很好的消除了二维的直角坐标的取向性 3、 将模型三转换为优化模型，并且利用遗传算法进行求解。 缺点：

1、 模型二算法虽然使用，但是比较复杂，且最优点可能掉进海里； 2、 模型三的优化模型在求解上可能陷入局部最优解；

3、 模型四比较复杂，有多个过程，且转换坐标可能不是很好符合实际。

8.2模型的改进

可以设计一种可实现的智能算法，防止算法得到的优化解陷入局部最优解。

[1]武琳，基于太阳阴影轨迹的经纬度估计技术研究，硕士学位论文，天津大学，2010年12月；

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

[3]侯洁，P2P 网络搜索算法研究[M]，天津：天津师范大学，2008； [4]七大洲经纬度 百度百科

http://wenku.baidu.com/link?url=PnYcSDQSFR1jg1JCKdRnR8oDTPyDU0ftJebyBFXpsXL0CLakssjJyaGV5WSb2iqEXQ3Yc9fRtiHCT3rOzQULxHmnzzQJltXQq2meSIFc-dW,9月13日 [5]电影平均一分钟多少个画格

http://zhidao.baidu.com/link?url=AxvEQlMw_--0IPKxlhFFzN_pCMm1aAI55P9MM01rPaRCxT9etS8MuzL6ZN8B74a_qz9xsDE8Whof6gUZL2nWG_,9月13日 [6]太阳高度角, 百度百科,

http://baike.baidu.com/link?url=4ZCmRKvAbUCIekYHSaUYwtoqzdEE_fwLzdtiOCP1Aex

全国大学生数学建模竞赛一等奖论文

0P7F0qiYzGwjRBdb_cZR4,9月11日 [7]太阳赤纬, 百度百科,

http://baike.baidu.com/link?url=OHC2VAYavh0Ggqw48aEnVJh_vnWfCVGaxCK2ntuISw7yqQ9L2CdRMvS0eRu_duSG,9月11日

[8]汪和平, 探究日影运动轨迹, 中学数学月刊,2010年第九期,29`30页； [9]北京时与地方时,

http://baike.baidu.com/link?url=OHC2VAYavh0Ggqw48aEnVJh_vnWfCVGaxCK2ntuISw7yqQ9L2CdRMvS0eRu_duSG,9月11日

全国大学生数学建模竞赛一等奖论文

附录

clc,clear;

beta=chiwei(2015,295); alpha=39.9; weidu=116.4; t1=9:30/60:15;

t=dangdishi(weidu,120,t1); n=length(t1); h=3;

for i=1:n

theta(i)=(12-t(i))*15; end

for i=1:n

x(i)=-cosd(beta)*sind(theta(i))/(cosd(alpha)*cosd(theta(i))*cosd(beta)+sind(alpha)*sind(beta))*h; end

for i=1:n

y(i)=(sind(alpha)*cosd(beta).*cosd(theta(i))-cosd(alpha)*sind(beta))/(cosd(alpha)*cosd(theta(i))*cosd(beta)+sind(alpha)*sind(beta))*h; end

for i=1:n

l(i)=(x(i)^2+y(i)^2)^0.5; end

% plot(t,y) plot(t1,l)

xlabel('时间/时'); ylabel('影子长度/m');

function t=dangdishi(beta1,beta2,t) if beta1>beta2

t=t-(beta2-beta1)/15; end

if beta1

t=t-(beta2-beta1)/15; end

function beta=chiwei(year,n)

n0=79.6764+0.2422*(year-1985)-floor((year-1985)/4); t=n-n0;

theta=2*pi*t/365.2422;

全国大学生数学建模竞赛一等奖论文

beta=0.3723+23.2567*sin(theta)+0.1149*sin(2*theta)-0.1712*sin(3*theta)-... 0.758*cos(theta)+0.3635*cos(2*theta)+0.0201*cos(3*theta);

function [theta,rough]=ji(x,y) theta=atan(y./x); a=length(x);

rough=zeros(a,1); for i=1:a

rough(i)=(x(i)^1/2+y(i)^1/2)^1/2; end

rough=rough;

第二问：

clc,clear;

global jingdu1 weidu1 beta=chiwei(2015,295); jingdu=[]; load Qa.mat x=Qa(:,2); y=Qa(:,3); weidu=[]; for i=0:5:160 jingdu(i/5+1)=i; end

for j=-80:5:80 weidu(j/5+19)=j; end

% l=xianxin(110,x,y); % h=tand(90-(beta+20))*l for i=1:33

l=xianxin(jingdu(i),x,y); for j=1:33

h(i,j)=tand(90-(beta+weidu(j)))*l; yizhi(jingdu(i),weidu(j),h(i,j)); end end

clc,clear; load Qa.mat x=Qa(:,2); y=Qa(:,3);

[theta0,rough0]=ji(x,y);

load jingdu1.mat load weidu1.mat g=length(jingdu1); e=length(x);

[theta,rough]=chaa(jingdu1,weidu1); theta1=theta; rough1=rough; for i=1:g

theta(i,:)=(theta(i,:)-min(theta(i,:)))/(max(theta(i,:))-min(theta(i,:)));

rough(i,:)=(rough(i,:)-min(rough(i,:)))/(max(rough(i,:))-min(rough(i,:))); end

for i=1:e

theta0(i)=(theta0(i)-min(theta0))/(max(theta0)-min(theta0)); rough0(i)=(rough0(i)-min(rough0))/(max(rough0)-min(rough0)); end

f=zeros(1,g); for i=1:g for j=1:e

f(i)=f(i)+((theta(i,j)-theta0(j))^2+(rough(i,j)-rough0(j))^2)^(1/2); end end f=f/21;

weidu=weidu1(find(f/21==min(f/21))); jingdu=jingdu1(find(f/21==min(f/21)));

clc,clear; for i=1:5

weidu(i)=15+(i-1)*1; jingdu(i)=105+(i-1)*1; end

for i=1:5

f(i)=ff(jingdu(i),weidu); end

i=find(f==min(f)); jingdu1=105+(i-1); j=ff1(jingdu1,weidu); weidu1=15+j-1;

function f=fun(x,y) f=0;

load A.mat

load B.mat load C.mat load D.mat load E.mat load F.mat C=0;D=0; for i=1:21

f=f+(A*x(i)^2*cos(y(i))^2+B*x(i)^2*sin(y(i))^2+... C*x(i)^2*sin(y(i))*cos(y(i))+D*x(i)*cos(y(i))+... E*x(i)*sin(y(i))+F)^2; end f=f/21; y=Qa(:,3);

[theta0,rough0]=ji(x,y); % load jingdu1.mat % load weidu1.mat g=length(weidu); e=length(x); for i=1:e

theta0(i)=(theta0(i)-min(theta0))/(max(theta0)-min(theta0)); rough0(i)=(rough0(i)-min(rough0))/(max(rough0)-min(rough0)); end

[theta,rough]=ch(jingdu,weidu); for j=1:g

theta(j,:)=(theta(j,:)-min(theta(j,:)))/(max(theta(j,:))-min(theta(j,:)));

rough(j,:)=(rough(j,:)-min(rough(j,:)))/(max(rough(j,:))-min(rough(j,:))); end

f=zeros(1,g); for i=1:g for j=1:e

f(i)=f(i)+((theta(i,j)-theta0(j))^2+(rough(i,j)-rough0(j))^2)^(1/2); end end f=f/21;

f2=f(find(f/21==min(f/21)));

function i=ff1(jingdu,weidu) load Qa.mat x=Qa(:,2); y=Qa(:,3);

[theta0,rough0]=ji(x,y);

% load jingdu1.mat % load weidu1.mat g=length(weidu); e=length(x); for i=1:e

theta0(i)=(theta0(i)-min(theta0))/(max(theta0)-min(theta0)); rough0(i)=(rough0(i)-min(rough0))/(max(rough0)-min(rough0)); end

[theta,rough]=ch(jingdu,weidu); for j=1:g

theta(j,:)=(theta(j,:)-min(theta(j,:)))/(max(theta(j,:))-min(theta(j,:)));

rough(j,:)=(rough(j,:)-min(rough(j,:)))/(max(rough(j,:))-min(rough(j,:))); end

f=zeros(1,g); for i=1:g for j=1:e

f(i)=f(i)+((theta(i,j)-theta0(j))^2+(rough(i,j)-rough0(j))^2)^(1/2); end end f=f/21;

i=find(f/21==min(f/21));

function [theta1,rough]=chaa(jingdu1,weidu1) load Qa.mat x=Qa(:,2); y=Qa(:,3);

g=length(jingdu1); t1=14.7:1/20:15.7; n=length(t1);

beta=chiwei(2015,108); for i=1:g

b(i)=xianxin(jingdu1(i),x,y);

h(i)=abs(tand(90-(beta+weidu1(i)))*b(i)); end

for i=1:g

t=dangdishi(jingdu1(i),120,t1); % b=xianxin(jingdu1(i),x,y); % h=tand(90-(beta+weidu1(i)))*b; % h=fun2(jingdu1(i),x,y); alpha=weidu1(i);

for j=1:n

theta(j)=(12-t(j))*15;

x(j)=-cosd(beta)*sind(theta(j))/(cosd(alpha)*cosd(theta(j))*cosd(beta)+sind(alpha)*sind(beta))*h(i);

y(j)=(sind(alpha)*cosd(beta).*cosd(theta(j))-cosd(alpha)*sind(beta))/(cosd(alpha)*cosd(theta(j))*cosd(beta)+sind(alpha)*sind(beta))*h(i); end

[theta1(i,:),rough(i,:)]=ji(x,y); end

function yizhi(jingdu,weidu,h) global jingdu1 weidu1 t=14.7:1/20:15.7; alpha=weidu;

beta=chiwei(2015,295);

t1=dangdishi(jingdu,120,t); n=length(t1); for i=1:n

theta(i)=(12-t(i))*15; % end

% for i=1:n

x(i)=-cosd(beta)*sind(theta(i))/(cosd(alpha)*cosd(theta(i))*cosd(beta)+sind(alpha)*sind(beta))*h; % end

% for i=1:n

y(i)=(sind(alpha)*cosd(beta).*cosd(theta(i))-cosd(alpha)*sind(beta))/(cosd(alpha)*cosd(theta(i))*cosd(beta)+sind(alpha)*sind(beta))*h; end

[theta,rough]=ji(x,y); f=fun(rough,theta); if f

jingdu1=[jingdu1,jingdu]; weidu1=[weidu1,weidu]; end

function l=xianxin(jingdu,x,y) p=polyfit(x,y,1); t0=14.7;

T0=t0-(120-jingdu)/15;

x0=x(1)-(T0-12)*(x(21)-x(1)); y0=x0*p(1)+p(2);

l=(x0^2+y0^2)^(1/2); end

%主程序：本程序采用遗传算法接力进化，

%将上次进化结束后得到的最终种群作为下次输入的初始种群 clc;

close all; clear all; %进化的代数 T=100;

optionsOrigin=gaoptimset('Generations',T/2); x=zeros(1,3);

[x,fval]=ga(@ch14_1f,3);

%,optionsOrigin||reason,output,finnal_pop %进行第二次接力进化

% options1=gaoptimset('Generations',T/2,'InitialPopulation',finnal_pop,... % 'PlotFcns',@gaplotbestf);

[x,fval]=ga(@ch14_1f,3,[],[],[],[],[0;-90;0],[180;90;365]); Bestx=x

BestFval=fval

function y=distance(address,order) nmb=size(address,1); y=0;

for i=1:nmb-1

y=y+sqrt((address(order(i+1),1)-address(order(i),1))^2+(address(order(i+1),2)-address(order(i),2))^2); end

y=y+sqrt((address(order(i+1),1)-address(order(1),1))^2+(address(order(i+1),2)-address(order(1),2))^2);

function y=exhgpath(order) while 1

b=size(order,1); r=unidrnd(b,1,2); if r(1)-r(2)~=0 break end end

b=order(r(2));

order(r(2))=order(r(1));

order(r(1))=b; y=order;

function [theta1,rough]=ch(jingdu1,weidu1,day) % load Qb.mat % x=Qb(:,2); % y=Qb(:,3);

[x,y]=yizhi(jingdu1,weidu1,day); g=length(weidu1); t1=12.7:1/20:13.7; n=length(t1);

beta0=chiwei(2015,day); beta=beta0; %beta=18; for i=1:g

b(i)=xianxin(jingdu1,x,y);

h(i)=abs(tand(90-(beta+weidu1(i)))*b(i)); end

for i=1:g

t=dangdishi(jingdu1,120,t1); % b=xianxin(jingdu1(i),x,y); % h=tand(90-(beta+weidu1(i)))*b; % h=fun2(jingdu1(i),x,y); alpha=weidu1(i); for j=1:n

theta(j)=(12-t(j))*15;

x(j)=-cosd(beta)*sind(theta(j))/(cosd(alpha)*cosd(theta(j))*cosd(beta)+sind(alpha)*sind(beta))*h(i);

y(j)=(sind(alpha)*cosd(beta).*cosd(theta(j))-cosd(alpha)*sind(beta))/(cosd(alpha)*cosd(theta(j))*cosd(beta)+sind(alpha)*sind(beta))*h(i); end

[theta1(i,:),rough(i,:)]=ji(x,y); end

%子函数：适应度函数同时也是目标函数, 函数存储名称为ch14_2f.m function f=ch14_1f(x) jingdu=x(1); weidu=x(2); day=x(3);

f=ff(jingdu,weidu,day); end

问题三

clc,clear; N=41;

for i=1:41

b=huidu(i);

data(:,:,i)=rgb2gray(b); end

for i=1:41

a(:,:,i)=data(200:1000,600:1800,i); end

% a(:,:,2)=a(:,:,2)-a(:,:,1); %

for i=1:41

for j=1:801

for k=1:1201

c(j,k,i)=a(j,k,i)-a(j,k,1); end end end

imshow(c(:,:,2))

load c.mat for i=1:41

d(:,:,i)=c(600:700,1000:1100,i); end

for i=2:41 if i

[j,k]=find(d(:,:,i)>=10); [A,B]=max(k); x(i-1)=A+1000; y(i-1)=k(B)+600; else

[j,k]=find(d(:,:,i)>=15); [A,B]=max(k); x(i-1)=A+1000; y(i-1)=k(B)+600; end end

x=x-277; y=y-678;

k1=2/((878-320)^2+9^2)^(1/2); x=k1*x; y=k1*y;

p4=k1*(960^2+540^2)^(1/2);

clc,clear; e11=[-62,0,1]; e12=[43,0,1]; e21=[43,0,1]; e22=[43,21,1]; e31=[43,0,1]; e32=[59,4,1]; a=cross(e11,e12); b=cross(e21,e22); c=cross(e31,e32);

m1=[a(1)*a(1),a(1)*a(2),a(1)*a(3); a(2)*a(1),a(2)*a(2),a(2)*a(3); a(3)*a(1),a(3)*a(2),a(3)*a(3)]; m2=[b(1)*b(1),b(1)*b(2),b(1)*b(3); b(2)*b(1),b(2)*b(2),b(2)*b(3); b(3)*b(1),b(3)*b(2),b(3)*b(3)]; m3=[c(1)*c(1),c(1)*c(2),c(1)*c(3); c(2)*c(1),c(2)*c(2),c(2)*c(3); c(3)*c(1),c(3)*c(2),c(3)*c(3)]; m=m1+m2+m3; [x,y]=eig(m);

eigenvalue=diag(y); y_lamda = x(:, 3);

问题二的实测照片：

问题三的实测xy 坐标：