太阳影子定时定位的技术研究

太阳影子定时定位的技术研究

摘要

本文的目的在于利用视频中的物体在太阳光下的阴影长度与角度的变 化分析物体所处的拍摄地点与拍摄时间，从而进行视频数据分析。

我们进行合理的假设，运用几何分析、向量分析、建立坐标系将研究 目标解析化的方法，利用太阳直角坐标系与地平坐标系之间的转换，建立 了影长变化与日影定位定时模型。整合函数形式后，对已知数据进行预处 理，编程进行曲线拟合求出待定参数。搜寻大量数据检验模型的准确性，

依据现实对模型进行修正，得到太阳影子定时定位分析视频的方法，并有 效地对具体视频进行了运用。

对于问题 1，我们建立直杆影长变化模型，结合地球公自转运动进行 几何分析，利用向量建立起各个参数与影长之间的函数关系，再利用建立 的模型对具体时空下的直杆影长进行函数分析，并作出变化曲线。

对于问题 2、3，我们建立日影定时定位模型，先创建太阳直角坐标 系与地平坐标系，研究两者之间的转换关系，再将直杆与光线解析化，得 到影子顶点的坐标函数。对不同时间测量的影子坐标数进行预处理，将坐 标函数整合为向量值函数，在 Matlab 中用最小二乘进行曲线拟合，求出 待定的参数值，再分析拍摄时间与地点的可能组合解。

对于问题 4，我们运用视频处理技术将原彩色视频处理为黑白可读点 的形式，建立像素坐标，读取大量杆与影的坐标数据，处理后代入模型中 进行运算与比对，再对时空进行有效分析。

最后，我们搜查大量真实的测量数据，对它们进行处理，代入模型后 进行检验，对产生的误差作图分析，并寻找误差来源。依据误差对模型进 行了优化与修正。

关键词：日影定时定位，解析化，赤纬角，太阳高度角，最小二乘拟合

一．问题提出

视频数据分析中，视频的拍摄地点与日期是很重要的信息，如何由视 频本身蕴涵的内容推断出视频拍摄的地点与时间是有必要研究的课题。

其中，利用视频中物体在太阳光下阴影的长度、形状、与物体夹角、

顶点位置坐标等信息，进行合理的假设，建立适当的数学模型，可以有效 地分析出视频的拍摄地点与拍摄时间。

这种利用物体太阳阴影变化的信息进行视频分析的方法有广泛的应 用，同时也操作简便，具有很多的优处。在日常生活中，通过观察周边物 体在日光下的影子变化也可以推测出所处的时间与地点，从而在勘测与地 理中具有很多的用途。

二．问题分析

为了推断视频的拍摄地点与拍摄时间，我们应该研究如何利用物体 在太阳下的阴影变化分析出物体所在的经度纬度，以及当下所处的时间。

我们不妨以垂直于地面的直杆为研究对象。

首先，对问题一，我们研究在一年任何时间中，在任何经度纬度下，

直杆在太阳下的阴影长度如何变化。影子的长度由物体的实际高度和当时 的太阳高度角决定，太阳高度角由当日的赤纬角和当时的地球时角决定，

建立太阳高度角与地球公自转运动之间的关系即可知道测量时的太阳高度 角。最后可以利用建立的模型对给定的具体时间地点下的直杆影子长度进 行具体的函数分析，并作出变化曲线。

其次，对问题二三，我们研究在测定了一组某天内不同时间的直杆 影子顶点位置后，如何建立合适的数学模型，推断出直杆所处的地点经纬 度和测量日期。我们可以利用空间解析几何或向量的方法将影子解析化，

先建立太阳直角坐标系，再转换成当地的地平坐标系。因为题目中没有给 出坐标轴的方向，我们可以先假定 X-Y 的标准定向，再分析与题中坐标相 差的旋转角度。

最后，我们可以利用以上建立好的模型对具体的视频进行分析，处理 视频从中得到可用的数据信息，用建立好的模型分析拍摄地点与拍摄日期 的组合解。

根据实际情况，我们可以对模型进行误差分析，并做出更加符合实际 的修正与优化。

三．基本假设

 不考虑大气层中的水蒸气，颗粒，尘埃等对太阳光线的折、散射，

假设太阳光线射到地球时没有经过弯折，是太阳中心与地心的连线

 假定太阳光射到地球时是平行光。因为地球与太阳之间的距离比地 球半径大许多数量级

 假设地球公转轨道在同一个平面内，不考虑地球太阳外的其它星体 对地球公转产生包括摄动、岁差等的影响

 模型一假定地球围绕太阳运动的轨迹是圆形，地球公转是匀速的。

因为公转的椭圆轨道的偏心率是 0.010643，接近于 0

 假设正午太阳高度角出现的时刻是当地时间中午十二点

 假设一天之内太阳直射点的纬度不发生改变

 不考虑日光衍射对影子长度产生的影响，影子边缘清晰可测，因为 测量的直杆都不太高，“表高景虚”现象不明显

