标签：

太阳影子定位问题研究

摘要

本文主要研究影长与日期、时间及地点的关系。通过其中某组数值，去刻画出其他变量的大致范围。

对于第一问，由于日地距离随公转而改变，故我们需要寻找出一个合适的方程去描述不同日期下的日地距离，从而建立合理的空间坐标系，确定当下太阳在坐标系下的空间坐标。又根据给定地点的经纬度，我们可以建立影长方程，从而用MATLAB绘画出影长图像。

对于第二问，所给数据的影长变化，可以推测出测量地处于下半天（正午12时至凌晨0时），由此可以推算出不符合要求的地区，从而缩小数据量。利用第一问所推导的影长公式，使用MATLAB算出符合上述条件的各地区影长，与所给数据的影长作标准差，来看出不同地区的影长数据的偏离程度。从而确定出更为合理的测量地。同时可以根据给定的影长坐标，绘画出影长变化曲线，并用二次函数进行拟合，求出影长最短的时刻，通过所给数据的时刻，来估测出测量地的时区或经度，并与我们推理的测量地进行比较。经过推算，得知测量地位于东经108.67°，北纬19.26°，也就是我国的海南省东方市附近。

对于第三问，与第二问比较起来缺少变量日期。与第二问类似，我们先推测出可行的经纬度及日期范围，再对此经纬度、日期所对应的影长进行求解，并与实际观测数值进行对比。由于相对多了一个自由量，所以相应解的个数会增加，对所有可能解我们再利用图像进行进一步检验。根据模型构建我们得到第一组数据测量地位于东经79°，北纬39°，此地为新疆喀什附近，日期为5月12日或7月31日。对于第二组数据测量地位于东经110°，北纬34°，此地为陕西省商洛市附近，日期为10月26日。

对于第四问，视频信息每间隔3min进行截图，根据图片使用Digimizer软件测量出影长。与第二问比较起来，增加了杆长。也就是说，可以求出不同时刻的太阳入射角。根据第一问的太阳入射角与地点经纬度的关系式，通过所测数据的非线性拟合，确定出测量地为东经111°，北纬43°，在内蒙古乌兰察布市附近。

关键词： 空间坐标系 二次函数拟合 标准差 非线性拟合

1、问题重述

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

1.建立影子长度变化的数学模型，分析影子长度关于各个参数的变化规律，并应用你们建立的模型画出2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线。

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

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

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

如果拍摄日期未知，你能否根据视频确定出拍摄地点与日期？

 附件一：影长数据（含3个表格）

 附件二：视频

2、问题假设

1、同一天中不考虑地球公转所引起的太阳直射点纬度变化

2、不考虑各地的海拔、高地等引起的高度变化

3、照射到地球表面的太阳光线为平行光

4、题中所给的所有数据准确无误

5、视频所截取的数据能够具有代表性

3、符号说明

4、模型建立与求解

4.1模型一的建立与求解

  4.1.1 模型一的建立

 以地心为坐标原点，指向赤道上东经0°点方向为x轴正方向，指向北极点方向为z轴正方向，垂直于x轴、z轴的方向为y轴正方向，构建右手直角坐标系。记测量地点的经纬度为,地球半径为R，那么测量地的位置坐标为，那么直杆的方向向量可以取为，日地距离记为，

直杆长度为，影子长度为，太阳直射点的经纬度为。

其中，  即

太阳直射点：

经度，是测量地点的时间，是地球的自转角速度，。

纬度，是从春分日开始经过的天数，                                               [1]

记日地平均距离为，又称为天文单位，其中1天文单位

可以精确至。

由于日地距离对于任何一年的任何一天都是精确已知的，所以这个距离可用一个数学表达式表述。为了避免日地距离用具体长度计量单位表示过于冗长，一般均以其与日地平均距离比值的平方表示，即。               [2]

我们得到的数学表达式为：   

式中称日角，即，这里又由两部分组成，即，

