張志鵬,崔克楠 ,李厚彪

(電子科技大學 a.光電信息學院； b.數(shù)學科學學院；四川 成都 611731)

太陽影子定位的時空分析

張志鵬a,崔克楠a,李厚彪b

(電子科技大學 a.光電信息學院； b.數(shù)學科學學院；四川 成都 611731)

該文研究了利用太陽影子進行定位的正演和反演問題。基于太陽下直桿影子長度的變化，利用球面三角函數(shù)，建立了影長與日期、時間、經(jīng)度、緯度、太陽高度角的關系模型，并對此模型進行反演，得到了已知影長變化求解時空定量的方法;并采用定值匹配搜索與最小二乘擬合兩種不同算法進行了對比檢驗,并對此數(shù)學模型進行了基于大氣折射的誤差分析及修正。

最小二乘擬合；匹配搜索；太陽影子定位；太陽高度角

太陽影子定位是一種基于球面三角函數(shù)[1]分析得到陽光下物體影子長度變化規(guī)律與其影響因素，即經(jīng)緯度、日期、太陽高度角[2]的關系，并利用這些關系進行反演獲得測試點的時空位置的一種方法。在軍事行動、文物探測、森林探險等領域，太陽影子定位具有現(xiàn)實意義。

在利用搜索算法時，其搜索精度是要著重考慮的一項內(nèi)容，本文針對搜索結果做出相應的靈敏度和誤差分析。相對于常用的GPS和北斗定位來說，由于其測量都是在地面上進行，所以不可避免地受到天氣、大氣折射的影響。衛(wèi)星定位是在太空中進行，其載體是電磁波而不是自然光，相對來說其受大氣折射的影響可以忽略不計。對于利用自然光為載體的太陽定位，本文將結合光學知識提出關于大氣折射的修正意見。

1 影長關系表達式求解

首先，在不考慮太陽折射的情況下，太陽直射時的影子成像原理如圖1所示。

圖1 太陽照射物體影子示意圖

利用三角函數(shù)公式，從圖1可以直接推出：

(1)

式中，α為太陽高度角，L為直竿長度，l為陰影長度。

另外，根據(jù)如圖2所示的太陽高度角相關因素球面函數(shù)圖，利用球面三角函數(shù)公式，可以推導出太陽高度角與經(jīng)緯度、時角的關系表達式為：

sinα=sinψsinδ+cosωcosψcosδ

(2)

式中,ψ為緯度，δ為太陽赤緯角,ω為時角。

圖2 太陽高度角相關因素球面圖(注：圖中ω0為日出時角)

查閱資料知太陽赤緯角[3](單位：度)，即太陽直射點緯度與日角的關系如下[4]：

δ=0.372 3+23.256 7sinθ+0.114 9sin2θ-

0.171 2sin3θ-0.758cosθ+

0.365 6cos2θ+0.020 1cos3θ

(3)

式中,θ為日角(單位：度),表示為：

(4)

式中，T為日期修正，由3部分構成：

T=N+ΔN+N0

(5)

式中：N為積日，其值等于測試影長的日期到1月1日的時間，平年最后一天積日為365，閏年最后一天積日為366；ΔN表示積日修正日。

N0=79.676 4+0.242 2(xNF-1985)-
floor((xNF-1985)/4),

(6)

(7)

式中：xNF為年份，在本文中取2015年，此時N0=79.942 4；floor為Matlab中命令，表示向下取整；Ld為經(jīng)度差修正值；W為時間差修正值，計算公式為：

(8)

(9)

式中，Tc為觀測時刻(單位：h)，λ為測試點經(jīng)度。

文獻[3]中同時可以查到時角(單位：rad)的表達式為：

(10)

(11)

Eq=0.002 8-1.985 7sinθ+9.905 9sin2θ-

7.092 4cosθ-0.682 2cos2θ

(12)