四．符号说明

：直杆所处地点的纬度

：直杆所处地点的经度

：太阳赤纬角，即地球中心与太阳直射点连线与地球赤道平面的夹角

：黄赤交角(23°26’21”)

：夏至日到当前日期地球的公转角度 ℎ：太阳高度角，即太阳光线与地面的夹角

：时角，即正午太阳高度角出现的时刻到测量时刻地球自转的角度 H：直杆的竖直高度

L：直杆的太阳影子长度

：测量影长时的当地时间

：测量影长时的北京时间 n：空气折射率

五．直杆影长变化与日影定位定时模型

5.1 直杆影长变化模型

我们先来研究在一年中的任何时间，任何一个经纬度坐标上，直杆在 太阳光照射下的影子长度如何变化。直杆影子的长度与直杆的高度，太阳 直射点的经纬度，直杆所处的经纬度有关系。在模型一里假定地球公转轨 道是圆形，地球匀速公转。

5.1.1 求某日期里太阳直射点的纬度

如图 1，A 点是夏至日时太阳的直射点，位于黄道面上，只有夏冬至日 时直射点 A 与地心 O 的连线和地球自转轴决定的平面 AON 与黄道面垂直。

夏至日到当前测量日期地球的公转角度为 ，等效地可认为太阳直射点与 地心的连线转过了 角，太阳直射点由 A 变为 B。

由符号设定，∠BON＝ － ，∠AON＝

－ ，∠AOB＝ ,由几何关系：

sin ＝sin ∗ cos （1）

假设地球公转轨道是圆形，公 转速率不变，已知某日期与夏至日 的时间间隔后，可以算出从夏至日 到这一天地球公转的角度，即 ， 又黄赤交角 已知，则可以得到这一 天里太阳直射点的纬度

5.1.2 求直杆所在位置的太阳高度角ℎ

如图 2，以地球自转 轴为 Z 轴，正午时的光线 与自转轴决定的平面为 ZOY 平面，从而建立三维 坐标系。用 ⃗表示光线向 量， ⃗表示直杆向量。直 杆在正午时位于 ZOY 平面 内，已知测量时刻后，由 于地球匀速自转，可以得 到正午时刻到测量时刻地 球自转的角度

在建立的三维直角坐标下分别写出光线与直杆的方向向量：

⃗＝（0，cos ，sin ）

⃗＝（cos sin , cos cos ，sin ） 则由内积公式< ⃗, ⃗ >= | ⃗| ∗ | ⃗| ∗ sin ℎ算出：

sin ℎ＝_cos cos cos ＋sin sin (2) 当地时间为中午十二点时，太阳直射点与直杆在同一经度线上，此时

＝0，测量影长时的当地时间 t＝(直杆所在经度－北京时间对应的经 度)/15°＋测量时刻的北京时间 t^o，则测量时刻的 ＝(t-12)∗ 15° (3)

5.1.3 求直杆在阳光下阴影的长度

给定具体的日期，可以求出夏至日到当前日期地球的公转角度

，由式(1)可以求出该天太阳直射点的纬度；给定当天中的具体时间，又 已知直杆所在的经纬度，可以由式(3)得出正午时刻到测量时刻地球自转

的角度 ，再将(1),(3)代入(2)，则可以得到太阳光线与地面的夹角ℎ，

则 H 表示直杆的高度，得到

= ∗ cot ℎ (4) 5.1.4 利用建立的模型求解问题 1

因为所求日期(10 月 22 日)更接近秋分日(9 月 23 日)，为了使结论更 准确，我们用秋分日到所求日期地球公转的角度+90°代替 ，则 ＝_{90° +}

∗ 360°，代入式(1)得到 ＝－10.9778°；又已知直杆所在经度 ^E、当地 时间 ，代入式(3)求得时角 ；将 、 、 代入式(2)求得太阳高度角ℎ；

再带入式(4)求得直杆影子长度 L，以 6 分钟即 0.1 小时为间隔，由不同 的北京时间 计算得到不同的影子长度 L，具体数据表格见附录一。

由表格作出影长随时间变化的函数曲线如图 3，程序见附录二。

从北京时间 9:00 到 15:00，影长曲线先从 7.40 米单调递减到 3.69 米，再单调递增到 6.08 米，在 12:14 达到最小值点。

5.2 日影定时定位模型

5.2.1 日影定位模型求解问题二

如左图 4，以地心为原点 O，地球自转轴为 Z 轴，直杆与 Z 轴决定的平面为 ZOY 平面，

建立太阳直角坐标系。以直杆 所在地点为 O’原点，直杆为 Z’轴，X’轴与 ZOY 平面垂直 指向正东，建立地平坐标系。

建立两个坐标系后，我们 将光线、直杆、影子解析化，

为了求影子在地平坐标系中的 横纵坐标表达式，应该先求光 线方向向量在地平坐标系中的表示。