其中为积日，所谓积日，就是日期在年内的顺序号，例如，1月1日其积日为1，平年12月31日的积日为365，闰年则为366，等等。

又有，

其中为数值向下取整为最接近的整数的函数。

因此我们可以得到了观测点所在地的日地距离,.

记太阳入射角，太阳的坐标为，其中

所以太阳入射角为与的夹角的余角，即，此时，若则舍去。

故可以得到影子长度:

。

  4.1.2 模型一的求解：

由于观测点为天安门广场，从而可以得到观测点的经度为东经116度23分29秒，纬度为北纬39度54分26秒，即，

又根据地球半径，从而可以得到北京天安门广场坐标为

根据观测时间2015年10月22日北京时间9:00-15:00可以得到

，，

，

那么此时的日地距离。

也就是说，在北京时间2015年10月22日太阳直射点位于南纬

接着根据公式，我们可以得到不同时刻太阳直射点的纬度，从而可以得出不同时刻太阳的位置坐标，（见下表）

由于已知杆长，所以根据入射角公式及影长公式可以得到不同时刻的太阳入射角与影子长度（见下表）。

为了绘制出影子长度变化曲线，我们对函数中的21组数据进行MATLAB制图。

程序：

R=6371;

time=9:0.3:15;

a=116+23/60+29/3600;

b=39+54/60+26/3600;

c=120-15.*(time-12);

d=asin(0.39775*sind(180+(214-186)));

r=1.48887659888219*10^8;

A=[R*cosd(b)*cosd(a) R*cosd(b)*sind(a) R*sind(b)];

B=[r*cos(d).*cosd(c); r*cos(d).*sind(c);r*sin(d).*ones(1,21)];

for i=1:21

ci(i)=pi/2-acos(dot(A,B(:,i))/norm(A)/norm(B(:,i)));

end

yz=3./tan(ci);

plot(time,yz)

运行结果：

    至此，我们得到了2015年10月22日北京时间9:00-15:00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线。

4.2模型二的建立与求解

  4.2.1 模型二的建立

根据所给定的21组不同时间的影子顶点坐标数据，绘画出影子长度随时间的变化曲线。

此时我们可以使用二次函数拟合这组数据，并求出最低点的时间。

程序1：

a=[14.7 1.149625826;14.75 1.182198976;14.8 1.215296955;14.85 1.249051052;14.9 1.28319534;14.95 1.317993149;15 1.353364049;15.05 1.389387091;15.1 1.426152856;15.15 1.463399853;15.2 1.501481622;15.25 1.540231817;15.3 1.579853316;15.35 1.620144515;15.4 1.661270613;15.45 1.703290633;15.5 1.74620591;15.55 1.790050915;15.6 1.835014272;15.65 1.880875001;15.7 1.927918447];
A=polyfit(a(:,1),a(:,2),2);
low=-A(2)/2/A(1)
low =

   12.5984

根据运行结果，可以估算出北京时间12.5984小时，测量地为影子长度最小时刻，也就是太阳高度最高时刻，即当地时间的12时。

       为测量地的经度。

数据的时间是北京时间14.7时到15.7时。此时的影长是递增的，可以知道测量地此时处于下午时刻。

，代表测量地在北京时间15.7时才刚到下午。

，代表测量地在北京时间14.7时才刚到凌晨。

我们可以得知，从西经到东经处不可能是测量地。

太阳直射点的经纬度是，则记，代表北京时间4月18日，以此纬度为界限向南全天均为黑夜。

由，

则，  

那么南纬到在北京时间4月18日一直处于黑夜。

   接着根据公式，我们可以得到不同时刻太阳直射点的纬度，从而可以得出不同时刻太阳的位置坐标，（见下表）

由于影长，这里分别代表第一小问的杆长和太阳入射角。

又有，对于每一固定的地点是确定的，所以杆长不确定的情况下，我们带入杆长时，与实际测量值相差倍。计算各地的影长与测量值之间的标准差，使得标准差最小的的地点，也就是每组数据与测量值相差的倍数越靠近同一个值得地点，为测量地的可能性就越大。

  4.2.2 模型二的求解