式中，Trs為真太陽時，Eq為平太陽時和真太陽時時差。

2 公式項的誤差分析

2.1 天安門廣場的影長變化

在已知北京天安門廣場經(jīng)緯度和日期的情況下，利用前文公式可以計算出對應的影長變化數(shù)據(jù)。北京天安門廣場經(jīng)緯度為：(116°23′29″E,39°54′26″N)，計算出在2015年10月22日北京時間9:00-15:00的影長變化數(shù)據(jù)并繪制如圖3所示(測試桿長為3m)。

圖3 天安門廣場影長變化

從圖可看出，在北京時間12∶15左右影長最短，即太陽高度角最小，北京達到正午時刻。北京的正午時刻不在12∶00的原因是北京時間是東八區(qū)中央經(jīng)線即120°E，而北京的經(jīng)度為116°23′，時差約為14min左右。

2.2 與專業(yè)軟件對比的誤差分析

利用專業(yè)計算太陽高度角的軟件“云算筆記”將對應時刻的天安門影長計算出來，與前文公式計算結果對比如圖4所示。

圖4 陰影長度對比

本文定義平均差異比率為每個時間點差異的比率均值，即：

(13)

式中，η為模型計算得到的影長， 為專業(yè)軟件計算得到的影長，N為影長總個數(shù)。計算結果為：

η值相對較小，表明模型較為可信。

3 已知影長變化反演特定時空

在實際的應用中，通常是在測出一段時間內(nèi)影子長度變化之后，利用這些數(shù)據(jù)求出測試點所在的經(jīng)緯度，在某些特殊情況下，也會需要求出測試的日期。下文將會分成已知日期和未知日期兩種情況反演時空。

某地在2015年4月18日測量的一段時間內(nèi)影長變化數(shù)據(jù)如表1所示[5]。

表1 某地直桿影長變化

3.1 已知日期反演空間位置

3.1.1 定值匹配搜索方法

1)確定定值及匹配度

定值匹配可以作為遍歷法的一個評判標準，其意為對于每一個搜索點建立統(tǒng)一的標準來判斷該點是否符合要求。以二維空間為例，設在遍歷法中需要搜尋的點為(xi,yj)，在遍歷的過程中可以對每個點(xi,yj)求出其匹配程度,從中尋找出匹配度最大的點即為目標點。

對于本模型來說，測試的直桿長度為定值，所以可以將桿長作為匹配對象，對于每個經(jīng)緯度計算出來不同時間的影長和太陽高度角，進而計算出對應的桿長，時間列下桿長波動最小的就是最優(yōu)解。

將時間列下桿長的波動程度定義為匹配度，對于上文公式組是一個非線性方程組，帶入不同的經(jīng)緯度不會出現(xiàn)極端值誤差，所以可以用方差衡量時間列數(shù)據(jù)的波動程度,即匹配度用方差衡量表示如下：

(14)

式中，ξ為搜索匹配度，li為計算得到的某個經(jīng)緯度時刻i的陰影長度，μ為計算得到的某個經(jīng)緯度所有測量時刻陰影長度的均值。

2)定值匹配方法計算結果

利用式(1)～式(12)，對經(jīng)緯度進行搜索計算出的匹配度如圖5表示。

圖5 定值匹配搜索匹配度

可以看出，在(109°E，19°N) 附近匹配程度最高?？紤]到測量的誤差及精準度問題，取1°的誤差限進行細精度搜索，即初步定位到：

(109°±1°E，19°±1°N)

3)增大精度小范圍搜索

前文已敘述了精度對于定位的重要性，之所以不直接使用細精度進行搜尋的原因是大范圍細精度搜索會導致程序時間復雜度無比龐大。太陽影子定位的前景是將其做到智能化和軟件化(APP)以使在特殊情況下，如登山、野外生存等情況下能夠快速反應，所以直接細精度搜索實際意義不大。

