根据视频中物体的太阳影子变化对物体进行定位的方法

摘要

本发明涉及一种根据视频中物体的太阳影子变化对物体进行定位的方法，包括步骤一：构建模型一与求解，构建试点物体的影子长度变化的数学模型：针对影响影子长度变化的相关变量：太阳高度角、所在地的经纬度、所在地的时间时刻，建立非线性模型，并针对太阳高度角、影长和物体高度的三角关系，建立等式，将模型应用于具体数据并绘制太阳影子长度的变化曲线；步骤二：用步骤一建立的模型，对视频中物体的太阳影子变化数据进行非线性拟合，得出视频中物体的经纬度。本发明技术方案，具有根据视频中物体的太阳影子变化能够对物体进行定位技术效果。

说明书

技术领域

本发明涉及一种根据视频中物体的太阳影子变化对物体进行定位的方法。

背景技术

如今的社会处于信息化时代，信息的交互与传递显得尤为重要。其中，太阳影子定 位技术将不断运用于大众之中。太阳影子定位技术根据物体的太阳影子的长度变化，准确 的判断物体的地理位置，既节省时间，又具有准确性，对于逐渐依赖电子产品的人们，这项 技术的发展具有巨大的潜力。

发明内容

针对现有技术存在的不足，本发明的目的在于提供一种根据视频中物体的太阳影 子变化对物体进行定位的方法。

为实现上述目的，本发明提供了如下技术方案：一种根据视频中物体的太阳影子 变化对物体进行定位的方法，包括步骤一：构建模型一与求解，构建试点物体的影子长度变 化的数学模型：针对影响影子长度变化的相关变量：太阳高度角、所在地的经纬度、所在地 的时间时刻，建立非线性模型，并针对太阳高度角、影长和物体高度的三角关系，建立等式， 将模型应用于具体数据并绘制太阳影子长度的变化曲线；步骤二：用步骤一建立的模型，对 视频中物体的太阳影子变化数据进行非线性拟合，得出视频中物体的经纬度。

本发明进一步设置为：所述步骤一与步骤二中用MATLAB进行求解，MATLAB程序设 置纬度的上下限为0°-90°，经度的上下限为0°-180°。

本发明进一步设置为：步骤一中所述太阳高度角用地理经度、时角以及太阳赤纬 近似表示：

$\alpha =\arcsin \left(sin\right(\phi ).sin(23.45sin\left(\frac{2\pi (284+n_1)}{365}\right))+cos(\phi ).$

$cos\left(23.45\sin \right(\frac{2\pi (284+n_1)}{365}\left)\right).\cos \left(15*\left(\left(t_1+\left(\lambda -120\right)\times \frac{360}{24}\right)-12\right)\right))$

其中，α为太阳高度角，φ为当地纬度，n1为该年距离1月1日的时间，t₁为具体的北 京时间，λ为当地经度；

物体影长和物体高度的关系可表述为

$l=\frac{h}{tan\left(\alpha \right)}$

其中1为影长，h为物体高度。

本发明进一步设置为：所述步骤二中根据最小二乘法建立模型如下目标函数和约 束条件：

目标函数：

$min=\underset{i=1}{\overset{21}{\Sigma }}\left(f_i\right(\phi ,\lambda ,l\left)\right)-L_{{\mathrm{yinzi}}_1}{)}^2$

约束条件：

$\left(\begin{array}{c}f(\phi ,\lambda ,l)=l/\tan \left(\arcsin \right(\sin \left(\phi \right).\sin \left(23.45\sin \left(\frac{2\pi (284+n_2)}{365}\right)\right)+\cos \left(\phi \right).\\ \cos \left(23.45\sin \right(\frac{2\pi (284+n_2)}{365}\left)\right).\cos (15*(\left(t_2+\left(\lambda -120\right)\times \frac{360}{24}\right)-12\left)\right)\left)\right)\end{array}\right)$

，φ∈(-90°，90°)，λ∈(0°，180°)，

其中，n₂为该年距离1月1日的时间，t₂为具体的北京时间，为已知测量点直杆 的太阳影子长度。

