太阳影子定位技术-2015高教社杯-数学建模-获奖论文

1太阳影子定位技术摘要本文以太阳影子定位技术为背景，结合直杆影子轨迹的变化规律建立数学模型。并运用视频数据分析的方法，确定拍摄地点及日期等地理信息条件。第一问给出了北京时间、拍摄日期，以及拍摄地点的经纬度。我们可以结合太阳赤纬、时角、直杆的经纬度与太阳高度角之间的关系建立模型，求出符合时间条件要求的太阳高度角，再根据已知的杆的高度和三角公式求出影长关于时间的变化曲线。第二、三问在第一问的基础上增加难度，使部分变量未知。通过文献查阅和方程推导，得出阴影运动轨迹形状是双曲线的一支，并且具体形状和当地的纬度以及赤纬有关，本文根据这点进行模型假设与建立。附件中给出的坐标并不一定是标准地理坐标，通过对其进行坐标变换，引入了实际坐标系与标准地理坐标系的偏角。在拟合多项高次变量组成的隐函数方程的过程中，为增加精确度，运用最小二乘法进行拟合求解未知参量时，可以利用直杆阴影顶点轨迹的形状，建立参量和变量之间的关系，简化需拟合的隐函数方程。这样就可以根据太阳影子顶点横纵坐标以及对应的时刻，把偏角、纬度、经度、日期作为未知参数进行拟合，得出要求的地理位置和相应的日期。如通过对附件1数据的拟合求解可得到一组地理坐标（东经104.425度，北纬15.6578度），对附件2数据的拟合求解可得一个可能的日期6月21日，坐标（东经116度，北纬26度），由附件3得到的可能的日期地点为：6月21日，（东经164.55度，北纬71.26度）。为了便于定位，根据一般工程的实际需求，对美国天文学家纽康（NewComb）提出的太阳公式作了综合、简化，舍去了一些高阶微小量。结合测量学的理论，用数学模型进行非线性拟合求得直杆所处的经纬度。第四问给出一段视频，实际是对前三问模型的实际应用。本问对一些已有的论文以及专利进行借鉴，创新与简化。首先对视频中的图像进行取帧，在灰度处理中因为技术限制，改为运用Matlab二值化处理。并根据简单测量画出运行轨迹。运用主元分析法求得阴影尖端坐标与杆底坐标的关系。确定影子的运动轨迹。之后借鉴已有成熟理论将2D图像去畸变，恢复仿射的度量属性，通过对3D图形转变2D过程的逆向推导，将坐标恢复为符合现实要求的坐标。之后回归前几问建立的的日晷数学模型进行求解，得到一个可能的地理坐标为（东经104.9度，北纬25.33度）。并在最后进行误差修正。关键词：日晷投影原理、杆影端点轨迹、非线性最小二乘法、主元分析法、二值化处理、Floodfill图论算法2一、问题重述如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用建立的模型画出2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线。2.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。将模型应用于附件1的影子顶点坐标数据，给出若干个可能的地点。3.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干个可能的地点与日期。4．附件4为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度为2米。请建立确定视频拍摄地点的数学模型，并应用模型给出若干个可能的拍摄地点。如果拍摄日期未知，能否根据视频确定出拍摄地点与日期？二、问题分析2.1问题一的分析本问以直杆影子长度为研究对象，寻找影响影子长度与各个参数的关系及其变化规律。为了使影长的计算科学严谨，我们应了解太阳与地球之间相互的运动轨迹，由此计算出太阳对地球上某一定点的相对位置。这主要由当地的地理纬度、季节(月、日)和时间三个因素决定，可以用地理纬度(φ)、太阳赤纬角(δ)、太阳高度角(h)、及时角t等参数进行定量表达。2.2问题二的分析由于附件中给出的直杆阴影顶点的坐标系的x轴和y轴并不一定垂直或重合，有可能坐标轴向与南北方向存在一定的偏角（θ）。本问根据太阳影子顶点横纵坐标的21组数据，对x和y的坐标进行旋转坐标变换，得到阴影顶点在新坐标系下的坐标，该新坐标系以正东方向为x’轴正向，正北方向为y’轴正向，直杆底端为坐标原点，由原坐标系旋转θ得到。这样就可以由x’和y’求出相应时刻的太阳方位角（A），再结合太阳高度角(h)的计算公式，经过一系列化简，可以得到经纬度之间的关系。将经纬度作为参数，利用Matlab进行非线性拟合，选取合适的初值，即可得直杆所处地点的经纬度。2.3问题三的分析该问题的求解可利用问题二建立起来的模型，将由日期确定的太阳赤纬作为未知参数，在Matlab中对时间和直杆影子长度进行非线性拟合，选取合适的初3值得到经纬坐标和日期的值。2.4问题四的分析本问考察基于视频数据分析方法进行太阳影子定位。从图像或视频中估算经纬度是目前计算机视觉领域的研究热点问题，估算经纬度不仅自身具有重要理论意义，而且它对计算机视觉问题也有积极的启示意义。第四问提供的视频，其中体现了标志物的影子在一段时间内的移动轨迹。在这种条件下求经纬度，实际上就是基于视频中太阳影子轨迹来估计经纬度的实际应用。首先我们把视频取帧处理得到影子的轨迹点，在拟合出地平线后，运用计算机作图的相关知识对坐标进行纠正后，把图片上的2D坐标恢复为实际中真实的3D坐标，最后把问题回归日晷模型算出经纬度，并对因为地方时标准时引起的误差进行修正。三、模型假设1.将太阳光近似地看成平行光投射到地球；2.忽略太阳光受地球大气层折射和漫反射的影响；3.地球和太阳在运行中不规则变化和周期性变化产生的误差忽略不计；4.假设直杆所在的地面为水平的；5.忽略地球形状对试验结果的影响；6.假设每天的时间为24小时整；7.假设所求日期均均为2015年日期。四、符号说明符号符号说明φ物体的地理纬度ω物体的地理经度δ太阳对应日期的赤纬角h太阳高度角A太阳方位角T时角H直杆长度五、模型建立及求解5.1问题一模型的建立及求解5.1.1模型的建立4影响影长的因素：hcotHL其中h为太阳光线与水平地面的夹角，即太阳高度角；通过查阅目前国内的大部分天文学文献，可得太阳高度角的计算公式[1]，即：Thcosδcosφcosδsinφsinsin其中T为时角，可用1215tT计算，t为24小时制的当地时间;δ为赤纬角，其较精确公式为[2]：θ0.0201cos3θ0.3656cos20.758cosθ-θ0.1712sin3-θ0.1149sin2θsin2567.233723.0δ式中θ称日角，即2422.365π2θn这里n又分两部分组成，即n=N-N0；式中N为积日，即日期在年内的顺序号；N0的计算公式如下：41985-19852422.06764.790年份年份INTN其中INT(X)为取不大于X的最大整数。5.1.2模型的求解根据问题一的题目已知，10月22日是2015年的第295天，即N=295。对于日照计算来说，地球和太阳在运行中不规则变化和周期性变化产生误差的数值相当小，可以忽略不计。所以此赤纬角δ计算公式的结果较为精确，可以满足计算影长变化曲线精度要求。由此可建立影长和北京时间的数学函数，并作出图形。Matlab程序见附录一，影子长度变化曲线见图1。(3)(1)(2)(4)(5)5图1直杆的太阳影子长度的变化曲线5.2问题二模型的建立及求解5.2.1模型的建立（1）坐标变换第二问的附件虽然在直杆所在的地平面处建立了正交分解的x-y轴，但并没有指明坐标轴和地理正北方向的夹角。据此分析，可设θ为地平面坐标系的x轴正向与正东方向的夹角。对x和y进行下式给出的变换，可得直杆顶点在新坐标系下的投影坐标。变换如下：sincos'ysincos'xxyyx（2）杆影端点移动轨迹关于杆影端点移动轨迹的周年变化规律如图2所示，该图描绘了一年内不同日期的轨迹。夏至日和冬至日是太阳直射点移动方向发生转换的日期，春分日和秋分日是直射点的南北半球位置转换的日期。因此，选择的这几个日期具有较强的代表性。在图中，AA′、CC′、EE′依次表示冬至日、两分日、夏至日的轨迹，BB′、DD′表示介于二分二至日的轨迹（如立冬和立夏），该图清晰的反(6)6映除两分日之外，其余日期的轨迹图都是双曲线中的一条。图2杆影端点移动轨迹的周年变化规律而对于具体的某一天中直杆阴影顶点的运动轨迹的描绘，则可以借助Analemmatic日晷的模型来描绘。在这种日晷的模型中，日晷采用了垂直的指时针，它的时间线是赤道日晷的时间线在地平面的投影水平方向指示东西轴，垂直方向指向南北轴，投影日晷的位置沿着短轴移动，如图3所示。因此具有时间一致对应关系的阴影轨迹如果投影到水平平面上将得到一个二次曲线，二次曲线的对称轴位于南北方向指示轴上。轨迹的形状与当地的纬度以及太阳直射点的纬度即太阳赤纬有关，而太阳赤纬又与日期相关，故轨迹的形状与当地的纬度和日期相关联。图3日晷模型及影子变化规律变化规律7（3）太阳方位角和太阳高度角有了（1）得到的新坐标，就可以用''yx来表示太阳方位角的正切值。由天文学的相关知识[3]，可以得到以下等式：φcoscosδsin-φsinsincoshhAA为太阳方位角，并且''tanyxA又因太阳高度角公式为Thcosδcosφcosδsinφsinsin以及°180°1512°15°120tT联立化简可得到如下等式：btaTxyyxcossincossincossincos其中cossinacossinb5.2.2模型的求解根据表格附件一做出杆影的轨迹图如下(7)(8)(9)(10)(11)(12)(13)8图4附件一杆影顶点轨迹图对上述结果进一步化简可得：TTYYcsccostancotsintantan1其中yxY上式建立了阴影顶点坐标与时间（t）、经度（ω）、纬度（φ）以及所用坐标系与地理坐标系之间偏角（θ）之间的关系模型。由于附件已给出直杆阴影顶点在每一时刻的实际坐标，可以将其它待求的量作为未知参量，基于非线性最小二乘法在Matlab中对Y和t进行拟合，结合地理资料选取合适的初值求得直杆可能的地点。我们得到了（104.4249°E，15.65784°N）、（33.0454°E，16.7577°N）等地点。5.3问题三模型的建立与求解5.3.1模型的建立基于问题二的模型，我们很自然地得出问题三的模型：(14)(15)9TtYYcsccostancotsintantan1其中δ为太阳赤纬，是与日期有关的量，问题三相对于问题二的变化是日期作为参量是未知数。结合问题一中赤纬的计算公式，可用日期N表示出δ。这样就得到了含纬度φ、经度ω、日期N、偏角θ四个未知参量的隐函数方程，也就是问题三的模型。5.3.2模型的求解类似于问题二的求解，先画出附件二、三中杆影端点轨迹的变化规律如下图：图5附件二杆影端点轨迹图图6附件三杆影端点轨迹图(16)10关于时间的已知参量在第三问中变为未知参量，因此可以在Matlab中调用非线性拟合函数拟合隐函数方程来求解纬度φ、经度ω、日期N、偏角θ四个未知参量。对于附件二日期和地理坐标的求解，通过选取合适的不同初值，我们得到了以下不同的地点和日期，如下所示：（1）6月1日，（21°N,101°E）（2）6月21日，（30°N,122°E）（3）6月21日，（26°N,116°E）其中括号内为该地点的经纬度。对于附件三，我们也求解出了一系列点，如下所示：（1）6月21日，（71.26°N，164.55°E）（2）6月21日，（52.85°N，140.11°E）5.4问题四模型的建立与求解5.4.1模型的建立问题四实质与问题三的原理相同