在太阳直角坐标系中，直杆的方向向量为 ′⃗＝（0, ,sin ） 当地时间中午十二点时，太阳光线 ⃗位于 ZOY 平面内， ⃗＝（0，

−cos ，− sin ）从当地十二点到测量时刻，地球自转了 角度，相当于 地球不动，太阳光线向西转动 角度，光线变为：

′ ⃗＝（− cos sin ，－cos cos ，－sin ） (5)

地平坐标系中的三个单位向量基在太阳坐标系下的坐标表示如下:

1′⃗ = （ − 1，0，0）

3′⃗ = （0，cos ，sin ）

2′⃗ = 3′⃗ × 1′⃗ =（0，－sin ， cos ） 那么光线 ′ ⃗在地平坐标系下的三个坐标投影为：

(6)

＝ < ′ ⃗, 1⃗ >＝cos sin

=< ⃗, 2⃗ >＝cos cos sin － sin cos

=< ⃗, 3⃗ >＝− cos cos cos − sin sin 现在转换视角，在地平坐标系中计算直杆的影子坐标，如图 5:

设影子的坐标为（ ’， ’，0）由相 似关系得：

= ′

= ′

(8) 将(7)代入(8)，求得直杆影子顶点的 坐标为：

= cos sin

cos cos cos + sin sin

= H(cos cos sin － sin cos ) cos cos cos + sin sin ) 由(9)得到影长 L、影子与 x’轴的夹角正切，它们可以等效地表示直 杆影子顶点的坐标位置：

L＝

’ + ＝

H (

cos cos cos + sin sin

)

⁻² − 1

(10)

′

′

^{＝ sin cot － cos tan}

1

sin (11)

其中 = ^{北京时间对应的经度}

° + − 12: 00 ∗ 15° (12) 先将(12)代入(10)，已知测量日期为 2015 年 4 月 18 日，由此算出测 量日期的赤纬角 ＝10°37 3′′ ，视测量时刻的北京时间 为自变量，影长

L

为因变量，

H、

、 为待求的参数。用附件一给出的数据算出每组

−

L的

值(数据表格见附录三)，以(10)的形式在 Matlab 中用

（9）

(7)

lsqcurvefit 函数进行最小二乘拟合，程序见附录四，求出各个待定的参 数值为：

H＝2.034 米

＝19°17′37′′

＝108°43′14′′

为了有效利用数据，再将(12)代入(11)，同样视测量时刻的北京时间 为自变量，夹角

tan

^{为因变量， 、} 为待求的参数。用附件 1 的 数据算出每组的夹角值（数据表格见附录三）进行最小二乘拟合，程序见 附录四。由于与我们假定的 x’-y’轴是正东与正北，是题目中的 x-y 坐 标轴旋转一定的角度得到，拟合得到的函数形式与(11)求反正切后相差固 定角度。由此求出待定的参数值为：

＝18°17’37”

＝109°28’28”E

Y’轴与正北的夹角为15.8°

比较由影长和夹角正切得到的(13)与(14)，为了减小误差，取数值的 平均值，得出结论：

直杆的高度 H 为 2.034 米，位置在（_{18°47 37} ，_{109°5′51′′} ）

5.2.2 日影定时定位模型求解问题三 A.求解附件二

依然使用 5.2.1 中的模型，同样得到式(10),(11)，测量时刻的北京

时间 为自变量，影长

L、

夹角

tan

为因变量，此时由于测量日期未 知， 、

H、

、 均是待求的参数。

(13) )

(14) )

先将(12)代入(10)，用附件二给出的数据算出每组 −

L的

^值（数据

表格见附录五），以(10)的形式在 Matlab 中用 lsqcurvefit 函数进行最 小二乘拟合，程序见附录六，求出各个待定的参数值为：

H=2.008 米

＝39°53′33′′

＝79°44 38

＝_{20°45’38”}

再将(12)代入(11)，用附件 1 的数据进行最小二乘拟合（数据表格见 附录五），程序见附录六。由于与我们假定的 X’-Y’轴是正东与正北，

是题目中的 X-Y 坐标轴旋转一定的角度得到，拟合得到的函数形式与(11) 求反正切后相差固定角度。由此求出待定的参数值为：

＝39°31’42”

＝79°44’55”E

＝− 20°19’13”

Y’轴与正北夹角为－110°

(15)(16)得到的数据中纬度与赤纬角恰好都相反，这是因为在 L 关于 时间 t 的函数(10)中，观察函数形式发现纬度与赤纬角都取相反数后却可 以得到相同的影长 L 值，但tan 关于时间 t 的函数得出的解是没有歧 义的，因此我们使用(16)得到的经纬度结论，即：

直杆的高度为 2.008 米，位置在（39°31’42” ，79°44’55”E），由 求 得测量日期为 1 月 20 日或 11 月 23 日。

B.求解附件三