本发明进一步设置为：所述步骤二中根据最小二乘法建立模型如下目标函数和约 束条件：

目标函数：

$min=\underset{i=1}{\overset{21}{\Sigma }}{\left(f_i\right(\phi ,\lambda ,l,n)-L_{{\mathrm{yinzi}}_1})}^2$

约束条件：

$\left(\begin{array}{c}f(\phi ,\lambda ,l,n)=l/\tan \left(\arcsin \right(\sin \left(\phi \right).\sin \left(23.45\sin \left(\frac{2\pi (284+n_3)}{365}\right)\right)+\cos \left(\phi \right).\\ \cos \left(23.45\sin \right(\frac{2\pi (284+n_3)}{365}\left)\right).\cos (15*(\left(t_3+\left(\lambda -120\right)\times \frac{360}{24}\right)-12\left)\right)\left)\right)\end{array}\right)$

，φ∈(-90°，90°)，λ∈(0°，180°)，n∈(1，365)，

其中，n₃为该年距离1月1日的时间，t3为具体的北京时间，为已知测量点直杆 的太阳影子长度。

通过采用上述技术方案，地理位置的估计不仅仅自身是一个非常重要的问题，特 也能给其他问题带来不可预计的帮助。例如，已知某绑架犯发了一张受害人的照片或视频 给受害人家属，警方可运用太阳影子定位技术，通过太阳影子长度的分析，估计受害人所在 的位置，更快速、更有效得营救受害人。本文研究太阳影子定位技术，尽管实验结果表明太 阳影子定位技术具有很好的效果，但是还有许多工作需要进行，本文希望本研究工作仅仅 是地理信息化计算机时代的一个开端。

下面结合附图对本发明作进一步描述。

附图说明

图1为本发明实施例模型一的程序求解流程图；

图2为本发明实施例的太阳影子长度变化曲线图；

图3为本发明实施例的模型二的程序求解流程图；

图4为本发明实施例的模型二数值代入模型后的拟合图；

图5为本发明实施例的模型三数值代入模型后的拟合图；

图6为本发明实施例的模型四数值代入模型后的拟合图。

具体实施方式

问题一：建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律， 并应用建立的模型画出某时段天安门广场3米高的直杆的太阳影子长度的变化曲线。

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

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

问题四：建立确定视频拍摄地点的数学模型，并应用模型给出若干个可能的拍摄 地点。若拍摄日期未知，根据视频确定拍摄地点与日期。

本实施例被测物体优选为直杆，

符号说明：

t：地方时间，

φ：当地纬度，

λ：当地经度，

h：直杆的长度，

x₁：某个点的x轴坐标，

y₁：某个点的y轴坐标，

t_i：具体问题中的北京时间，

L_yinzi：直杆的太阳影子长度，

：已知测量点直杆的太阳影子长度，

G_i：影子坐标产生的角度，i＝1，2，

δ：纬角，又称为太阳赤纬，即太阳直射纬度，

n：该年距离1月1日的时间；即1月1号时值为1，1月2号值为2，

ω：时角，以正午12点为0度开始算，每一小时为15度，上午为负下午为正，即10点 和14点分别为-30度和30度，

α：太阳高度角，是太阳相对于地平线的高度角，是以太阳视盘面的几何中心和理 想地平线所夹的角度，

α_i：具体的太阳高度角，i＝1，2，

一、问题一模型构建与求解

要构建影子长度变化的数学模型，我们要寻找影响影子长度的各个参数及变化规 律。首先根据维基百科关于太阳高度角的计算，得出太阳高度角α，太阳时角，太阳赤纬δ， 当地纬度φ，存在以下函数关系：

方程(1)：sinα＝sinφ.sinδ+cosφ.cosδ.cosω

方程(2)：太阳赤纬：

$\delta =23.45\sin \left(\frac{2\pi (284+n)}{365}\right),(n=1,2,...,365)$

方程(3)：时角：

ω＝15*(t-12)(t＝0，2...24，t∈Z)

方程(4)：影子长度：

$l=\frac{h}{\tan \alpha }$

把方程(1)、方程(2)和方程(3)代入方程(4)，即可得到影