利用MATLAB编写程序可以得到各地的影长与测量值之间的标准差。

程序：

z=zeros(19700,17000);

for i=64:0.01:260

for j=-80:0.01:89

for k=1:21

c(k)=120-15*(a(k,1)-12);

dt=asin(0.39775*sin(pi*(28+(a(k,1)-6.75)/24)/186));

l=[cosd(j)*cosd(i),cosd(j)*sind(i),sind(j)];

m=[cos(dt)*cosd(c(k)),cos(dt)*sind(c(k)),sin(dt)];

ac=acos(dot(l,m));

if(ac>pi/2)

z((i-64)*100+1,(j+80)*100+1)=1;

end

t(k)=1/(a(k,2)*tan(pi/2-ac));

end

if(z((i-64)*100+1,(j+80)*100+1))==0)

st((i-64)*100+1,(j+80)*100+1)=std(t);

else

st((i-64)*100+1,(j+80)*100+1)=1000;

end

end

end

运行结果：

由于运行结果过于庞大，我们截取了部分结果，同时利用EXCEL求出了标准差的最小值以及相应的点。

因此得到了最小的标准差为1.34012E-05，又由于在该点附近的标准差与最小值相差甚小，故可推测出直杆所处的地点为东经108.67°，北纬19.26°。也就是我国的海南省东方市附近。

  4.2.3 模型二的检验

计算值与测量值的影长曲线图，如下：

根据上表的数据绘制误差变化图，如下：

4.3模型三的建立与求解

4.3.1 模型三的建立

对于这两组数据，以第一组数据为例。

影子长度由四个变量所控制：日期，时间，地点(经度、纬度)。

我们使用相同的符号：观测点经纬度,太阳坐标(方向)。利用蒋洪力论文的结论，根据太阳运行进度来求太阳坐标的纬度.得到：

其中为每年的积日。

,

利用之前得到的公式，得到太阳入射角与观测点，时间的函数关系：

其中

.

因此影长.

因此应为常值.(为测量值)

我们再次利用的循环变量，即对经度、纬度、日期作为三轴，建立三维网格，我们对网格所有格点进行一次遍历，并选择出最优解(即的波动范围较小者，即方差较小者).

在循环中，我们利用二次拟合，求出测量地的大致精度范围.

4.3.2 模型的求解

对于第一个测量地，我们利用二次拟合得到曲线

所以观测地的大致范围在

对于第二个测量地，我们利用二次拟合得到曲线

所以观测地的大致范围在

程序：

(1)

for i=61:120

for j=-90:89

for k=1:365

yz=pnt(i,j,k,a(:,1));

t=yz‘./a(:,2);

st((i-60)+(k-1)*60,j+91)=std(t);

end

end

end

(2)

for i=78:137

for j=-90:89

for k=1:365

yz=pnt(i,j,k,a(:,1));

t=yz‘./a(:,2);

st((i-77)+(k-1)*60,j+91)=std(t);

end

end

end

对于第一个测量地，得到如下的运行结果：

相对应的结果为11月25日 东经80° 南纬39°.

在世界地图上寻找对应点，发现为海上，所以此值应当舍去，选取次优解。

次优解有两个：5月12日 东经79° 北纬39°。

7月31日 东经79° 北纬39°，此地为新疆喀什附近。

对于第二个测量地，得到如下的运行结果。

   相应的结果为10月26日 东经110° 北纬34°，此地为陕西省商洛市附近。

4.4模型四的建立与求解

  4.4.1模型四的建立

坐标系与各变量的表示均与问题一中相同。

根据附件4中直杆在太阳下的影子变化的视频，以3分钟为一节点截取视频中的图像，利用Digimizer测量软件测量出影子长度，共测得14组数据（如下表）。

记观测地A的经纬度为（a,b）,那么直杆的方向向量可表示为.

接下来考虑太阳直射点的经纬度：

    经度，是测量地点的时间，是地球的自转角速度，。