处理附件三的数据后(L 与 见附录七)，利用 5.2.1 的方法，由 L 与 两个自由度分别计算经纬度与赤纬角，得到的两个结论差别较大。因此我

(15) )

(16) )

们将（L， ）视为关于时间 t 的向量值函数， 、

H、

、 为待求的参 数，再次利用 lsqcurvefit 函数进行最小二乘法拟合，求出参数：

H=2.887 米

＝31°42′33′′

＝110°22 30

＝− 18°30’37”

Y’轴与正北的夹角是9.9°

结论:直杆高度为 2.887 米，位置在(31°42′33′′ ，110°22 30 )，

由 求得测量日期为 1 月 27 日或 11 月 16 日。

5.3 用日影定时定位模型分析视频(问题四)

在对问题四中的视频进行分析时，我们的目的是利用日影变化求出视 频拍摄的地点与时间。首先对视频进行预处理，使得更易读出日影数据。

先放大视频中的像素点，找到直杆的最高点并用黑色十字进行标记，

再找到直杆底点也用黑色十字标记（见附录八图 1）以整个画面的左上角 的像素点为(1，1)，横向向右每增加一个像素点，y 坐标加 1，纵向向下 每增加一个像素点，x 坐标加 1，则由此读出直杆顶点与底点的坐标值。

再从邻近直杆底部又不会成为阴影顶点的某处向阴影顶端右侧截取小 区域，读出此小区域左上角像素点在整个图中的坐标。将画面由彩色处理 为黑白，用 sobel 函数求出小区域内的阴影边界，阴影边界用白色像素点 表示，其它则是黑色(见附录八图 2)在截取的小区域中按照如上建立新的 坐标系。

最后在小区域内从右向左找出阴影顶点，为了排除噪点的影响，规定 如果向左连续 18 个点都是白色像素点，则此点为阴影顶点。在小区域中 读出不同时间帧上的阴影顶点坐标。分别从第 5，55，500，250，125 帧 开始，每隔 1 分钟找到阴影顶点并读出它在小区域中坐标值（坐标数据见 附录九）

由顶点在小区域中的坐标、小区域左上角坐标、直杆顶点坐标、直杆 底点坐标、直杆实际高度(2 米)可以得到的一系列杆高 H 与影长 L 的数

(17) )

据，将数据代入前面建立的日影定时定位模型中，计算出视频拍摄的地点 为（29°49’ ，_{112°41′ ）}

六．模型检验与误差分析

6.1 模型检验

在我们的模型中最关键的两个数学量就是赤纬角与太阳高度角，由于 很多假设都会对这两个量产生直接影响，它们在计算过程中也有着关键的 作用，下面我们通过验算这两个量来对模型进行检验。

我们从《2015 天文年历》中抽取每个月的月中与月末两天的赤纬角 作为实际测量值，然后用我们的模型计算出这 24 天的赤纬角，并进行对 比，结果见下图。

可以看出除了图中某一特殊点，其它点基本上完全符合，我们又计算 了相对误差，除特殊点外，最大相对误差为 1.31%，大部分均不到 1%，说 明此模型对赤纬角计算较为精确。

我们又找到了一组实地测量的影长数据，来对太阳高度角进行检 验，数据与程序见附录十，结果见下图。

发现也基本符合实际情况 6.2 误差分析

由以上的图可以看出误差不是随机误差，我们做出了相对误差图如下

可以发现相对误差与时间存在着明显的线性关系，我们对相对误差进 行线性拟合，得到关系 error=0.0168t-0.2025，并对模型进行相应修正，即 在原来值基础上乘以因子（1-error），重新计算误差，得到如下结果：

可以发现误差基本上是随机分布，得到较好修正。

在计算赤纬角的时候，只有一个点相差较大，我们可以猜测这一天 是特殊的一天，有可能因为太阳或地球的特殊活动导致这一天赤纬角的异 常，不影响我们模型对赤纬角的精确度。

在计算太阳高度角的时候，我们发现计算值曲线可以看成实际值曲 线向右平移了一小短距离，这可能是由于一回归年为 365.2422 天，并不 是 365 天，而且闰年的时候还要进一步修正，再加上其他星球对地球运动 的影响导致了时差的产生。

七．模型优化

7.1 基于公转椭圆轨道的模型优化

为了消除时差以及太阳公转轨道为椭圆形引起的误差来源，我们查考 了王炳忠研究员的《太阳辐射计算讲座》第一讲的内容，对赤纬角和太阳 高度角的计算进行优化。

太阳赤纬角 = 0.3723 + 23.2567 sin + 0.1149 sin 2 − 0.1712 sin 3 − 0.7580 cos + 0.3656 cos 2 + 0.0201 cos 3

其中 Q=2 ∗ 57.3 ∗ ( − + ∆ )/365.2422，N 为积日， ∆ 为积日 订正值。