子长度模型：

$l=\frac{h}{\tan \left(arcsin\right(sin\phi sin\left(23.45sin\right(\frac{2\pi (284+n)}{365}))+cos\phi cos\left(23.45sin\right(\frac{2\pi (284+n)}{365}))cos\left(15\right(t-12)\left)\right))}$

由公式(1)可知，太阳高度角由地理经纬度、时角以及太阳赤纬决定。已知北京时 间9：00整，而北京时间测量位置为东经120°。已知直杆测量位置为东经116度23分29秒，则 直杆所处位置的时间与北京时间存在一定时间差，所以我们首先要对时间进行修正。由于 是地球自转影响时差，每24小时自转360°，则每小时自转15°。已知直杆测量位置的经纬度 为北纬39度54分26秒，与北京时间测量位置相差3度26分31秒。对方程(3)进行时间修订，该 测试点的地方时表达式可表示为：

方程(5)：$t=t_1+(\lambda -120)\times \frac{360}{24}$

(-180≤λ≤180)

t₁为该测试点北京时间；

方程(4)带入方程(3)可得出修正后的时角，修正后的该测试点的时角表达式为：

方程(6)$\omega =15*\left(\right(t_1+\left(\lambda -120\right)\times \frac{360}{24})-12)$

由方程(4)可知，直杆所处位置的时间与北京时间相差时，即地方时比北京时 间早小时，则直杆所处位置的实际时间为时。

方程(1)、(2)和(5)可得出太阳高度角可以用地理经度、时角以及太阳赤纬近似表 示：

$\alpha =\arcsin \left(sin\right(\phi ).sin(23.45sin\left(\frac{2\pi (284+n_1)}{365}\right))+cos(\phi ).$

方程(7)：$cos\left(23.45\sin \right(\frac{2\pi (284+n_1)}{365}\left)\right).cos(15*(\left(t_1+\left(\lambda -120\right)\times \frac{360}{24})-12\right)\left)\right)$

n₁为该问题测量日期距离1月1号的具体时间，

影长与太阳高度角和杆长成三角函数，则影长1和杆高h的关系可表述为：

方程(8)：$l=\frac{h}{tan\left(\alpha \right)}$

由问题资料已知得：

地理纬度：φ＝39.91°

地理经度：λ＝116.39°

北京时间：9：00～15：00

时间：n₁＝295

由上述具体数值，利用如下matlab软件程序联合方程(7)和方程(8)，模型一的程 序求解流程参见附图1，直接计算直杆的太阳影子长度，太阳影子长度的变化曲线绘制成图 形成“浴缸曲线”，参见附图2。

建立直杆的太阳影子长度的变化曲线模型的程序：

t＝9：05：15

phi＝39.907222222222222222224；

jd＝116.3913888888888；

n＝295；

delta＝23.45*sin(2*pi*(284+n)/365)；

w＝15*(t-(120-jd)/15-12)；

alpha＝asind(sind(phi).*sind(delta)+cosd(phi).*cosd(delta).*cosd(w))

yinzi＝3./tand(alpha)

plot(t，yinzi)

程序运行结果

t＝9.00009.500010.000010.500011.000011.500012.000012.5000 13.000013.500014.000014.500015.0000

alpha＝21.187825.507529.351232.616035.192736.976437.883137.8655 36.924635.109332.504729.216125.3526

1＝7.73946.28755.33484.68814.25393.98453.85603.85853.9920 4.26714.70825.36436.3315

“浴缸曲线”(又称U型曲线)是一个常用术语，用于描述半导体产品随时间变化的 瞬时故障率的通用曲线。尽管特定的半导体产品的实际故障率曲线差异很大，此通用曲线 仍然是一些通用属性的有用示例。这些是基于通用数据在最坏情况下的样本。由图2可知在 北京时间9点整至15点整时，9点整的太阳影子长度最长。上午时间，太阳影子由长逐渐变 短，下午由短逐渐变长，符合一天中影子客观变化规律。根据该图表的绘制，可观察该测量 点的各个时刻的影长。其中9点的太阳影子长度为7.7394米，15点的太阳影子长度为6.3315 米影子长度最短的时间大约在12点之后，12点半之前。以30分钟时间间隔观察，可得出出某 时段天安门广场3米高的直杆的太阳高度角和影子长度，如下表：

