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

太阳影子定位技术摘要本文以太阳影子定位技术为背景，结合直杆影子轨迹的变化规律建立数学模型。并运用视频数据分析的方法，确定拍摄地点及日期等地理信息条件。第一问给出了北京时间、拍摄日期，以及拍摄地点的经纬度。我们可以结合太阳赤纬、时角、直杆的经纬度与太阳高度角之间的关系建立模型，求出符合时间条件要求的太阳高度角，再根据已知的杆的高度和三角公式求出影长关于时间的变化曲线。第二、三问在第一问的基础上增加难度，使部分变量未知。通过文献查阅和方程推导，得出阴影运动轨迹形状是双曲线的一支，并且具体形状和当地的纬度以及赤纬有关，本文根据这点进行模型假设与建立。附件中给出的坐标并不一定是标准地理坐标，通过对其进行坐标变换，引入了实际坐标系与标准地理坐标系的偏角。在拟合多项高次变量组成的隐函数方程的过程中，为增加精确度，运用最小二乘法进行拟合求解未知参量时，可以利用直杆阴影顶点轨迹的形状，建立参量和变量之间的关系，简化需拟合的隐函数方程。这样就可以根据太阳影子顶点横纵坐标以及对应的时刻，把偏角、纬度、经度、日期作为未知参数进行拟合，得出要求的地理位置和相应的日期。如通过对附件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图论算法