我们用模型检验中的数据对新的模型进行优化，发现误差曲线和我们 的模型形状相同，虽然相对误差变小了，但依然存在相对误差与时间线性 相关的问题，需要进一步优化。（具体程序见附录十一）

7.2 基于大气折射率的模型优化

我们发现实际的太阳高度比计算出来的太阳高度要大，这可能是由于 太阳光穿过大气层时产生折射导致的，原理见下图。

根据斯涅尔折射定律，我们可以推导出如下公式，这也解释了为什么 比实际情况偏大。

ℎ = 90° − sin sin(90° − )

八．参考文献

[1]杜春旭,王普,马重芳,吴玉庭,申少青. 一种高精度太阳位置算法.能源工 程,2010,01:41-48.

[2]陈晓勇,郑科科. 对建筑日照计算中太阳赤纬角公式的探讨.浙江建筑,2011,09:6-8 [3]张文华,司德亮,徐淑通,祁东婷. 太阳影子倍率的计算方法及其对光伏阵列布局的 影响.太阳能,2011,09:28-31

[4]谈小生,葛成辉. 太阳角的计算方法及其在遥感中的应用.国土资源遥 感,1995,06:48-57

[5]王炳忠. 第一讲太阳能中天文参数的计算.太阳能,1999,02

[6]中国科学院紫金山天文台. 2015 年中国天文年历. 北京:科学出版社,2014.

九．附录

附录一 北京时间

（h) 10.3 10.4 10.5 10.6 10.7 10.8

时角(°) -29.109 -27.609 -26.109 -24.609 -23.109 -21.609 太阳高度

角(°) 32.396 33.034 33.645 34.228 34.782 35.307 影长(m) 4.728 4.614 4.508 4.410 4.319 4.236

10.9 11.0 11.1 11.2 11.3 11.4 11.5

-20.109 -18.609 -17.109 -15.609 -14.109 -12.609 -11.109 35.801 36.264 36.695 37.093 37.457 37.786 38.081

4.159 4.089 4.026 3.968 3.916 3.869 3.829 北京时间

（h) 11.6 11.7 11.8 11.9 12.0 12.1

时角(°) -9.609 -8.109 -6.609 -5.109 -3.609 -2.109 太阳高度

角(°) 38.339 38.561 38.746 38.894 39.005 39.077 影长(m) 3.793 3.763 3.738 3.719 3.704 3.694

12.2 12.3 12.4 12.5 12.6 12.7 12.8

-0.609 0.891 2.391 3.891 5.391 6.891 8.391 39.112 39.108 39.067 38.987 38.869 38.714 38.522

3.690 3.690 3.696 3.706 3.722 3.743 3.769 北京时间

（h) 12.9 13.0 13.1 13.2 13.3 13.4

时角(°) 9.891 11.391 12.891 14.391 15.891 17.391 太阳高度

角(°) 38.293 38.028 37.727 37.391 37.020 36.616 影长(m) 3.800 3.836 3.878 3.925 3.978 4.037

13.5 13.6 13.7 13.8 13.9 14.0 14.1

18.891 20.391 21.891 23.391 24.891 26.391 27.891 36.179 35.710 35.210 34.680 34.120 33.531 32.916

4.102 4.173 4.251 4.336 4.428 4.527 4.635 北京时间

（h) 14.2 14.3 14.4 14.5 14.6 14.7

时角(°) 29.391 30.891 32.391 33.891 35.391 36.891 太阳高度

角(°) 32.273 31.605 30.911 30.194 29.454 28.692 影长(m) 4.751 4.876 5.010 5.156 5.312 5.481

14.8 14.9 15.0 38.391 39.891 41.391 27.908 27.105 26.281

5.664 5.861 6.075 附录二

function a = declination1( x )

%计算赤纬角,x 为夏至日至现在的公转角度

% 此处显示详细说明

b=23+26/60+21/3600;%黄赤交角 a=asind(cosd(x)*sind(b));

end

function h = elevation1( a,x,y,t )

%已知赤纬角 a,纬度 x,经度 y,北京时间 t,求太阳高度角 h

% 此处显示详细说明

t=15*(timetransform(t,y)-12);%计算时角

h=asind(cosd(a)*cosd(x)*cosd(t)+sind(a)*sind(x));

end

function t = timetransform( t0,x )

%北京时间 t0 转换为当地时间 t,x 为当地经度

% 此处显示详细说明 t=t0+(x-120)/15;

end

H=3;

a=declination1(118.60274);

x=39+54/60+26/3600;%纬度 y=116+23/60+29/3600;%经度 T0=9:0.1:15;

h=elevation1(a,x,y,T0);%太阳高度角 L=H*cotd(h);%影长

plot(T0,L);

附录三:

北京时

间 影长 角度/rad