表1：具体时间的太阳高度角和影子长度

二、问题二模型的构建与求解

根据题目二给出的问题一的太阳影子顶点坐标数据，由于坐标系以直杆底端为原 点，水平地面为xy平面。由光影的形成原理可得垂直于地面的直杆的影子长度与xy坐标轴 满足勾股定理，影长l可表示为

方程(9)：

为该题由顶点数据推出的太阳影子长度。

由问题一可知阳高度角由地理经纬度、时角以及太阳赤纬决定，而太阳高度角和 杆长决定影长。该问题中有三个未知自变量，分别为该测试点的经度φ、纬度λ和杆长1。根 据最小二乘法建立模型如下目标函数和约束条件。

目标函数：

方程(10)：$\min =\underset{i=1}{\overset{21}{\Sigma }}{\left(f_i\right(\phi ,\lambda ,l)-L_{{\mathrm{yinzi}}_1})}^2$

约束条件：

$f(\phi ,\lambda ,l)=l/\tan \left(arcsin\right(\sin \left(\phi \right).\sin \left(23.45\sin \left(\frac{2\pi (284+n_2)}{365}\right)\right)+\cos \left(\phi \right).$

$\cos \left(23.45\sin \right(\frac{2\pi (284+n_2)}{365}\left)\right).\cos (15*(\left(t_2+\left(\lambda -120\right)\times \frac{360}{24}\right)-12\left)\right)\left)\right)$

φ∈(-90°，90°)，

λ∈(0°，180°)，

又由方程(7)和方程(8)可得三个未知参数之间的关系式为

$arctan\frac{1}{L_{{\mathrm{yingzi}}_1}}=arcsin(sin\phi .sin\left(23.45sin\left(\frac{2\pi (284+n_2)}{365}\right)\right)+cos\phi .$

方程(11)：$cos\left(23.45sin\right(\frac{2\pi (284+n_2)}{365}\left)\right).cos(15*(\left(t_2+\left(\lambda -120\right)\times \frac{360}{24}\right)-12)\left)\right)$

t₂为该地的拍摄时间，由该测试时间为4月18日，则n＝108。对模型进行最小二乘 法估计，拟合出未知参数的具体值，即确定问题二中的直杆所处的若干个可能的地点，并利 用如下matlab软件程序求解具体值。其中程序的解题程序流程图参见附图3。

对问题二中三个未知量的求解程序：

clc

x0＝[(rand()-0.5)*180，(rand()-0.5)*360rand()*10]％Startingguess

[x，resnorm]＝lsqnonlin(myfunB2，x0，[0，0，0.1]，[90，180，10])；％Invoke optimizer

x，resnorm

fid＝fopen(′xy.txt′，′r′)；

xy＝fscanf(fid，′％f′，[2，21])′；

fclose(fid)；

t＝14.7：0.05：15.7；

yinzi＝sqrt(xy(：，1).^2+xy(：，2).^2)′；

phi＝x(1)/180*pi；

delta＝23.45*sin(2*pi*(108+284)/365)/180*pi；

w＝15*((t-(120-x(2))/15)-12)/180*pi；

alpha＝asin(sin(phi).*sin(delta)+cos(phi).*cos(delta).*cos(w))；

F＝x(3)./tan(alpha)；

plot(t，F)

holdon

plot(t，yinzi，′r0′)

legend(′计算值′，′精确值′，′Location′，′SouthEast′)

xlabel(′时间t′)

ylabel(′影子长度′)

函数文件

functionF＝myfunB2(a)

fid＝fopen(′xy.txt′，′r′)；

xy＝fscanf(fid，′％f′，[2，21])′；

fclose(fid)；

t＝14.7：0.05：15.7；

phi＝a(1)/180*pi；

delta＝23.45*sin(2*pi*(108+284)/365)/180*pi；

w＝15*((t-(120-a(2))/15)-12)/180*pi；

