基于参数拟合的太阳影子定位模型 摘要
太阳影子定位技术就是通过分析视频中物体太阳影子的变化，确定视频拍摄的地 点和日期的一种方法。本文针对四个问题，分别建立了影子长度变化的数学数学模型，
并结合附件数据给出了各问题的直杆所在的可能地点。
针对问题 1，建立了影子长度变化的数学模型，并给出了 2015 年 10 月 22 日北京 时间 9:00-15:00 之间天安门广场一 3 米高的直杆的太阳影子长度的变化曲线，该曲线 呈先减后增的变化趋势。
针对问题 2，首先建立了关于影子长度的函数关系，然后结合附件 1，利用麦夸特 算法，和 1stopt 软件对函数中的未知参数纬度、经度、高度进行了拟合，得到了最优 的参数值为
19 .35,
108. 95 , h
2.03
，将经、纬度进行全球定位，得到最优的可能结 果位于海南省西部。针对问题 3，当建立了影子长度函数后，结合附件 2 和附件 3，利用麦夸特算法，
和 1stopt 软件对函数中的未知参数纬度、经度、高度、测量日期进行了拟合。从而得 到附件 2 对应的最优的参数值
39.92,
79.8, h
2, n
145, nf
2008
，将经、纬度进 行全球定位，得到附件 2 数据所对应的最优可能结果位于新疆境内；另外，得到附件 3 对应的最优的参数值
32.85,
110.18, h
3.07, n
38, nf
2000
，将经、纬度进行全 球定位，得到附件 3 数据所对应的最优可能结果位于湖北省境内。针对问题 4，首先以两分钟为时间间隔，对给出的视频做了画面分割，然后对画面 中直杆的底端、顶端和影子顶点的像素做了提取，并根据像素值得到了 20 个影子长度 值。当测量日期已知时，利用麦夸特算法，和 1stopt 软件对函数中的未知参数纬度、
经度进行了拟合，得到了最后的结果为
41.55,
111.39
，经过全球定位得到对应的 位置为内蒙古境内；当假定测量日期未知时，得到
41.29,
111.25, n
194, nf
2015
， 即积日为 194 天，为 7 月 13 日，恰好是视频拍摄日期。本文建立的模型误差较小，实用性强，有进一步推广的价值。
关键词：
太阳影子；定位；影子长度；参数拟合；麦夸特法一、 问题重述：
太阳影子定位技术是确定视频的拍摄地点和拍摄日期的重要手段，基于此方法，本 文主要解决以下问题：
1.建立影子长度变化的数学模型，并以此模型刻画出 2015 年 10 月 22 日北京时间 9:00-15:00 之间天安门广场（北纬 39 度 54 分 26 秒,东经 116 度 23 分 29 秒）3 米高的 直杆的太阳影子长度的变化曲线。
2.根据固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型，以确定直杆 所处的地点，并将此模型应用于附件 1 中的数据，给出若干可能的地点。
3.根据固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型，以确定直杆 所处的地点和日期，并将此模型应用于附件 2 和附件 3 中的数据，给出若干可能的地点 和日期。
4.在已经估算出直杆高度为 2 米的情况下，根据附件 4 中直杆在太阳下影子的变化，
建立视频拍摄地点的数学模型，并依据此模型给出若干可能的拍摄地点。另外，考虑若 拍摄日期未知，如何根据视频确定拍摄地点与日期。
二、 问题分析：
2.1 问题 1 的分析
建立影子长度变化的数学模型，并给出天安门广场 2015 年 10 月 22 日 9:00-15:00 之间 3 米高的直杆的太阳影子长度的变化曲线。我们知道，影子的长度与直杆的关系如 下图：
所以影子的长度应为直杆高度除以太阳光线与水平地面夹角（太阳高度角）的正切值，
也就是说，本题主要是求太阳高度角随时间的变化关系。而太阳高度角与太阳赤纬角、
观测地地理位置、时角等有关[3]。所以，确定上述几个因素的表达式成为本题的关键。
2.2 问题 2 的分析
问题 2 要求根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型 确定直杆所处的地点，并将模型应用于附件 1 的影子顶点坐标数据。问题 1 已经解决了 影子长度的函数关系，而在坐标系中，影子长度的平方等于影子顶点横坐标与纵坐标的 平方和。所以，把影子顶点坐标带入问题 1 的影长公式即可。接下来要根据附件 1 的数 据确定该直杆所在的位置。把附件 1 的数据代入公式，可以知道当地的经度、纬度和直 杆高度都是未知量，这些量可以利用非线性函数的参数估计来确定。
第一步：确定直杆所处地点的数学模型：
2.3 问题 3 的分析
问题 3 与问题 2 比较相似，但难度更大，因为附件 2 和附件 3 中未给出测量时间，
即本题的未知量为当地的经度、纬度、直杆高度和测量时间。我们可以利用非线性函数 的参数估计来确定这些量。但值得注意的是，问题 3 给出的附件 2 中的数据表明，在北 京时间 12:41 到 13:41 时，观测地的影子长度在减小，这说明此观测地在北京时间正午 时分却处于上午，所以，可以知道，观测地与北京时间存在较长时差，并且它的地方时 要小于北京时间。而附件 3 中给出的数据显示，在北京时间 12:41 到 13:41 时，观测地 的影子长度在增大，说明此观测地的地理经度与北京的地理经度相似或位于北京东部。
2.4 问题 4 的分析
本题要根据高清视频中直杆在太阳下的影子变化情况，定视频拍摄地点的数学模型，
并给出若干个可能的拍摄地点。由上面三个问题的分析可以知道，要确定拍摄地点，关 键要知道影子长度，或者影子顶点坐标。我们可以对视频的画面进行等时间间隔的提取，
然后利用 Matlab 软件提取画面的像素信息和磅值，进而确定直杆底端、顶端和影子顶 点的信息。先利用直杆底端、顶端的信息和直杆高度，可以计算出来单位像素所表示的 真实距离，然后用该距离和影子顶点的信息就可以计算出来每个画面中直杆影子的长度，
然后将其代入问题 3 中的影子长度函数，利用非线性参数拟合就可以确定拍摄地点和日 期。
三、模型假设
1、假定问题求解直杆所处的时间为北京时间。
2、问题四中不考虑摄像机的拍摄位置对数据的影响。
3、杆子最下端与地面相接。
四、符号说明
：太阳高度角； ：
观测地地理纬度； :
观测地地理经度；
：太阳赤纬角；:
T 时角；
:
n 积日（从 1 月 1 日开始计算）；
:
nf
年份；:
h 直杆高度；
:
l
影长；:
x 直杆的太阳影子顶点横坐标；
:
y
直杆的太阳影子顶点横坐标；五、 模型建立与求解
5.1 问题 1 模型的建立与求解：
5.1.1 直杆影子长度模型的建立：
由问题分析可以知道，影子长度与太阳高度角与太阳赤纬角、观测地地理位置、时 角等因素有关。因此在建立模型之前，先做一些准备工作。
1.预备知识
⑴太阳赤纬角
太阳赤纬角表示地球的赤道平面和太阳与地球中心的连线之间的夹角，根据文献[1]， 可以得到太阳赤纬角的计算公式为：
0.3723 23.2567 sin +0.1149sin 2
0.1712sin 3 0.758cos 0.3656 cos 2 0.0202 cos 3 ,
(1)
其中
为日角，其计算表达式[1]为：2 (
0) 365.2422 ,
n N
(2)
式中 n 为积日，表示在以元旦日为起始日的前提下，观测日期在本年内的序号，即 1 月 1 日为起始数 1，其他日期在此基础上累积算得。并且，可以得到N 的公式0 [1]为：
N
0 79.6704 0.2422(
nf
1985)
INT nf ((
1985) / 4),
(3)
式中nf
为当年年份。⑵时角
时角即太阳的位置和当地子午面的偏角，根据文献[2]，可以得到时角的计算公式为：
T 15( t (120 - ) /15 12),
(4)
其中， t 为北京时间，为观测地的经度。⑶太阳高度角
太阳高度角太阳入射光线与地平面的夹角。根据文献[3]，结合太阳赤纬角和时角的计 算公式，可以得到太阳高度角的计算公式为：
sin sin sin cos cos cos( ), T
(5)
2.模型确立：
由图 1 可知，影长等于直杆高度与太阳高度角的正切值之比，即可得到刻画影子长
度变化的数学公式为：
= cot ,
l h
(6)
上式结合公式（6），得到影子长度变化的数学模型为：cot(arcsin(sin sin cos cos cos( ))),
l
h
T
(7)
其中各参数表达式为：0
0
2 ( )
365.2422 ,
79.6704 0.2422( 1985) (( 1985) / 4),
0.3723 23.2567 sin +0.1149sin 2 0.1712sin 3 0.758cos 0.3656 cos 2 0.0202 cos 3 15( (120 - ) /15 12),
n N
N nf INT nf
T t
5.1.2 模型的求解
对于 2015 年 10 月 22 日北京时间 9:00-15:00 之间天安门广场（北纬 39 度 54 分 26 秒,东经 116 度 23 分 29 秒）中 3 米高的直杆，可以得到换算后的各变量值为：
积日：n
295， 年份：nf
2015北京时间：
t
[9,15]
，观测地北京天安门广场的纬度：
39.907
， 观测地北京天安门广场的经度：
116.39
， 直杆高度：h
3。首先，将积日值、年份值分别代入公式（2）和公式（3）中，得到日角为：
2 216.06 365.2422
212.96
，再将θ
212.96
代入公式（1），利用MATLAB
计算得到太阳赤纬角：-11.09
，其次将天安门广场经度
代入公式（4）后，得到时角的表达式为：15* 183.61
T t
，最后将太阳赤纬角
、北京天安门广场纬度值
、直杆高度 h 及时角T
的结果代入公式（7）中，即可得到 2015 年 10 月 22 日北京时间 9:00-15:00 之间天安门广场（北纬 39 度 54 分 26 秒,东经 116 度 23 分 29 秒）中 3 米高的直杆影长随时间
t
变化的数学表达式 为：3cot(arcsin(-0.1234+0.7528*cos(15t-183.61)))
l
在
MATLAB
中运行附录程序，得到直杆影子长度随时间t
的变化曲线图为：图 2 直杆影子长度的变化曲线图
5.2 问题 2 模型的建立与求解：
5.2.1 模型建立：
因为在坐标系中，直杆的影长与横、纵坐标的关系为：
l x
2 y
2(8)
将上式与模型公式（7）结合，得到根据直杆所处地点的模型为：2 2
( , , , , , ) x y cot(arcsin s ( in sin cos cos cos( ) ))
l t n nf h
h
T
(9)
式中各参数表达式为：
0
0
2 ( ) 365.2422 ,
79.6704 0.2422( 1985) (( 1985) / 4),
0.3723 23.2567 sin +0.1149sin 2 0.1712sin 3 0.758cos 0.3656 cos 2 0.0202 cos 3 15( (120 - ) /15 12),
n N
N nf INT nf
T t
，
根据题意可知，与问题一不同的是，在本问中，直杆影子长度可以由直杆影子顶点坐标 表示，最后需要求得直杆所处的地理位置，所以未知参数只有观测地地理纬度
和经度
。5.2.2 模型的求解
求解直杆所处的地点的经、纬度，我们主要利用非线性函数的参数拟合的麦夸特 法，所以下面简单介绍这个算法。
1. 麦夸特法：
麦夸特法也称
LM
算法，主要用于解决非线性最小二乘问题，多用于做曲线拟合。根据文献【4】其主要原理是利用迭代程序，通过计算残差平方和来评估是否达到最佳拟合 效果，当残差平方和最小时迭代结束，得到最优的拟合结果。
其目标函数：
22
1 2
( ) 1
ij i,
i, , 1, 2,
i j
a y f x x a a
N P
式中
2为残差平方和； N 为总实验点数；P
为参数个数。由最小二乘法可知，目标函数应满足误差平方和最小，即：
a
y
min
2对
2展开成二阶泰勒级数，并略去高次项可得：
a D aa a
k
22 1
2
式中
D
为 Hessian 矩阵， a 为带求参数，k
1 , 2 m
。麦夸特法算法通过多次迭代 使本次迭代的参数acur无限接近最佳参数amin，即：
a
min a
curD
1
2 a
cur
5．2.3.模型求解：
针对本问题，利用麦夸特算法，就是要求下面式子的最小值，
22 1 '
( , , )
21 , , , ,
3 i i ,
i
h l l t n nf h
即得到优化问题：
min
2 , , h 。利用
1stOpt
软件对以上问题进行求解。表 1 附件 1 所给直杆影子坐标数据
北 京 时 间（ t ）
x 坐标
y
坐标 影长（ l ） 北 京 时 间（ t ）
x 坐标
y
坐标 影长（l ）14.7 1.0365 0.4973 1.149626 15.25 1.4349 0.5598 1.540232 14.75 1.0699 0.5029 1.182199 15.3 1.4751 0.5657 1.579853 14.8 1.1038 0.5085 1.215297 15.35 1.516 0.5715 1.620145 14.85 1.1383 0.5142 1.249051 15.4 1.5577 0.5774 1.661271 14.9 1.1732 0.5198 1.283195 15.45 1.6003 0.5833 1.703291 14.95 1.2087 0.5255 1.317993 15.5 1.6438 0.5892 1.746206 15 1.2448 0.5311 1.353364 15.55 1.6882 0.5952 1.790051 15.05 1.2815 0.5368 1.389387 15.6 1.7337 0.6013 1.835014 15.1 1.3189 0.5426 1.426153 15.65 1.7801 0.6074 1.880875 15.15 1.3568 0.5483 1.4634 15.7 1.8277 0.6135 1.927918
15.2 1.3955 0.5541 1.501482
北京时间 t 和影子长度 l 的数据如表 1 所示，再由附件 1，可以得到：积日n
108， 年份nf
2015，将以上已知量及公式（9）利用麦夸特算法在1stOpt
软件中实现。（1）按 0.4 的距离给出直杆高度，得到观测地地理纬度
、经度
可能的三组参数值 分别为：第一组：
14.73 114.36, 1.7 h
第二组： 101.
23 73, .5
.69 2 h
⑵.不给定直杆高度，利用
1stOpt
软件全局优化，得到观测地地理纬度
、经度
和直 杆高度 h 的最优化结果为：108.9 19.
5, 2 3
35 .0 h
图 3 基于麦夸特算法的地点定位图
图 4 附件 2 参数拟合图
定位其地理位置，如图 3，参数拟合图形如图 4。
（3）对最优结果进行误差分析：
表 2 影长实测值与计算值列表 影长实测
值
l i ( )
1.1496 25826
1.182 198976
1.2152 96955
1.24905 1052
1.2831 9534
1.31 7993 149
1.353 36404 9
1.38 9387 091
1.42 6152 856
1.46 3399 853 影长计算
值
l i
'( )
1.149 6217
1.182 2144
1.2153 291
1.24898 41
1.2831 987
1.31 7993 4
1.353 3903
1.38 9412 4
1.42 6084 2
1.46 3431 8 影长实测
值
l i ( )
1.5014 81622
1.5402 31817
1.5798 53316
1.62014 4515
1.6612 70613
1.70 3290 633
1.746 20591
1.79 0050 915
1.83 5014 272
1.88 0875 001
1.92 7918 447 影长计算
值
l i
'( )
1.5014 828
1.5402 661
1.5798 126
1.62015 48
1.6612 848
1.70 3313 6
1.746 2463
1.79 0123 5
1.83 4988 6
1.88 0887 5
1.92 7868 9
从
1stOpt
软件计算的结果中，提取出影长的实测值与计算值，如表 2。根据表 2，求出1stOpt
软件计算的影长误差值为' 2
( ( ) ( )) 0.000171
r
l i
l i
，从计算结果可以看到，误差值精确到小数点后 6 位数，即误差非常小。并且，从软件的 计算结果得到相关系数的平方
R
2 0.99999995
，说明麦夸特算法拟合参数的效果非常好，其得到的可能的参数值可信。
5.3 问题 3 模型的建立和求解：
5.3.1 直杆所处地点和时间的模型确定：
相较于问题 2 中求直杆所处的地理位置，问题三要求通过直杆影子的顶点坐标数据，
确定直杆所处的地点和时间，所以，此模型适用于当直杆所处的位置和时间都未知的情 况，模型中包含观测地地理纬度
、经度
、观测年份nf
和观测积日 n 等四个未知参数，那么此问所求模型可以应用模型二的公式（9），得到求解直杆所处地点和时间的模型三 为：
2 2
( , , , , , ) x y cot(arcsin s ( in sin cos cos cos( ) ))
l t n nf h
h
T
(10)
其中，
n nf h , , , ,
均为未知参数。5.3.2.模型求解：
针对本问题，利用麦夸特算法，就是要求下面式子的最小值，
22 '
, , , ( , , , , ) 1
5 , ,
21 i i i
n nf h l l t n nf h
即得到优化问题：
min
2 , , ,n nf h, 。 对附件 2 数据做简单处理得到影长如下表：
表 3 附件 2 所给直杆影子坐标数据
北 京 时 间
（ t ）
x 坐标
y
坐标 影长（ l ） 北 京 时 间（ t ）
x 坐标
y
坐标 影长（l ）12.68 -1.2352 0.173 1.247256 13.23 -0.947 0.3354 1.00464 12.73 -1.2081 0.189 1.222795 13.28 -0.9217 0.3488 0.985491 12.78 -1.1813 0.2048 1.198921 13.33 -0.8965 0.3619 0.96679 12.83 -1.1546 0.2203 1.175429 13.38 -0.8714 0.3748 0.948585 12.88 -1.1281 0.2356 1.15244 13.43 -0.8464 0.3876 0.930928 12.93 -1.1018 0.2505 1.129917 13.48 -0.8215 0.4001 0.913752 12.98 -1.0756 0.2653 1.107835 13.53 -0.7967 0.4124 0.897109
13.03 -1.0496 0.2798 1.086254 13.58 -0.7719 0.4246 0.880974 13.08 -1.0237 0.294 1.065081 13.63 -0.7473 0.4366 0.865492 13.13 -0.998 0.308 1.044446 13.68 -0.7227 0.4484 0.850504 13.18 -0.9724 0.3218 1.024264
利用麦夸特法在
1stOpt
软件中实现：（1） 按 0.2 的距离给出直杆高度，得到观测地地理纬度
、经度
可能的三组参数 值分别为：表 4 优化参数值
参数 第一组 第二组 第三组
地理纬度
36.2 40.35 -41.54地理经度
78.38 79.54 81.82积日n 123 148.34 361
年份
nf
2014 2004 2000直杆高度h 1.8 2 2.2
⑵.不给定直杆高度，利用
1stOpt
软件全局优化，得到观测地地理纬度
、经度
和直 杆高度 h 的最优化结果为：2 79.8 39.92, 145
2008 h
n nf
图 5 附件 2 参数拟合图
图 6 基于
LM
算法的地点定位图 参数拟合图形如图 5，定位其地理位置，如图 6。（3）对最优结果进行误差分析：
表 5 基于附件 2 的影长实测值与计算值列表 影长实测
值
l i ( )
1.247 256
1.2227 95
1.1989 21
1.17542 9
1.1524 4
1.12 9917
1.107 835
1.08 6254
1.06 5081
1.04 4446
影长计算 值
l i
'( )
1.247 261
1.2228 37
1.1988 9
1.17541 7
1.1524 16
1.12 9886
1.107 827
1.08 6239
1.06 5125
1.04 4486
影长实测 值
l i ( )
1.004 64
0.9854 91
0.9667 9
0.94858 5
0.9309 28
0.91 3752
0.897 109
0.88 0974
0.86 5492
0.85 0504
1.00 464
影长计算 值
l i
'( )
1.004 649
0.9854 6
0.9667 67
0.94857 7
0.9308 96
0.91 3736
0.897 105
0.88 1016
0.86 548
0.85 051
1.00 4649
从
1stOpt
软件计算的结果中，提取出影长的实测值与计算值，如表 5。根据表 5，求出1stOpt
软件计算的影长误差值为' 2
( ( ) ( )) 0.00013
r
l i
l i
，从计算结果可以看到，误差值精确到小数点后 5 位数，即误差非常小。并且，从表 2 的 数据可以看出，观测地的直杆影子在北京时间 12:41 到 13:41 之间逐渐减小，说明观测 地的地方时要小于北京时间，观测地地理经度相比北京的地理经度较小。而从图 7 的定 位来看，观测地在新疆境内，而从实际知识可知，新疆的地方时刚好比北京时间晚两小 时，这与分析是相符合的。并且，从软件的计算结果得到相关系数的平方
R
2 0.99999999
， 说明麦夸特法算法拟合参数的效果非常好，其得到的可能的参数值可信。2. 基于附件 3 的模型求解：
表 6 附件 3 所给直杆影子坐标数据
北 京 时 间
（ t ）
x 坐标
y
坐标 影长（ l ） 北 京 时 间（ t ）
x 坐标
y
坐标 影长（l ） 13.15 1.1637 3.336 3.533142 13.7 1.8215 3.2881 3.758918 13.2 1.2212 3.3299 3.546768 13.75 1.8848 3.2859 3.788088 13.25 1.2791 3.3242 3.561798 13.8 1.9488 3.284 3.818701 13.3 1.3373 3.3188 3.578101 13.85 2.0136 3.2824 3.85081 13.35 1.396 3.3137 3.595751 13.9 2.0792 3.2813 3.884585 13.4 1.4552 3.3091 3.614934 13.95 2.1457 3.2805 3.919912 13.45 1.5148 3.3048 3.635426 14 2.2131 3.2801 3.956876 13.5 1.575 3.3007 3.657218 14.05 2.2815 3.2801 3.995535 13.55 1.6357 3.2971 3.680541 14.1 2.3508 3.2804 4.035751 13.6 1.697 3.2937 3.705168 14.15 2.4213 3.2812 4.077863 13.65 1.7589 3.2907 3.731278北京时间 t 和影子长度 l 及坐标值
x 、 y
的数据如表 6 所示，将以上已知量代入公式（10），并利用模麦夸特法在
1stOpt
软件中实现。（1） 按 0.2 的距离给出直杆高度，得到观测地地理纬度
、经度
可能的三组参数 值分别为：表 7 优化参数值
参数 第一组 第二组 第三组
地理纬度
-31 32.64 29.65地理经度
110.58 110.28 109.51积日n 146 35 33
年份
nf
2008 2004 2005直杆高度h 2.7 3 3.3
⑵.不给定直杆高度，利用
1stOpt
软件全局优化，得到观测地地理纬度
、经度
和直 杆高度 h 的最优化结果为：3.07 110.18
32.85, 38
2000 h
n nf
图 7 附件 3 参数拟合图
图 8 基于
LM
算法的地点定位图参数拟合图形如图 7，定位其地理位置，如图 8。
（2） 对最优结果进行误差分析：
表 8 基于附件 3 的影长实测值与计算值列表 影长实测
值
l i ( )
3.533 142
3.5467 68
3.5617 98
3.57810 1
3.5957 51
3.61 4934
3.635 426
3.65 7218
3.68 0541
3.70 5168
影长计算 值
l i
'( )
3.533 202
3.5468 29
3.5617 98
3.57812 3.5958 05
3.61 4864
3.635 312
3.65 7162
3.68 0429
3.70 5131
影长实测 值
l i ( )
3.758 918
3.7880 88
3.8187 01
3.85081 3.8845 85
3.91 9912
3.956 876
3.99 5535
4.03 5751
4.07 7863
3.75 8918
影长计算 值
l i
'( )
3.758 911
3.7880 32
3.8186 68
3.85084 6
3.8845 92
3.91 9935
3.956 905
3.99 5536
4.03 5862
4.07 7923
3.75 8911
从
1stOpt
软件计算的结果中，提取出影长的实测值与计算值，如表 8。根据表 8，求出1stOpt
软件计算的影长误差值为' 2
( ( ) ( )) 0.000262 r
l i
l i
从计算结果可以看到，误差值精确到小数点后 6 位数，即误差非常小。从表 5 的数据可 以看出，观测地的直杆影子在北京时间 13:09 到 14:09 之间逐渐增大，说明观测地的地 方时与北京时间相差不大，其地理经度相比北京的地理经度差异也较小。而从图 11 的 定位来看，观测地在湖北境内，而从实际知识可知，湖北和北京属于同一个时区，这与 分析是相符合的。并且，从软件的计算结果得到相关系数的平方
R
2 0.99999999
，说明 麦夸特算法拟合参数的效果非常好，其得到的可能的参数值可信。5.4 问题 4 模型的建立和求解：
5.4.1 视频画面等时间间隔的提取
我们队附件 4 长达 40 分钟的视频利用
Potplayer
视频处理软件进行每隔两分钟的视 频画面提取，并保证提取的画面大小一样，如下图：图 9‘Appendix4_01.jpg’第一张画面
图 10‘Appendix4_20.jpg’最后一张画面 5.4.2 20 张画面信息的提取
我们首先将 20 张图片进行二值化处理，如第一张画面处理结果为下图：
图 11 Appendix4_01.jpg 二值图像
然后利用
MATLAB
中提取图谱像素的函数impixlinfo
对 20 张图片分别提取直杆底端点、顶点和影子顶点的像素坐标 ；（程序见附件），其中杆底为（ ， ），杆头为a1 a2 （ ， ），a1 a3 影子顶点（b ，b ）。 1 2
影子顶点坐标如下表：
表 9 影子顶点坐标列表 北京时间
( ) t
x 轴坐标 ( )b 1
y 轴坐标
( ) b
2北京时间
( ) t
x 轴坐标 ( )b 1
y 轴坐标
( ) b
28:56 973 770 9:16 880 782 8:58 961 771 9:18 870 782 9:00 951 773 9:20 860 783
9:02 941 775 9:22 851 784 9:04 934 775 9:24 842 785 9:06 924 776 9:26 835 785 9:08 916 776 9:28 822 786 9:10 903 778 9:30 815 787 9:12 895 779 9:32 806 788 9:14 888 781 9:34 798 789
5.4.3 各时间点影长的确定
对于高清图片来讲，其像素的长宽比为 1:1,也就是说，长方向的像素间隔如果和宽 方向的像素间隔相同，则代表的实际距离也是一样的。所以在每个时间点上，我们可以 先 计 算 出 来 直 杆 的 像 素 间 隔 距 离
dx
a
3a
2 和 影 子 的 像 素 间 隔 距 离2 2
1 1 2 2
( ) ( )
dy
b
a
b
a ，则根据直杆的实际长度 2 米就可以计算影子的长度：
l
2* d
y/ d
x （11）式中相应参数的表达式为：
3 2
2 2
1 1 2 2
( ) ( )
dx a a
dy b a b a
。根据公式（11），计算得到 20 个时刻影子的长度如下表：
表 9 20 个时刻的影子长度表 序
号
1 2 3 4 5 6 7 8 9 10
影 长
2.345 6
2.310 1
2.280 5
2.250 9
2.230 3
2.200 8
2.177 2
2.138 7
2.118 1
2.094 4 序
号
11 12 13 14 15 16 17 18 19 20 影
长
2.070 8
2.041 3
2.011 8
1.985 3
1.958 7
1.938
1 1.8997 1.879 1
1.852 5
1.828 9 5.4.3 模型的建立和求解
1．日期已知时的模型建立：
（1）针对本问题，利用麦夸特算法，就是要求下面式子的最小值，
22 1 '
( , )
2 , , , ,
2 ,
0 i i i
t
l l n nf h
即得到优化问题：
min
2 , 。 利用 软件对以上问题进行求解，得到直杆的位置：41.55 111.39
，
将所得经纬度进行全球定位，得到直杆可能的所处位置，位于内蒙古境内，如图 12 所 示：
图 12 问题 4 的地点定位图
（3） 误差分析：
表 10 视频直杆影子实测值与计算值列表 影长实测
值
l i ( )
2.3455 501
2.3101 037
2.2805 107
2.25093 08
2.2302 836
2.20 0750 1
2.177 153
2.13 8742 8
2.11 8068 1
2.10 0328 1 影长计算
值
l i
'( )
2.3483 274
2.3177 012
2.2874 636
2.25778 33
2.2284 876
2.19 9662 8
2.171 2713
2.14 3302 6
2.11 5746 5
2.08 8593 1 影长实测
值
l i ( )
2.0708 154
2.0560 662
2.0118 081
1.98525 29
1.9587 021
1.93 8053 1
1.899 7073
1.87 9065 3
1.85 2528 5
1.82 8946 6
2.07 0815 4 影长计算
值
l i
'( )
2.0618 33
2.0354 571
2.0094 565
1.98382 27
1.9585 476
1.93 3623 3
1.909 0422
1.88 4962 9
1.86 0880 3
1.83 7285 6
2.06 1833
从
1stOpt
软件计算的结果中，提取出影长的实测值与计算值，如表 10。图 13 视频直杆的参数拟合图
从
1stOpt
软件计算的结果中，提取参数拟合的曲线图，如图 13。根据表 10，求出
1stOpt
软件计算的影长误差值为：' 2
( ( ) ( )) 0.034022
r
l i
l i
，从结果可以看到，影长实测值与计算值的误差控制在 0.05 以内，精确度是不错的。并 且，从软件的计算结果得到相关系数的平方
R
2 0.9976
，说明此题应用麦夸特算法拟合 参数的效果较好，其得到的可能的参数值可信。3. 日期未知时的模型建立：
针对本问题，利用麦夸特算法，就是要求下面式子的最小值，
22 '
, , , ( , , , ) 1
4 ,
20 i i ,
i
t n n
n nf l l f h
即得到优化问题：
min
2 , , ,n nf 。利用
1stOpt
软件对以上问题进行求解，得到直杆的位置及测量日期为：41.29 111.25, 194
2015 n
nf
将经、纬度计算的结果进行全球定位，得到直杆可能所处的地点为内蒙古，如图 14。
图 14 视频中直杆的所测位置图
从所得结果看到，积日为 194，换算到日期为 7 月 13 日，恰好是视频拍摄日期，并且，
所得的地理位置与上面的问题结果相同，所以，我们可以说，此模型是合理的。
六、模型评价
本文在正确、清楚地分析题意地基础上，建立了对应问题的各个模型。这些模型都 有效的考虑了全局优化的问题，从而同时得到多个未知参数的最优组合解并且通过软件 实现的参数解精度很高，因此得到的结果可信度较高。不仅如此，本文建立的模型能结 合实际情况对问题进行求解，实用性强，具有很好的推广性。
本文存在的不足是模型较依赖软件的算法实现，并且，第四问中的影子坐标取法并 不十分精确，这对模型求解造成了一定的误差。
参考文献
[1].王炳忠.太阳辐射计算讲座（第一讲）[J]. 太阳能.1999 年第 02 期
[2]. .樊峰鸣，马良涛.单轴 太阳能跟踪系统的 研究[J].河南城建高等 专科学校学 报.2000 年 9 月第 9 卷第 3 期
[3].李玉海.太阳高度角及其计算[J].气象.1977 年第 04 期
[4].程荣兰，陈广志.利用麦夸特理论高效求解Van Genuchten 模型参数[J].2008 年 9 月 9 日
附录 1 问题一：
t=[9:0.1:15];
a=39.907;
n=295;
nf=2015;
r=n-(79.6764+0.2422*(nf-1985)-round((nf-1985)/4));
m=2*180*r/365.2422;b=0.3723+23.2567*sind(m)+0.1149*sind(2*m)-0.1712*sind(3*
m)-0.758*cosd(m)+0.3656*cosd(2*m)+0.02/cosd(3*m);
c=116.39;
l=3./tand(asind(sind(a)*sind(b)+cosd(a)*cosd(b)*cosd((15*t-120+c-12*15))));
plot(t,l)
问题二：
a,b,c;
Variable t,l;
Function
l=c/tan(arcsin(cos(a*pi/180)*cos((t-(120-b)/15-12)*15*pi/180)*cos(10.9813*p i/180)+sin(a*pi/180)*sin(pi/180*10.9813)));
Data;
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
问题三，附件 2
Parameters a[-90,90],b[-180,180],c[0,25],n[0,366],nf[2000,2015];
Variable t,l;
Function
l=c/tan(arcsin(cos(a*pi/180)*cos((t-(120-b)/15-12)*15*pi/180)*cos((0.3723+2 3.2567*sin((pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc(
(trunc(nf)-1985)/4)))/365.2422)))+0.1149*sin(2*(pi/180*(2*180*(trunc(n)-(79 .6704+0.2422*(trunc(nf)-1985)-trunc((trunc(nf)-1985)/4)))/365.2422)))-0.171 2*sin(3*(pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc((tr unc(nf)-1985)/4)))/365.2422)))-0.758*cos((pi/180*(2*180*(trunc(n)-(79.6704+
0.2422*(trunc(nf)-1985)-trunc((trunc(nf)-1985)/4)))/365.2422)))+0.3656*cos(
2*(pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc((trunc(nf )-1985)/4)))/365.2422)))+0.0201*cos(3*(pi/180*(2*180*(trunc(n)-(79.6704+0.2 422*(trunc(nf)-1985)-trunc((trunc(nf)-1985)/4)))/365.2422))))*pi/180)+sin(a
*pi/180)*sin(pi/180*(0.3723+23.2567*sin((pi/180*(2*180*(trunc(n)-(79.6704+0 .2422*(trunc(nf)-1985)-trunc((trunc(nf)-1985)/4)))/365.2422)))
+0.1149*sin(2*(pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-tru nc((trunc(nf)-1985)/4)))/365.2422)))-0.1712*sin(3*(pi/180*(2*180*(trunc(n)- (79.6704+0.2422*(trunc(nf)-1985)-trunc((trunc(nf)-1985)/4)))/365.2422)))-0.
758*cos((pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc((tr unc(nf)-1985)/4)))/365.2422)))+0.3656*cos(2*(pi/180*(2*180*(trunc(n)-(79.67 04+0.2422*(trunc(nf)-1985)-trunc((trunc(nf)-1985)/4)))/365.2422)))+0.0201*c os(3*(pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc((trunc (nf)-1985)/4)))/365.2422)))))));
Data;
12.68 1.247256205 12.73 1.22279459 12.78 1.198921486 12.83 1.175428964 12.88 1.152439573 12.93 1.12991747 12.98 1.10783548 13.03 1.086254206 13.08 1.065081072 13.13 1.044446265 13.18 1.024264126 13.23 1.004640314 13.28 0.985490908 13.33 0.966790494 13.38 0.948584735 13.43 0.930927881 13.48 0.91375175 13.53 0.897109051 13.58 0.880973762
13.63 0.865492259 13.68 0.850504468 问题三 附件 3
a[-90,90],b[-180,180],c[0,25],n[0,366],nf[2000,2015];
Variable t,l;
Function
l=c/tan(arcsin(cos(a*pi/180)*cos((t-(120-b)/15-12)*15*pi/180)*cos((0.3723+2 3.2567*sin((pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc(
(trunc(nf)-1985)/4)))/365.2422)))+0.1149*sin(2*(pi/180*(2*180*(trunc(n)-(79 .6704+0.2422*(trunc(nf)-1985)-trunc((trunc(nf)-1985)/4)))/365.2422)))-0.171 2*sin(3*(pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc((tr unc(nf)-1985)/4)))/365.2422)))-0.758*cos((pi/180*(2*180*(trunc(n)-(79.6704+
0.2422*(trunc(nf)-1985)-trunc((trunc(nf)-1985)/4)))/365.2422)))+0.3656*cos(
2*(pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc((trunc(nf )-1985)/4)))/365.2422)))+0.0201*cos(3*(pi/180*(2*180*(trunc(n)-(79.6704+0.2 422*(trunc(nf)-1985)-trunc((trunc(nf)-1985)/4)))/365.2422))))*pi/180)+sin(a
*pi/180)*sin(pi/180*(0.3723+23.2567*sin((pi/180*(2*180*(trunc(n)-(79.6704+0 .2422*(trunc(nf)-1985)-trunc((trunc(nf)-1985)/4)))/365.2422)))
+0.1149*sin(2*(pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-tru nc((trunc(nf)-1985)/4)))/365.2422)))-0.1712*sin(3*(pi/180*(2*180*(trunc(n)- (79.6704+0.2422*(trunc(nf)-1985)-trunc((trunc(nf)-1985)/4)))/365.2422)))-0.
758*cos((pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc((tr unc(nf)-1985)/4)))/365.2422)))+0.3656*cos(2*(pi/180*(2*180*(trunc(n)-(79.67 04+0.2422*(trunc(nf)-1985)-trunc((trunc(nf)-1985)/4)))/365.2422)))+0.0201*c os(3*(pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc((trunc (nf)-1985)/4)))/365.2422)))))));
Data;
13.15 3.533142184 13.2 3.546768029 13.25 3.561797643 13.3 3.578100715 13.35 3.595750783 13.4 3.61493428 13.45 3.635425983 13.5 3.657218272 13.55 3.680541115 13.6 3.705167836 13.65 3.731278025 13.7 3.758917911 13.75 3.788087888 13.8 3.818701015 13.85 3.850809619
13.9 3.88458522 13.95 3.919911828 14 3.956875992 14.05 3.99553479 14.1 4.035750835 14.15 4.077863059 问题四：
1.图片处理 MATLAB 程序：
I = imread('Appendix4_01.jpg'); %读入图像 if ndims(I)==3%
I=rgb2gray(I);
end %图象灰度化 [~, threshold] = edge(I, 'sobel');
fudgeFactor = .5;
BWs = edge(I,'sobel', threshold * fudgeFactor);
figure, imshow(BWs), title('binary gradient mask');
length(BWs); %二值化矩阵
impixelinfo %图像像素信息
2.已知日期时：
Parameters a[-90,90],b[-180,180];
Variable t,l;
Function
l=2/tan(arcsin(cos(a*pi/180)*cos((t-(120-b)/15-12)*15*pi/180)*cos(21.7583*pi/180)+sin(a*pi/180)*sin(pi /180*21.7583)));
Data;
8.9333 2.345550138 8.9666 2.310103686 9 2.280510732 9.0333 2.250930761 9.06668 2.230283583 9.10003 2.200750111 9.13338 2.177153027 9.16673 2.13874275 9.20008 2.118068051 9.23343 2.100328129 9.26678 2.070815369 9.30013 2.056066243 9.33348 2.011808061 9.36683 1.985252929 9.40018 1.958702065
9.43353 1.938053097 9.46688 1.899707305 9.5 1.879065309 9.53358 1.852528512 9.56693 1.828946617
3.未知日期时
Parameters a[40,45],b[-180,180],n[194,195],nf[2015,2016];
Variable t,l;
Function
l=2/tan(arcsin(cos(a*pi/180)*cos((t-(120-b)/15-12)*15*pi/180)*cos((0.3723+23.2567*sin((pi/180*(2*180
*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc((trunc(nf)-1985)/4)))/365.2422)))+0.1149*sin(2*(pi/1 80*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc((trunc(nf)-1985)/4)))/365.2422)))-0.1712*si n(3*(pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc((trunc(nf)-1985)/4)))/365.2422)))- 0.758*cos((pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc((trunc(nf)-1985)/4)))/365.2 422)))+0.3656*cos(2*(pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc((trunc(nf)-1985)/
4)))/365.2422)))+0.0201*cos(3*(pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc((trunc(
nf)-1985)/4)))/365.2422))))*pi/180)+sin(a*pi/180)*sin(pi/180*(0.3723+23.2567*sin((pi/180*(2*180*(trun c(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc((trunc(nf)-1985)/4)))/365.2422)))
+0.1149*sin(2*(pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc((trunc(nf)-1985)/4)))/3 65.2422)))-0.1712*sin(3*(pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc((trunc(nf)-19 85)/4)))/365.2422)))-0.758*cos((pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-trunc((trunc(
nf)-1985)/4)))/365.2422)))+0.3656*cos(2*(pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(nf)-1985)-tr unc((trunc(nf)-1985)/4)))/365.2422)))+0.0201*cos(3*(pi/180*(2*180*(trunc(n)-(79.6704+0.2422*(trunc(n f)-1985)-trunc((trunc(nf)-1985)/4)))/365.2422)))))));
Data;
8.9333 2.3456 8.9666 2.3101 9 2.2805 9.0333 2.2509 9.0666 2.2303 9.1 2.2008 9.1333 2.1772 9.1666 2.1387 9.2 2.1181 9.2333 2.0944 9.3 2.0413 9.3333 2.0118 9.3666 1.9853 9.4 1.9587 9.4666 1.8997 9.5 1.8791 9.5333 1.8525 9.5666 1.8289
4.20 个二值图片：