14:42 1.149626 0.447347 14:45 1.182199 0.439397 14:48 1.215297 0.431701 14:51 1.249051 0.424289 14:54 1.283195 0.417069 14:57 1.317993 0.410112 15:00 1.353364 0.403272 15:03 1.389387 0.396679 15:06 1.426153 0.390298 15:09 1.4634 0.384047 15:12 1.501482 0.377971 15:15 1.540232 0.37197 15:18 1.579853 0.366201 15:21 1.620145 0.360504 15:24 1.661271 0.354973 15:27 1.703291 0.349528 15:30 1.746206 0.344172 15:33 1.790051 0.338958 15:36 1.835014 0.333848 15:39 1.880875 0.328829 15:42 1.927918 0.32385

附录四

Appendix1size.m 文件：

clear;

%从北京时间 12:00 到各数据点记录时刻所经历的小时数 t = [162:3:222]/60;

%该时间段内地球自转过的弧度值 b = t*15*pi/180;

%由附录一中所给坐标计算出的影子长度 p = [1.149625826

1.182198976 1.215296955 1.249051052 1.28319534 1.317993149 1.353364049

1.389387091 1.426152856 1.463399853 1.501481622 1.540231817 1.579853316 1.620144515 1.661270613 1.703290633 1.74620591 1.790050915 1.835014272 1.880875001 1.927918447 ]';

%待确定参数的函数

fun = @(a,x)a(1)*sqrt(1./(cos(a(2))*cos(0.18545).*cos(x+a(3)) + sin(a(2))*sin(0.18545)).^2 - 1) ;

%用 lsqcurvefit 进行拟合来确定参数

[aa res] = lsqcurvefit(fun,[2 0.924 -0.187],b,p);

%各参数值 aa

%拟合误差 res

%计算出拟合值 p1 = fun(aa,b);

%画出拟合值和数据点 plot(b,p,b,p1,'o');

%经度

cita = aa(2)*180/pi

%纬度

beita = 120 + aa(3)*180/pi

appendix1angle.m 文件：

clear;

%从北京时间 12:00 到各数据点记录时刻所经历的小时数 t = [162:3:222]/60;

%该时间段内地球自转过的弧度值 b = t*15*pi/180;

%由附录一中所给坐标计算出的影子与坐标轴ｘ夹角 d = [0.447347454

0.439396868 0.431700895 0.424288556 0.417069084 0.410112214 0.403271514 0.396679059 0.390298113 0.384046702 0.377971042 0.371970396 0.366201319 0.360504459 0.354973256 0.349528408 0.344171895 0.338957945 0.333848406 0.328828803 0.323850082

]';

%待确定参数的函数

fun = @(a,x)atan(sin(a(1)).*cot(x+a(2)) - cos(a(1))*tan(0.18545)./sin(x+a(2))) + a(3) ;

%用 lsqcurvefit 进行拟合来确定参数 [aa res] = lsqcurvefit(fun,[0.5 -0.5 0],b,d);

%各参数值 aa

%拟合误差 res

%计算出拟合值 p1 = fun(aa,b);

%画出拟合值和数据点 plot(b,d,b,p1,'o');

%经度

cita = aa(1)*180/pi

%纬度

beita = 120 + aa(2)*180/pi

%y 轴与正北方向的夹角 C = aa(3)*180/pi()

附录五

北京时间 影长 角度/rad

12:41 1.247256205 -1.431643216 12:44 1.22279459 -1.415610211 12:47 1.198921486 -1.399134289 12:50 1.175428964 -1.382260429 12:53 1.152439573 -1.364908944 12:56 1.12991747 -1.34724127 12:59 1.10783548 -1.328970225 13:02 1.086254206 -1.310276956

13:05 1.065081072 -1.291129615 13:08 1.044446265 -1.271452682 13:11 1.024264126 -1.25120693 13:14 1.004640314 -1.230410479 13:17 0.985490908 -1.209020915 13:20 0.966790494 -1.187120723 13:23 0.948584735 -1.164603359 13:26 0.930927881 -1.141359553 13:29 0.91375175 -1.117573674 13:32 0.897109051 -1.093140261 13:35 0.880973762 -1.067898413 13:38 0.865492259 -1.042048215 13:41 0.850504468 -1.01547483

附录六：

Appendix2size.m 文件：

clear;

t = [41:3:101]/60;

b = t*15*pi/180;

p = [1.247256205 1.22279459 1.198921486 1.175428964 1.152439573 1.12991747 1.10783548 1.086254206 1.065081072 1.044446265 1.024264126 1.004640314

0.985490908 0.966790494 0.948584735 0.930927881 0.91375175 0.897109051 0.880973762 0.865492259 0.850504468 ]';

fun = @(a,x)a(1)*sqrt(1./(cos(a(2))*cos(a(3)).*cos(x+a(4)) + sin(a(2))*sin(a(3))).^2 - 1) ; [aa res] = lsqcurvefit(fun,[2 -0.924 -0.1 -0.187],b,p);

aa res

p1 = fun(aa,b);

plot(b,p,b,p1,'o');

