赵珮瑶

摘要：该文针对太阳影子定位问题建立了依赖于拍摄地点经度、纬度、日期的影子长度变化的数学模型，分析确定拍摄地点的日期。

关键词：超定方程组；最小二乘；数据矫正

中图分类号：TP311 文献标识码：A 文章编号：1009-3044（2016）17-0260-03

1 问题重述

1.1 问题背景

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

1.2 目标任务

问题一：建立影子长度变化的数学模型，分析影子长度的变化规律，并画出2015年10月22日北京时间9：00-15：00之间天安门广场（北纬39度54分26秒，东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线。

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

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

如果日期已知，该如何处理？

2 模型的建立与求解

2.1 问题一模型建立与求解

2.1.1 问题一模型的建立

我们采用太阳高度角和三角函数相结合的思路。直杆与太阳照射下成的影子可构成直角三角形，其三角函数关系式可以表示为

[sinθ=HH2+L2] ⑴

其中，[H]表示直杆长度，[L]表示直杆影长，[sinθ]表示太阳高度角的正弦值。

由参考文献[1]，太阳高度角随着地方时和太阳赤纬的变化而变化。太阳高度角的计算公式为：

[sinθ=sinφsinδ+cosφcosδcost] ⑵

上式中，[θ]，。表示太阳高度角，[φ]表示观测地理纬度（太阳赤纬和地理纬度都是北纬为正，南纬为负），[δ]表示太阳赤纬（与太阳直射点维度相等），[t]表示时角。

2.1.2 问题一模型的求解

（1）太阳赤纬[δ]与地理纬度[φ]正弦值的确定

根据太阳赤纬的公式：

[sinδ=0.39795cos0.98563N-173] ⑶

式中[N]表示观测日期，自每年1月1日开始计算。即代入2015年10月22日，[N]=295。

计算可得：

[sinδ=0.39795cos0.98563295-173=-0.20046] ⑷

因为天安门广场地处（北纬39度54分26秒，东经116度23分29秒），与代表北京时间的东八区经度差为0，所以天不考虑安门广场的经度位置，得

[sinφ=sin[（39+54/60+26/3600）]=0.641546] ⑸

（2）太阳高度角的正弦值的求解

由公式⑷、⑸计算太阳赤纬[σ]与地理纬度[φ]的余弦值如下：

[cosδ=0.979702cosφ=0.767084] ⑹

由参考文献[3]可知，时角公式为：

[t=（T-12）×15]° ⑺

我们在计算时角时发现，北京时间并不是北京（东经116.4°）地方的时间，而是东经120°地方的地方时间。根据经度变化15°，时间变化1小时这一地理规则，我们推算出北京当地的真正时间比代表东八区的北京时间晚约14分钟。计算[sinθ]如下表1所示：

（3）计算直杆的太阳影子的长度变化

结合公式⑴的三角函数公式以及表1数据，计算不同时刻影长如下表2所示：

通过MATLAB拟合曲线，我们可知直杆的太阳影子长度[L]最短的时间并不是在12时，而是在11时至12时之间。

2.2 问题二模型建立与求解

2.2.1 问题二模型的建立

（1）直杆直角坐标系的建立

假设可能地区在太阳直射点的以北地区数据为正。2015年4月18日的太阳直射点在赤道和北回归线之间，且相对靠近赤道的位置。我们假设，在附表中所给14：42-15：42的时间段内，太阳在直杆的南面。并且在14：42-15；42这一时间段，已经过正午太阳高度角所处的12时，太阳从直杆的西面进行照射。综上可得到假设，直杆的影子在它的东北方向。以直杆底端为坐标原点，分别以东向、北向为[x]轴、[y]轴的方向，建立直角坐标系，如图3所示：

图1 直杆直角坐标系的建立

（2）经度定位模型的建立

由参考文献[5]，太阳方位角与观测地点经度、纬度以及观测时刻的关系如下：