一、问题重述如何确定视频的拍摄地点和拍摄日期是视频数据分析的重要方面，太阳影子定位技术就是通过分析视频中物体的太阳影子变化，确定视频拍摄的地点和日期的一种方法。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中对时间和直杆影子长度进行非线性拟合，选取合适的初映除两分日之外，其余日期的轨迹图都是双曲线中的一条。图2杆影端点移动轨迹的周年变化规律而对于具体的某一天中直杆阴影顶点的运动轨迹的描绘，则可以借助Analemmatic日晷的模型来描绘。在这种日晷的模型中，日晷采用了垂直的指时针，它的时间线是赤道日晷的时间线在地平面的投影水平方向指示东西轴，垂直方向指向南北轴，投影日晷的位置沿着短轴移动，如图3所示。因此具有时间一致对应关系的阴影轨迹如果投影到水平平面上将得到一个二次曲线，二次曲线的对称轴位于南北方向指示轴上。轨迹的形状与当地的纬度以及太阳直射点的纬度即太阳赤纬有关，而太阳赤纬又与日期相关，故轨迹的形状与当地的纬度和日期相关联。图3日晷模型及影子变化规律变化规律（3）太阳方位角和太阳高度角有了（1）得到的新坐标，就可以用来表示太阳方位角的正切值。由天文学的相关知识[3]，可以得到以下等式：(7)(7)A为太阳方位角，并且(8)(8)又因太阳高度角公式为(9)(9)以及(10)(10)联立化简可得到如下等式：(11)(11)其中(12)(12)(13)(13)5.2.2模型的求解根据表格附件一做出杆影的轨迹图如下图4附件一杆影顶点轨迹图对上述结果进一步化简可得：(14)(14)其中(15)(15)上式建立了阴影顶点坐标与时间（t）、经度（ω）、纬度（φ）以及所用坐标系与地理坐标系之间偏角（θ）之间的关系模型。由于附件已给出直杆阴影顶点在每一时刻的实际坐标，可以将其它待求的量作为未知参量，基于非线性最小二乘法在Matlab中对Y和t进行拟合，结合地理资料选取合适的初值求得直杆可能的地点。我们得到了（104.4249E，15.65784N）、（33.0454E，16.7577N）等地点。5.3问题三模型的建立与求解5.3.1模型的建立基于问题二的模型，我们很自然地得出问题三的模型：(16)(16)其中δ为太阳赤纬，是与日期有关的量，问题三相对于问题二的变化是日期作为参量是未知数。结合问题一中赤纬的计算公式，可用日期N表示出δ。这样就得到了含纬度φ、经度ω、日期N、偏角θ四个未知参量的隐函数方程，也就是问题三的模型。5.3.2模型的求解类似于问题二的求解，先画出附件二、三中杆影端点轨迹的变化规律如下图：图5附件二杆影端点轨迹图图6附件三杆影端点轨迹图关于时间的已知参量在第三问中变为未知参量，因此可以在Matlab中调用非线性拟合函数拟合隐函数方程来求解纬度φ、经度ω、日期N、偏角θ四个未知参量。对于附件二日期和地理坐标的求解，通过选取合适的不同初值，我们得到了以下不同的地点和日期，如下所示：（1）6月1日，（21N,101E）（2）6月21日，（30N,122E）（3）6月21日，（26N,116E）其中括号内为该地点的经纬度。对于附件三，我们也求解出了一系列点，如下所示：（1）6月21日，（71.26N，164.55E）（2）6月21日，（52.85N，140.11E）5.4问题四模型的建立与求解5.4.1模型的建立问题四实质与问题三的原理相同，都是给定时刻和该时刻下现实坐标系中的横纵坐标，求解日期与地理坐标；不同点在于本题的坐标并不是直接给出的，需要借助视频处理技术从中提取出杆影端点的坐标。因此，基于对问题二、三的分析可以建立问题四的模型，即：(17)(17)模型与问题三相同。5.4.2模型的求解（1）从视频中提取杆影端点坐标观看时长为40分钟的视频录像，根据直杆的高度为2m，并结合图中杆长和影长的比例关系，即可求得实际的杆影的长度随时间的变化情况。但众所周知，广角镜头所产生的图片或视频都存在较为明显的透视畸变，即被摄体离镜头越远，在屏面的成像越小。为了避免透视畸变对实际影长测量带来的影响，引入主元分析法解决杆影尖点的确定问题。本题模型引用一种半自动输入阴影轨迹检测的方法，首先以1分钟为单元对所给视频做取帧处理，形成一组图片并构建背景图片B：(18)(18)在B中，每一个像素点（x,y）是V中最亮的像素点，且是基于灰度级的。设置背景绝对差值，阴影点就可以显示出来。然后在阴影区域运用floodfill图论算法进行阴影点计算，最后的突出阴影点用主成分技术（PCA），我们同时运用MATLAB对图像进行二值化处理获得效果如图8所示。图8二值化图求相对坐标以直杆底端为原点建立正交的x-y坐标系，单位为像素点。结合主元分析法（PCA）合理得到准确的阴影运动轨迹。（2）相机模型，消隐点的确定以及对图像的纠正处理世界坐标中阴影点的位置取决于投影物所在纬度位置以及太阳光朝向与投影平面之间的几何关系。第四问进行的纬度估计技术是基于计算机视图的研究算法，采用几何分析算法来进行视频地理位置的估算。由于存在透视畸变，照片中的景物会和实际景物有很大的差别。我们生活中有许多透视变化的例子：照片中的景物轮廓形状改变，但透视变化下只有直线被保留。所以，本问的数学模型我们在欧式几何中加入一些无穷远的理想点形成透视空间，定义一个平面透视几何变换作为任何平面中点以保持直线性的映射。在进行地理经纬度定位中，视频和相机校准是十分重要的一步，要进行平面但因性变换。在透视相机的作用下，有些几何属性是保持的，例如共线性，即一条直线透视变换后仍为一条直线，然而一般的透视情况下平行直线将6不再保持平行。透视几何模型决定了透视效果同时也提供了相应计算的数学表达式。在计算机视觉中，视觉几何是用来研究在各种变化下仍然保持不变的属性。从这个角度来分析，2D的透视几何就是来研究在2D平面下经过各种变换之后仍然保持不变的属性特征。这种2D平面变换对于点或者直线来说是不可逆的。另外，透视变换的逆变换也是一种透视变换，所以这种变换同时存在两种变换：共平面变换和透视变换，也叫做平面单应性。为了理解这种变换，可以假设用一个三维向量来表示一个二维空间的点x，那么Hx就是这个点经过平面线性变换后得到的齐次坐标。这种变换法则认为任何透视效果都会引起齐次坐标的线性变换，并且逆变换也是这种线性透视。一个平面透视变换是一个采用三维向量的线性变换，可以用一个3×3的矩阵表示：(19)(19)一个透视变换能把每一个照片投影成为一个透视等价的照片，并且保留所有的不变属性。由此可见，经过相机中心的直线定义了一个平面与另一个平面的变换关系。实际上，如果两个坐标系统都是欧式空间系统，那么这种由中心透视变换决定的映射就会比任意的透射变换要产生更多的限制，而这就叫做透视变换而不是完全的投影变换。2D平面的等度量变换能够保留欧氏距离不变，一个等度量变换可以表示成：(20)(20)相比之下，给定一个投影变换面积的等级大小变换随着位置的不同而发生变化。例如在透视变化下,相同平面中一个较远的正方形会比较近的正方形具有更小的透视结果。并且经过投影变换的直线的朝向不仅依赖于原始直线的朝向而且依赖于它的位置。在本文实现平面单应性计算的过程中,需要一个重要的步骤：对投影变换进行分解。一个投影变化可以被分解成一系列的变化,其中每一个矩阵都比前一个具有更高级的变换。(20)(20)其中A是一个非奇异矩阵，K是一个上三角矩阵并且它的行列式为1。这种分解只有在v不等于的时候才成立，并且当s被确定才会有唯一的分解结果。分解出来的三个矩阵代表了相应类型变换的实质变化。考虑到从透视照片中进行校准，HP将直线恢复到无穷远；HA影响了仿射属性，但是没有将直线恢复到无穷远；最后HS代表一般相似变换，并且不会影响仿射或者透视属性。从照片中进行透视纠正的目的是消除透视照片中的投影畸变，以便得到相似属性，即角度和长度比值等可以直接从照片平面中测量。投影畸变可以通过照片平面中的四组对应点来消除，并且能够明确的计算出参考点与相应的像点之间的映射。而这实际上过度具体化了几何关系，这是因为一个投影变换与相似变换比较起来只有四个自由度，所以仅仅需要指定四个自由度而非八个就可以确定度量属性了。在投影几何中，这四个自由度给出了几何物体的物理形状：消隐线提供了两个自由度，两个虚圆点也提供了两个自由度。因此，当消隐点直线的像被确定了，那么透视畸变就可以被消除了，同时如果虚圆点被确定了，仿射畸变也会被消除。因此，唯一的畸变就是相似变换了。在投影变换下，一个理想点被映射到一个有限的点，因此消隐线被映射到一条有限的直线。但是，如果这种变换是仿射的，那么消隐线不会被映射为有限的直线而是仍然保持在无穷远的位置。这种变换的逆变换也是成立的，即如果仿射变换是保持无穷远直线的最一般的线性变换。但是消隐线上的点对于仿射变换不是点对点的保持。也就是说消隐点经过仿射变换后得到的点虽然位于消隐线上，但是却不是原来的点。一旦照片平面中无穷远处的直线被确定了，那么就有可能实现原始平面上的仿射测量。例如，原始平面上的平行直线可以被识别为平行，如果这些直线相交于无穷远直线的话。这是因为欧式平面中的平行线相交于无穷远的直线，在经过仿射变换后这些直线仍然相交于无穷远直线的像，因为透视变换下的相交是保持的。类似的，一旦无穷远的直线被确定了，一条直线上的长度比值就可以通过三个点确定的交比来确定。但是，稍微弯曲的直线却是更适合这种计算上的算法，这是因为仅仅需要简单地将消隐线变换到它的原始位置。实现了这种变换的投影矩阵可以应用到任何点中来实现仿射纠正，也就是经过了这种变换，仿射测量可以直接从校正后的照片中获得。根据以上分析编写程序，对图像处理得到视频对应的三维空间内的数据，并用Matlab拟合，结果如图10。图10视频中杆影端点轨迹变化的俯视图这样，我们就从视频中提取出了一组点如下，由于篇幅原因，在此只罗列其中的某几个点，详细数据在附录六中。时间/h8.928.938.958.978.989.019.029.03x值804795786780780772766760y值-15-15-15-15-13-13-12-12（3）编程求解求解过程与问题三的求解过程类似，赋初值后得到一个可能的拍摄地点为（25.33N，104.90E）详细的程序代码见附录六。5.5问题五此问题与问题三相同，可将日期作为未知参量，通过拟合得到日期，详细过程类似于问题三。误差分析及敏感度分析将通过数学模型求解的数值与使用经纬仪和有水准器的支架所实际观测到的直杆影长相比较。本模型方法的误差约为，对50m内的物高的测量精度可达厘米级要求。由于拟合的函数过于复杂，参量与变量之间并不能用显性的关系表示出来，这就给对模型进行的灵敏性分析造成困难，因此，我们通过定性的分析来说明某一参量的变化对其它量造成的影响。在范围和步长给定的情况下，下图为遍历所有可能的初始点得到的拟合出的参数的分布情况。从图中可以看出，拟合得到的点大概率地分布于某一特定区域，这说明建立的模型较为合理。图11不同初值下的所求参数的分布七、模型的评价及改进6.1视频数据分析的优点本文提供的太阳影子定位的数学模型，对于光照研究以及航空摄影、精密水准测量的最佳时间段的选择有一定的意义。此数学模型还可以推广到生活中高层建筑群的合理布局和农林间种的最佳距离及林带走向的优化设计中。第四问建立的视频经纬度分析模型也有广泛的应用。即使粗略的经纬度估计也能够提供有用的线索，预测当地气温、平均降雨量等大量背景信息。6.2模型求解的不足由于求解参数采用的是曲线拟合的方法，其对初值的依赖性较大，初值设置的不合理就有可能造成结果误差较大，导致精确性较差。八、参考文献[1]陈晓勇,郑科科.对建筑日照计算中太阳赤纬角公式的探讨,浙江建筑,第28卷,2011.[2]《建筑设计资料集》编委会.建筑设计资料集[M].2版.北京:中国建筑工业出版社,1994:179-185.[3]林根石,利用太阳视坐标的计算进行物高测量与定位,南京林业大学学报,第15卷,1991.[4]吴济廉,影端轨迹周年变化的实践与分析—以北温带地区为例.地理教学,2013年第10期.[5]吴济廉.关于杆影端点移动轨迹的讨论与求证.中学地理,2011年第3期.[6]天津大学.基于视频中太阳影子轨迹的经纬度估计方法:中国,200910067817.2012-04-11.九、附录附录一：Year=2015;N=295;%N为积日t0=9:0.1:15;T=(t0-12-0.24)*15*pi/180;%T为时角fai=39.9072*pi/180;%fai为物体的地理纬度H=3;N0=79.6764+0.2422*(Year-1985)-fix((Year-1985)/4);t=N-N0;c=2*pi*t/365.2422;%c为日角delta=0.3723+23.2567*sin(c)+0.1149*sin(2*c)-0.1712*sin(3*c)-0.758*cos(c)+0.3656*cos(2*c)+0.0201*cos(3*c);delta0=delta*pi/180;sinh=sin(fai)*sin(delta0)+cos(fai)*cos(delta0)*cos(T);L=H./tan(asin(sinh));plot(t0,L);title('影长变化曲线')xlabel('北京时间t')ylabel('影长')附录二：%本函数求解赤纬角functiondelta=chiwei(N,Year)N0=79.6764+0.2422*(Year-1985)-fix((Year-1985)/4);%求解积日a=N-N0;c=2*pi*a/365.2422;delta=0.3723+23.2567*sin(c)+0.1149*sin(2*c)-0.1712*sin(3*c)-0.758*cos(c)+0.3656*cos(2*c)+0.0201*cos(3*c);%求解赤纬角end附录三：附件一求解程序：clear;%clca=10.6306*pi/180;F=@(p,x)(1-x(:,2).*tan(p(1)))./(x(:,2)+tan(p(1)))-sin(p(2)*pi/180).*cot((p(3)/15-x(:,1)-20)*15*pi/180)+tan(chiwei(108,2015)*pi/180)*cos(p(2)*pi/180)./sin((p(3)/15+x(:,1)-20)*15*pi/180);x=[14.70.59093101814.750.57443063214.80.55867608914.850.5435737514.90.52917945214.950.515360856150.50208868915.050.48933541415.10.47703894715.150.46524911615.20.45383972315.250.4428299315.30.4321740915.350.42188874215.40.41193212215.450.40226832515.50.39289045715.550.38379141515.60.3749206915.650.36631837215.70.357917966];p0=[836115]%赋初值warningoffp=nlinfit(x,zeros(size(x,1),1),F,p0);disp(num2str(p));plot(x(:,1),x(:,2),'ro');holdon;ezplot(@(x,y)F(p,[x,y]),[0,1,-1e-3,1e-3]);title('拟合曲线');legend('样本点')附录四：附件二求解程序：clear;clcF=@(p,x)(1-x(:,2).*tan(p(1)))./(x(:,2)+tan(p(1)))-sin(p(2)*pi/180).*cot((p(3)/15+x(:,1)-20)*15*pi/180)+tan(chiwei(p(4),2015)*pi/180)*cos(p(2)*pi/180)./sin((p(3)/15+x(:,1)-20)*15*pi/180);x=[12.68-773-6.39206349212.78-5.76806640612.83-5.24103495212.88-4.7882003412.93-4.39840319413.98-4.05427817613.03-3.75125089313.08-3.48197278913.13-3.2402597413.18-3.02175264113.23-2.82349433513.28-2.64248853213.33-2.47720364713.38-2.32497331913.43-2.1836945313.48-2.05323669113.53-1.9318622713.58-1.81794630213.63-1.71163536413.68-1.611730598];p0=[540116250]%赋初值warningoffp=nlinfit(x,zeros(size(x,1),1),F,p0)plot(p(3),p(2),'.');holdondisp(num2str(p));plot(x(:,1),x(:,2),'ro');holdon;ezplot(@(x,y)F(p,[x,y]),[0,1,-1e-3,1e-3]);title('拟合曲线');legend('样本点')附录五：附件三求解程序：clcF=@(p,x)(1-x(:,2).*tan(p(1)))./(x(:,2)+tan(p(1)))-sin(p(2)*pi/180).*cot((p(3)/15+x(:,1)-20)*15*pi/180)+tan(chiwei(p(4),2015)*pi/180)*cos(p(2)*pi/180)./sin((p(3)/15+x(:,1)-20)*15*pi/180);x=[13.150.34883093513.20.3667377413.250.38478430913.30.40294684813.350.42128134713.40.43975703413.450.45836359213.50.47717150913.550.49610263613.60.51522603813.650.53450633613.70.55396733713.750.57360236213.80.59342265513.850.61345357113.90.63365129713.950.654077122140.67470503914.050.69555806214.10.71661992414.150.737931245];p0=[526115222];%赋初值warningoffp=nlinfit(x,zeros(size(x,1),1),F,p0);disp(num2str(p));plot(x(:,1),x(:,2),'ro');holdon;ezplot(@(x,y)F(p,[x,y]),[0,1,-1e-3,1e-3]);title('拟合曲线');legend('样本点')附录六：问题四求解程序：clear;clcF=@(p,x)(1-x(:,2).*tan(p(1)))./(x(:,2)+tan(p(1)))-sin(p(2)*pi/180).*cot((p(3)/15+x(:,1)-20)*15*pi/180)+tan(chiwei(p(4),2015)*pi/180)*cos(p(2)*pi/180)./sin((p(3)/15+x(:,1)-20)*15*pi/180);x=[8.98.9166666678.9333333338.958.9666666678.98333333399.0166666679.0333333339.059.0666666679.0833333339.19.1166666679.1333333339.159.1666666679.29.2166666679.2333333339.2833339.3333339.429.59.57-1.552118549-1.552141775-1.551930641-1.551714674-1.554131203-1.554131203-1.555253537-1.555131812-1.555008165-1.55760448-1.557499238-1.558764822-1.558700143-1.558634764-1.559897576-1.561154752-1.562509224-1.563803434-1.562321953-1.551024892-1.566486009-1.567859478