cita = aa(2)*180/pi arfa = aa(3)*180/pi beita = aa(4)*180/pi

appendix2angle.m 文件：

clear;

t = [41:3:101]/60;

b = t*15*pi/180;

d = [-1.431643216 -1.415610211 -1.399134289 -1.382260429 -1.364908944 -1.34724127 -1.328970225

-1.310276956 -1.291129615 -1.271452682 -1.25120693 -1.230410479 -1.209020915 -1.187120723 -1.164603359 -1.141359553 -1.117573674 -1.093140261 -1.067898413 -1.042048215 -1.01547483 ]';

fun = @(a,x)atan(sin(a(1)).*cot(x+a(2)) - cos(a(1))*tan(a(3))./sin(x+a(2))) + a(4) ; [aa res] = lsqcurvefit(fun,[0.689898 -0.702518 0.354657 -0.9],b,d);

aa res

p1 = fun(aa,b);

plot(b,d,b,p1,'o');

cita = aa(1)*180/pi arfa = aa(3)*180/pi beita = aa(2)*180/pi C = aa(4)*180/pi()

附录七：

北京时

间 影长 角度/rad

13:09 3.533142 1.235163 13:12 3.546768 1.219289 13:15 3.561798 1.203475

13:18 3.578101 1.187752 13:21 3.595751 1.17208 13:24 3.614934 1.156493 13:27 3.635426 1.141009 13:30 3.657218 1.125578 13:33 3.680541 1.110271 13:36 3.705168 1.095042 13:39 3.731278 1.079926 13:42 3.758918 1.064912 13:45 3.788088 1.050013 13:48 3.818701 1.035227 13:51 3.85081 1.020543 13:54 3.884585 1.006 13:57 3.919912 0.99156 14:00 3.956876 0.977249 14:03 3.995535 0.963058 14:06 4.035751 0.949003 14:09 4.077863 0.935064

Appendix3all.m 文件：

clear;

t = [69:3:129]/60;

b = t*15*pi/180;

p = [3.533142184 3.546768029 3.561797643 3.578100715 3.595750783 3.61493428 3.635425983 3.657218272 3.680541115 3.705167836 3.731278025 3.758917911 3.788087888

3.818701015 3.850809619 3.88458522 3.919911828 3.956875992 3.99553479 4.035750835 4.077863059 ]';

d = [1.23516337 1.219288884 1.203475346 1.187752148 1.172079623 1.156493028 1.141009042 1.125577726 1.110271467 1.095042216 1.079926208 1.06491228 1.050013044 1.035227146 1.020543181 1.006000029 0.991560277 0.977249323 0.963057754 0.949002921 0.935064071 ]';

y = cat(1,p,d);

fun = @(a,x)[a(1)*sqrt(1./(cos(a(2))*cos(a(3)).*cos(x+a(4)) + sin(a(2))*sin(a(3))).^2 - 1) ; atan(sin(a(2)).*cot(x+a(4)) - cos(a(2))*tan(a(3))./sin(x+a(4))) + a(5) ];

[aa res] = lsqcurvefit(fun,[3 0.504 -0.3554 -0.1735 -0.1788],b,y);

aa res

y1 = fun(aa,b);

plot(b,p,b,y1(1,:),'o');

hold on;

plot(b,d,b,y1(2,:),'o');

hold off;

l = aa(1)

arfa = aa(3)*180/pi cita = aa(2)*180/pi beita = aa(4)*180/pi C = aa(5)*180/pi()

附录八：

图一，二 附录九：

x5 y x56 y x500 y x250 x x125 y

17 679 21 679 21 677 21 679 21 679 17 674 21 675 21 673 21 674 21 674 19 669 19 669 19 666 19 668 19 668 19 664 19 664 23 662 19 664 19 664 23 659 23 659 23 659 23 659 23 659 23 654 23 653 24 651 23 653 23 653 23 649 23 648 21 645 23 647 23 647 21 646 21 645 21 643 21 643 21 643 22 641 21 641 25 639 21 641 21 641 25 634 25 634 25 633 25 634 25 634 27 621 22 627 23 631 23 630 23 630 23 626 23 626 23 624 23 625 23 625 27 622 27 621 27 619 27 620 27 620 27 616 27 616 27 614 27 616 27 616 24 602 27 611 28 609 27 610 27 610 25 605 25 605 29 605 29 604 29 604 29 602 29 601 29 601 29 601 29 601 27 599 27 600 27 599 27 598 27 598 27 596 27 594 27 591 27 593 27 593 27 589 27 589 28 588 28 589 28 589 31 584 27 583 31 581 27 582 27 582 29 580 29 579 29 579 29 580 29 580 29 576 29 577 33 576 29 575 29 575 33 574 33 573 33 571 33 573 33 573 33 574 33 572 33 572 33 572 33 572 33 562 33 562 33 561 33 561 33 561 33 558 33 557 33 557 33 557 33 557 35 552 35 551 35 551 35 552 35 552 31 546 35 547 31 544 31 545 31 545 34 543 35 544 35 543 35 543 35 543 35 540 32 539 35 537 35 538 35 538 35 534 35 535 35 533 35 535 35 535 33 528 33 527 33 526 33 526 33 526 33 524 33 523 33 523 33 523 33 523 33 520 33 519 33 518 33 519 33 519 37 516 34 516 34 515 33 515 33 515 37 513 34 513 37 511 37 511 37 511 37 508 37 508 37 508 37 507 37 507 35 505 35 504 35 503 35 504 35 504 36 500 35 499 36 498 35 499 35 499 39 496 35 495 39 494 35 495 35 495