[A=arcsin（cosδsintcosθ）] ⑻

根据太阳高度角的余弦值为

[cosθ=LH2+L2] ⑼

结合公式⑻和公式⑼，化简为

[LsinA=H2+L2cosδsint] ⑽

在上式⑽中，我们可根据附件一的数据求解太阳方位角的正弦值，由表3和公式⑶可知不同时间的影长和太阳赤纬的余弦值。最终我们得到了未知量为直杆高度[H]、时角[t]的经度的定位模型。

（3）纬度定位模型的建立

结合问题一中有关太阳高度角的公式⑴⑵⑺，得到下列方程：

[H=H2+L2C1sinφ+C2cosφcos[（T-12）×15]] ⑾

上式中，[C1]表示太阳赤纬的正弦值，[C2]表示太阳赤纬的余弦值，因为时间是确定的，所以太阳赤纬的正弦值、余弦值都是常数。从而我们得到了未知量为直杆高度[H]、地理纬度[φ]的纬度的定位模型。

2.2.2 问题二模型的求解

（1）太阳赤纬的正弦值、余弦值的确定

太阳方位角是指太阳光线在地平线上的投影与当地经线的夹角。根据所给数据，由勾股定理[L=x2+y2]及太阳方位角计算公式：

[A=arctanxy] ⑿

我们计算出不同时间的影子长度，见下表3。

（2）利用解超定方程组的最小二乘解思想求解经纬度的定位模型

对于经度的定位模型，我们所给数据的21个不同组的数据，得到21组不同长度的直杆影长，建立了含有21个不同的方程，只有直杆高度[H]、时角[t]两个未知量。利用[MATLAB]中的[fsolve]函数，假定直杆高度[H]的初值为0米，时间初值为上午9时。计算得直杆的高度[H1]为1.0134米，推测测量地点的时间为14时55分5秒。已知代表东八区的北京时间（东经120）为14时42分，对比14时55分5秒和14时42分，根据地区经度与时间差的地理知识，计算出该地的经度为东经1231612。

同样，对于纬度的定位模型，我们可以根据附件一的21个组不同的数据，建立了含有21个不同的方程，只有直杆高度[H]、地理纬度[φ]两个未知量。计算可得，直杆的高度[H3]为1.0440米，地理纬度[φ]的弧度角为-0.1831[rad]。根据纬度算法，得到以下的计算公式⒀：

当地理纬度[φ]的弧度角为0.3025[rad]时，将其数值代入上述公式⑿，算的第一组结果的纬度值为北纬17°18；当地理纬度[φ]的弧度角为-0.1831[rad]，表示为弧度制2.9585[rad]，将其数值代入公式⒀，算的第二组结果的纬度值为南纬10°30。

结合所求经度，我们得到符合条件的两个地区：菲律宾的吕宋岛附近（北纬17°18，东经1231612）和印度尼西亚努沙登加拉群岛附近（南纬10°30，东经1231612）。

（3）将直杆高度H作为检验因子进行模型检验

前面我们分别求出直杆的高度[H1]、[H2]、[H3]依次为1.0134米、1.2442米以及1.0440米。[H2]和[H3]的平均值为[E（H）=12（H2+H3）=1.1441]，将[E（H）]记为[H4]。则[H1]与[H4]平均值为[E（H）=12（H1+H4）=1.07875]，则[H1]与[H4]方差为[DH=12i=12（EH2-Hi）2=4.27×10-3]，可知，方差值很小，认为两种方法下计算的直杆高度近似相等，模型误差很小。

2.3 问题四模型建立与求解

2.3.1 问题四模型的建立

（1）对视频进行处理

首先利用[KMPlayer]软件对视频进行取帧，得到一组像素值为[1920×1080]的图片。因为视频时长40分钟，为保证截取的图片能准确反映整个视频中直杆影子的变化，我们对视频一分钟截取一帧，从8时56分3秒至9时34分4秒共截取39张图片。然后，对这一组图片进行灰度处理，这样可以使颜色细节更加细致，更容易辨认物体边缘。