alpha＝asin(sin(phi).*sin(delta)+cos(phi).*cos(delta).*cos(w))；

F＝(a(3)./tan(alpha)-sqrt(xy(：，1).^2+xy(：，2).^2)′)*10；

TXT文档

1.16373.336

1.22123.3299

1.27913.3242

1.33733.3188

1.3963.3137

1.45523.3091

1.51483.3048

1.5753.3007

1.63573.2971

1.6973.2937

1.75893.2907

1.82153.2881

1.88483.2859

1.94883.284

2.01363.2824

2.07923.2813

2.14573.2805

2.21313.2801

2.28153.2801

2.35083.2804

2.42133.2812

问题二程序运行结果：

x0＝-2.5232140.86767.6210

Maximumnumberoffunctionevaluationsexceeded；

increaseoptions.MsxFunEvals

x＝25.133195.66782.9730

resnorm＝0.0024

根据题目给程序设置不同的经纬度上下限会改变未知参数拟合后的具体值，同时 影响拟合的优度。经过不断地尝试，发现设纬度的上下限为0°-90°、经度的上下限为0°- 180°时，拟合程度最好。得出拟合图形参见附图4。

由于使用最小二乘法，当总体误差大于一定值时，程序自动输出满足最小的总体 误差值，所以我们通过改变误差的上限，求出可能满足问题一模型的影子顶点坐标数据变 化的具体测量点及杆长。

表2：可能地点及杆长

由于线性拟合误差不超过0.05，所以保留第1、2组的地理位置，因此，表一中的第 1、2组地理位置为问题二中直杆所处的可能位置。第1组所在地为海南省乐东黎族自治县， 第2组所在地为云南省德宏傣族景颇族自治州梁河县。第3组所在地为印度尼西亚。

三、问题三模型的构建与求解

根据固定直杆在水平地面上的太阳影子长短的变化，我们可以沿用问题二的模 型，并在此模型基础上上引入未知量日期n。该问题中就有四个未知自变量，分别为该测试 点的经度、纬度、杆长1和日期n。由方程(7)可得到

方程(12)：

$arctan\frac{1}{L_{{\mathrm{yinzi}}_2}}=arcsin\left(\sin \right(\phi ).\sin (\frac{2\pi (284+n)}{365}\left)\right)+\cos \left(\phi \right).$

$\cos \left(23.45\sin \right(\frac{2\pi (284+n)}{365}\left)\right).\cos (15*(\left(t_3+\left(\lambda -120\right)\times \frac{360}{24}\right)-12\left)\right)),$

t₃为该地的拍摄时间。

方程(13)：

方程(14)：$A_1-A_2=arcsin\left(\frac{sin\omega ·cos\delta }{{\mathrm{cos\alpha }}_1}\right)-arcsin\left(\frac{sin\omega ·cos\delta }{{\mathrm{cos\alpha }}_2}\right),$

方程(15)：

方程(16)：

在同一象限中，由(13)式得出(14)，且由影子坐标产生相应的，同时代入式 子(16)中，从而排除存在较大误差的太阳高度角，从而缩小待定系数可能值的范围。

问题三用问题二里仅给出太阳影子顶点坐标数据。为了估算影子长度，需要对上 式在MATLAB中进行非线性拟合，拟合出未知参数的具体值，即确定问题三中的直杆所处的 若干个可能的地点和日期，其中程序的解题程序流程图如图3。

方程(17)：$\cos \left(23.45\sin \right(\frac{2\pi (284+n)}{365}\left)\right).\cos (15*(\left(t_4+\left(\lambda -120\right)\times \frac{360}{24}\right)-12\left)\right))$

对该测试点的经度、纬度、杆长h和日期n四个变量假设为待定系数，进行在MATLAB 中进行非线性拟合。t₄为该地的拍摄时间。根据最小二乘法建立模型如下目标函数和约束 条件。

目标函数：

方程(18)：$\min =\underset{i=1}{\overset{21}{\Sigma }}{\left(f_i\right(\phi ,\lambda ,l,n)-L_{{\mathrm{yinzi}}_1})}^2,$