附录十

实际值 -21.0341 -17.5289 -13.2008 -8.16847 1.956778 3.944083 9.570806 14.60481 计算值 -21.1927 -17.8872 -13.0276 -8.0874 -1.95959 3.906941 9.552383 14.64085 实际值 18.97175 21.83394 23.28017 23.19797 21.44458 18.39186 13.90125 8.831694 计算值 18.64978 21.6903 23.25916 23.2042 21.36501 18.14664 14.0068 8.824438 实际值 3.232861 -2.58783 -8.68275 -13.9257 -18.3383 -21.5473 -23.2839 -23.1317 计算值 3.130255 -2.74069 -8.82444 -14.0068 -18.6498 -21.6903 -23.3069 -23.142

%从天文年历中每月取月中与月末两个时间对两个模型中赤纬角的计算进行检验

%取我国某地实地观测影长对太阳高度角计算进行检验 clear;

M=[31 28 31 30 31 30 31 31 30 31 30 31];

for k=2:2:24,

D(k-1,:)=[2015 k/2 floor(((M(k/2)+1))/2)];

D(k,:)=[2015 k/2 M(k/2)];

end

for k=1:1:24,

d(k)=day([D(k,1),D(k,2),D(k,3)]);

end

A=[-21.03411111 -17.52888889 -13.20080556 -8.168472222 1.956777778 3.944083333 9.570805556 14.60480556 18.97175 21.83394444

23.28016667 23.19797222 21.44458333 18.39186111 13.90125 8.831694444 3.232861111 -2.587833333 -8.68275 -13.92572222 -18.33833333 -21.54727778 -23.28394444 -23.13166667]';

a2=declination2(2015,d);

p2_1=plot(d,A,'o',d,a2,'*');

%hold on;

p2_2=plot(d,a2-A,'o-');

b(1,24)=0;

for k=day([2015 1 1]):1:day([2015 1 31]), b(k)=-180+(k+9)*360/365;

end

for k=day([2015 2 1]):1:day([2015 4 30]), b(k)=-90+(k-day([2015 3 21]))*360/365;

end

for k=day([2015 5 1]):1:day([2015 7 31]), b(k)=0+(k-day([2015 6 22]))*360/365;

end

for k=day([2015 8 1]):1:day([2015 10 31]), b(k)=90+(k-day([2015 9 23]))*360/365;

end

for k=day([2015 11 1]):1:day([2015 12 31]), b(k)=180+(k-day([2015 12 22]))*360/365;

end

for k=1:1:24,

a1(k)=declination1( b(d(k)) );

end hold off;

p1_1=plot(d,A,'diamond',d,a1,'v');

p1_2=plot(d,a1-A,'o-');

附录十一

优化模型的检验 function d = day( date )

%积日的计算 date(1)年 date(2)月 date(3)日 M=[31 28 31 30 31 30 31 31 30 31 30 31];

d=0;

for i=1:1:date(2), d=d+M(i);

end;

d=d+date(3)-M(date(2));

if((mod(date(1),4)==0&&mod(date(1),100)~=0)||(mod(date(1),400)==0));

d=d+1;

end;

end

function a = declination2( N,D )

%精细计算赤纬角

% 此处显示详细说明年 N，累计日 D

N0=79.6764+0.2422*(N-1985)-floor((N-1985)/4);

t=D-N0;

b=2*pi*t/365.2422;

a=0.3723+23.2567*sin(b)+0.1149*sin(2*b)-0.1712*sin(3*b)- 0.758*cos(b)+0.3656*cos(2*b)+0.0201*cos(3*b);

function h = elevation2( N,D,a,x,y,t0 )

%精确计算太阳高度角

% N 年，D 累计日，a 赤纬角，x 纬度，y 经度，t0 北京时间 t=timetransform(t0,y);%北京时换成地方时

N0=79.6764+0.2422*(N-1985)-floor((N-1985)/4);

b=D-N0;

b=2*pi*b/365.2422;

Et=0.0028-1.9875*sin(b)+9.9059*sin(2*b)-7.0924*cos(b)-0.6882*cos(2*b);

t=t+Et/60;%进行时差订正 t=(t-12)*15;

h=asind(cosd(a)*cosd(x)*cosd(t)+sind(a)*sind(x));

end