纬度，是从春分日开始经过的天数，。

根据视频我们可以得到测量地点的日期和时间，从而和可知，即不同时刻太阳直射点的经纬度可知。

从而太阳位置可表示为。

所以太阳入射角，影子长度。

即

也就是，

  4.4.2 模型四的求解

首先根据上述模型，得到不同时刻太阳直射点的经纬度，并且求出太阳的坐标。

程序：

a=[8.9 2.4097;8.95 2.3645;9 2.3195;9.05 2.2757;9.1 2.2416;9.15 2.1917;9.2 2.1555;9.25 2.1062;9.3 2.0759;9.35 2.0315;9.4 1.9947;9.45 1.9525;9.5 1.9191;9.55 1.8792];

for i=78:137

for j=-90:89

yz=pnt(i,j,194,a(:,1));%pnt为求影子长度的函数，见附录。

t=yz‘./a(:,2);

st(i-77,j+91)=std(t);

end

end

end 

运行结果：

利用EXCEL查找出了标准差的最小值，所以可以得知测量地为东经111°，北纬43°，在内蒙古乌兰察布市附近。

  4.4.3 模型四的检验

a=[8.9 2.4097;8.95 2.3645;9 2.3195;9.05 2.2757;9.1 2.2416;9.15 2.1917;9.2 2.1555;9.25 2.1062;9.3 2.0759;9.35 2.0315;9.4 1.9947;9.45 1.9525;9.5 1.9191;9.55 1.8792];

st=zeros(60,180);

for i=78:137

for j=-90:89

yz=pnt(i,j,194,a(:,1));

t=2.*yz‘./a(:,2);

for k=1:14

st(i-77,j+91)=st(i-77,j+91)+abs(t(k)-1);

end

end

End

容易得到：影长的误差和为   ，

其中为东经111°，北纬43°在对应时间下的影长；为对应时间下的测量值。

则，

即

 4.4.4 模型四的检验模型四的推广问题

在拍摄日期未知的情形下，那么这个问题就类似于问题三。我们可以通过同样的方法，估算出一些可能的地点，但是问题三利用的是公式算得的影长与测量值之间存在一个杆长倍数的关系，然后利用计算的数值，比较其标准差，考察每个数据与真实值比值的离散程度。但是问题四已经给出了一个杆长，所以我们可以继续进行检验，观察在这些可能的地点中这个比值是否在杆长附近波动。波动程度越小，那么地点是测量地的可能性就越高。这样也就可以估测出大致的拍摄地点及日期。

5、模型的评价

  6.1 模型的优点

  （1）数学推导，具有较强的理论依据，逻辑的推理增强了数据间的连贯性，方便观察出其中明显的失误。

  （2）使用拟合数据，进行测量地的位置估算，可以对结果起到引导性作用，更容易做出精确的结果。

  6.2 模型的缺点

  （1）由于地理因素复杂，无妨用简单的数学模型完美刻画，得到的结果也只能是在不考虑众多地理复杂因素下的估计值。

  （2）由于不同的拟合方式可以得到不同的精确度，所以引导性的精确程度也是可以不断优化的，然而我们选择的拟合方式并不一定是最好的。

6、参考文献

[1]  [0583-1458] 蒋洪力，太阳直射点纬度的数学推导和分析，数学通报，第9期:P39-40，2007。

[2]  [TK511;P182.1] 王炳忠，太阳辐射计算讲座 第一讲 太阳能中天文参数的计算，太阳能，第2期:P8-10，1999。

[3]  [978-7-5643-2992-1] 龚涛，摄影测量学，成都：西南交通大学出版社,2014.04。

[4]  [9787040196382] 姜启源，谢金星，数学建模案例选集，北京：高等教育出版社，2006。

[5]  [9787302265535] 李学文，数学建模优秀论文精选与点评，北京：清华大学出版社，2011。

太阳影子定位问题研究

标签：

原文地址：http://blog.csdn.net/ycy0312/article/details/51191214