约束条件：

$f(\phi ,\lambda ,l,n)=l/tan\left(arcsin\right(sin\left(\phi \right).sin\left(23.45sin\right(\frac{2\pi (284+n_3)}{365}\left)\right)+cos\left(\phi \right).$

$\cos \left(23.45\sin \right(\frac{2\pi (284+n_3)}{365}\left)\right).\cos (15*(\left(t_3+\left(\lambda -120\right)\times \frac{360}{24}\right)-12\left)\right)\left)\right),$

φ∈(-90°，90°），

λ∈(0°，180°)，

n∈(1.365)，

根据问题二中的数据，在MATLAB中，进行非线性拟合，可知固定直杆所处的地点和 日期可能为

表3：可能地点及杆长

由于需要筛选线性拟合误差较小的数据，所以保留以上6组的地理位置，因此，表 二中的第1-6组地理位置为问题二中直杆所处的可能位置。第1组所在地为新疆维吾尔自治 区和田地区墨玉县，测量日期为7月27日。第2组所在地为新疆维吾尔自治区克孜勒苏柯尔 克孜自治州乌恰县，测量日期为5月18日。第3组所在地为新疆维吾尔自治区和田地区于田 县，测量日期为7月22日。第4组所在地为古尔伯加印度卡纳塔克邦，测量日期为3月2日。第5 组所在地为新疆维吾尔自治区阿克苏地区沙雅县，测量日期为7月2日。第6组所在地为哈萨 克斯坦，测量日期为6月10日。

选取上述第1组的结果，做代码修正，得出拟合图形参见附图5。

通过附图5发现拟合程度良好，预测得出的具体地点位置角，可信度较高。

四、问题四模型的构建与求解

求解视频的拍摄地点，要明确题中所给数据以及未知量个数。首先要分解视频，得 出影子长度的坐标，主要步骤为：

1、在MATLAB软件利用VideoReader函数读取该时刻直杆在太阳下的影子变化视 频：VideoReader(‘Appendix.avi’)，对视频每三分钟提取一张画面，该画面反映了影子的 具体位置。

2、保存图片，imwrite(img，filename)。

3、调用imread函数读取每一帧图片，imgRgb＝Imread(filename)。

4、用img2gray函数将彩色图片转化成灰度图片，rgb2gray(imgRgb)，灰度表明了 图像像素的明暗程度。

5、对每一页的灰度图像，将其放在(x，y)平面上，将图像分成均匀的坐标，以便于 找到图片中影子的顶点及坐标。

由于现实直杆长度为2米，根据实际直杆长度与图片直杆长度的比例关系：

$\frac{L_{yinzi}}{h}=\frac{L_{{\mathrm{yinzi}}_3}}{h_1}$

得出实际影子长度，经过处理得到时间与太阳影子长度的相关数据：

表5：时间和处理后的实际太阳影子长度

已知视频的拍摄时间，确定包含两个未知量，求解拍摄的地理位置，与问题二求解 相似，将数据代入问题二中的模型，进行非线性拟合，估算出影子长度，若拟合效果越好，则 表示拟合后的影子长度与实际影子长度误差越小，则地理位置越优。可知视频中直杆所处 的地点可能为：

表6：模型四的具体地点的经纬度与误差

未知视频的拍摄时间，确定包含三个未知量，求解拍摄的日期和地理位置，与问题 三求解相似，将数据代入问题三中的模型，进行非线性拟合。可知视频的拍摄时间和地点可 能为

表7：可能地点与时间

由于线性拟合误差不超过0.05，所以保留第1、2、3、4组的地理位置。因此，表三中 的第1、2、3组地理位置为问题二中直杆所处的可能位置。第1组所在地为内蒙古自治区包头 市固阳县，拍摄日期为5月21日。第2组所在地为内蒙古自治区呼和浩特市武川县，拍摄日期 为6月7日。第3组所在地为河北省保定市易县，拍摄日期为7月4日。第4组所在地为山东省潍 坊市临朐县，拍摄日期为8月22日。

同时对其进行代码修正，绘制图表，拟合程度较好。拟合图形参见附图6。