在已經(jīng)得到整數(shù)位級的定位位置后，根據(jù)已經(jīng)設定的1°的誤差限，在此范圍內(nèi)進行細精度搜索。

如果此時將誤差限設為0.1°,則空間復雜度為400；如果將誤差限設為0.01°,則空間復雜度為40 000。相對于現(xiàn)代最普通的計算機來說，這種級別的空間復雜度也可以忽略不計。

4)細精度定位結果及誤差分析

根據(jù)上述的細精度搜索方法，計算出定位結果為(109.6°E，18.7°N)，利用谷歌地圖定位在海南省三亞市附近。

實驗測試地點經(jīng)緯度為(109.5°E，18.3°N)，其經(jīng)度誤差為：

(16)

緯度誤差為：

(17)

其誤差相對來說比較小，認為模型較為可信。

3.1.2 最小二乘擬合方法[6]

式(1)～式(12)實際上是一個非線性方程組，也就是說這實際上是一個非線性優(yōu)化問題。對于已知部分數(shù)據(jù)的非線性方程組來說，通常有兩種解法：1)將已知的數(shù)據(jù)點帶入，組成一個n維的超定方程組(n為已知數(shù)據(jù)點個數(shù))，對于無約束的方程組來說可以用牛頓法解其他參量的最優(yōu)解，但是對于有約束的超定方程組來說該方法不再適用；2)利用最小二乘法根據(jù)已知的數(shù)據(jù)點對其他參量進行擬合，該方法對有無約束條件不作要求。對于已知影長反演空間位置的問題來說，非線性方程組收到經(jīng)緯度范圍和桿長的約束，用最小二乘擬合更為合適。

1)最小二乘擬合使用

非線性擬合是已知輸入向量x和輸出向量y，在已經(jīng)確定擬合函數(shù)的情況下， 得到系數(shù)向量最優(yōu)解的方法。而在進行曲線擬合時，常采用最小二乘法進行擬合優(yōu)化。其求解原理如下：

(18)

當輸入一組實驗觀測數(shù)據(jù)并設置好初始點后，經(jīng)過數(shù)次迭代，就可以得到確定擬合函數(shù)局部最優(yōu)的參數(shù)向量值。

對于本文討論問題，最小二乘擬合的對象是由式(1)～式(12)組成的非線性方程組，由于日期確定時太陽赤緯角和積日及其修正值可以求出，所以模型對象是桿長與經(jīng)緯度、影長的關系式：

l=f(L,ψ,λ)

(19)

2)最小二乘方法結果

表2 最小二乘法結果(已知日期)

3.2 未知日期反演空間位置和時間

對于未知日期的情況，相對于已知日期的情況下多了一個時間決策變量，由式(5)～式(9)可知，時間決策變量可以視為積日N，所以此時影長的關系式變?yōu)椋?/p>

l=f(L,ψ,λ，N)

(20)

如果此時采用搜索方法，其算法空間復雜度最小為2.365 2×107，在如此大的空間復雜度下經(jīng)緯度尚只能精確到1° ，所以此時采用搜索方式不合適。采用最小二乘法擬合時相對于已知日期時只需多增加一個擬合對象N。

擬合結果為如表所示。

表3 最小二乘法結果(未知日期)

4 定值搜索和最小二乘擬合算法步驟

圖6 算法流程圖

5 大氣折射產(chǎn)生誤差的修正意見

地球表面向太空方向大氣層被分為了對流層、平流層、中間層、暖層和稀薄的散逸層[7]，不同層的大氣折射率不同，這也就導致了太陽光的折射，前文建立的模型并沒有考慮折射帶來的影響。太陽光的折射示意圖如圖7所示(稀薄的散逸層忽略)。

由斯涅爾定理(折射定理)可得[8]：

n0sinα=n1sinα1=n2sinα2=
n3sinα3=n4sinα4

(20)

圖7 大氣折射太陽光路圖

所以實際太陽高度角為：