（2）对直杆顶端、直杆影子顶端、直杆底端坐标的求解

利用[MATLAB]中的[Data Cursor Tools]分别选中你所观察到的直杆顶端、影子顶端、直杆底端记录它的像素点坐标，我们得到39组灰色矩阵，每一组灰色矩阵中存储了图片中的所有像素点的灰度值。以直杆顶端为例，在所记录的直杆顶端的像素点坐标附近寻找灰度值最小的像素点，记录它的坐标，即为直杆顶端的真实坐标点。求解直杆影子顶端、直杆底端的坐标点时也采用相同的方法。但对于直杆底端，我们认为底端相邻两个直角的y坐标的平均值为直杆底端的坐标。已知杆长可估计为2m，则可根据直杆顶端、直杆影子顶端、直杆底端的像素坐标得到前五幅图影子的坐标如下表所示：

表4 前五幅帧图中影子的坐标

[图片＼&1＼&2＼&3＼&4＼&5＼&[x]坐标（米）＼&2.307463＼&2.298507＼&2.265672＼&2.241791＼&2.229851＼&[y]坐标（米）＼&0.016418＼&0.01791＼&0.016418＼&0.016418＼&0.013433＼&]

以直杆底端为原点，以图片右侧方位为[x]轴，垂直图片向里的方向为[y]轴，建立直角坐标系。

图2 直杆底端坐标系的建立

（3）影子[y]轴数据加权确定

因为摄影机在录制视频的过程中，与地面存在一个固定的角度，因为镜像成像的原理，[y]轴方面人眼看到的影子长度小于实际中影子的[y]轴长度。因此引入[y]轴加权系数[C]，当摄像机角度不变时，人眼看到的影子的[y]轴长度与实际影子的[y]轴长度形成一个直角三角形的斜边和直角边，且比例保持不变，则在对本问题运用问题二的数学模型进行求解时，[y][坐标轴=通过灰色图像得出的坐标值×设定的权值C]。由于[x]轴坐标与摄像机观测方向垂直，[x]轴坐标可看作实际中的[x]轴坐标。

2.3.2 问题四模型的求解

当日期已知时：可求出当地时间为9时9分48秒，与之对应的北京时间为8时56分3秒。[经度=116+9+9/60+48/3600-8-56/60-3/360015=119.4o；][北纬=0.5537×180°π=31.7°]，所以该地点北纬314329，东经11924，约为江苏常州。

当日期未知时：可求出当地时间为9时06分45秒。此时，[N=173]，为2015年6月22日。地处北纬515313，东经1184030。

3 模型的评价与改进

问题一中，突出优点是考虑到北京当地时间和代表东八区的北京时间之间存在约14分的时差，并且应用在时角的计算中。

问题二中，能够运用解超定方程组的最小二乘法的思想，利用[MATLAB]的[fsolve]函数求得最优解。缺点是在寻求最优解的过程中，运行程序对初值十分敏感，当初值稍作变化时，所求的数值变化很大，这点不利于寻求最优解。在对此模型进行改进的过程中，考虑遗传算法和模拟退火模型，减小初值对结果的影响。

问题四中，突出的优点是考虑到人眼成像的特点，在直杆影子的[y]轴方向引入加权系数[C]，使数据更加准确，缺点是运行程序的过程中，不同的初值对结果几乎没有影响，观察不到某些影响因子的作用力。

参考文献：

[1] 谈小生， 葛成辉. 太阳角的计算方法及其在遥感中的应用[J]. 国土资源遥感，1995（6）.

[2] 王昌明. 可照时数和太阳高度角计算公式的简化证明[J].山东气象，1989（2）.

[3] 闫利， 胡修兵， 陈长军， 等. 中波红外几何标定场靶标定制问题研究[J].测绘地理信息，2013（8）.