(21)

其中，n0～n4分別是太空、暖層、中間層、平流層和對流層的折射率。

6 結束語

本中利用計算機算法實現(xiàn)利用太陽影子進行定位，整體過程中盡量避免了主觀性因素，使得定位結果相對準確。定位過程采用了匹配思想和擬合思想，兩種方式適合不同需要情況，決策變量較少時選用搜索算法較為準確，決策變量較多時選用最小二乘擬合方法更契合實際。同時本文提出了定位過程中大氣折射影響的修正意見，在已知實時對流層和暖層折射率時可使結果更準確。

目前，國內(nèi)和國際對太陽影子的精確定位還沒有很多研究，此類研究對于衛(wèi)星信號相對較差的極端天氣或惡劣地帶的定位有實際意義。

[1]王昌明.可照時數(shù)和太陽高度角計算公式的簡化證明[J].山東氣象, 1989(2): 46-48.

[2]王國安, 米鴻濤, 鄧天宏,等.太陽高度角和日出日落時刻太陽方位角一年變化范圍的計算[J].氣象與環(huán)境科學, 2007(9):161-164.

[3]張闖, 呂東輝, 頊超靜.太陽實時位置計算及在圖像光照方向中的應用[J].電子測量技術,2010(11):87-89.

[4]全國大學生數(shù)學建模組委會.2015高教社杯全國大學生數(shù)學建模競賽A題—太陽影子定位[EB/OL].(2015-09-12)[2016-10-20].http://www.mcm.edu.cn/html_cn/block/8579f5fce999cdc896f78bca5d4f8237.html.

[5]龍輝平, 習勝豐, 侯新華.實驗數(shù)據(jù)的最小二乘擬合算法與分析[J].計算技術與自動化,2008(27):20-23.

[6]李延興, 胡新康, 趙承坤,等.光電信號在大氣層傳播中由于折射產(chǎn)生的路徑彎曲[J].中國科學：D輯, 1998(22):149-152.

[7]陳琳.斯涅爾定律數(shù)字模型的設計[J].教育技術導刊, 2009(2):9-10.

[8]葛耀林.便攜式實時太陽夾角計算裝置的實現(xiàn)[J].工業(yè)設計, 2011(3):96-97.

[9]談小生, 葛成輝.太陽角的計算方法及其在遙感中的應用[J].國土資源遙感, 1995 (2): 48-57.

[10]李小文,高峰,王錦地.單一太陽角BRDF數(shù)據(jù)反演過程中誤差傳播的估計[J].中國科學:E輯, 2000(30): 6-11.

Temporal and Spatial Analysis of the Sun Shadow Positioning

ZHANG Zhipenga, CUI Kenana, LI Houbiaob

(a.School of Optoelectronic Information; b.School of Mathematics Sciences, University of Electronic Science and Technology of China, Chengdu 611731, China)

The problem of forward and inversion of location methods with sun shadow are studied in this paper.Based on the change of shadow length, by using spherical triangle function, a relational model between shadow length with data, time, longitude, latitude and solar zenith angle is established.And then using this expression to carry out inversion, specific time and testing place can be got by the shadow length change.In the time and space location, two different methods are compared by using the matching search and the least square, respectively.Finally, the error analysis and correction based on the atmospheric refraction are isproposed.

least squares fitting; matching search; sun shadow positioning; solar elevation

2015-10-23；修改日期:2016-10-11

四川省科技支撐計劃項目(M110474012015GZ0002)；電子科技大學教改項目(2015XBJX0041,2016XJYYB037)；高等學校大學數(shù)學教學研究與發(fā)展中心項目(CMC20160403)。

張志鵬(1995-),本科，主要從事數(shù)學建模和微電子器件以及FPGA電路設計的研究。

李厚彪，副教授，lihoubiao 0189@163.com。

O151.2

A

10.3969/j.issn.1672-4550.2017.01.011