張志鵬,崔克楠 ,李厚彪

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

太陽(yáng)影子定位的時(shí)空分析

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

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

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

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

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

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

1 影長(zhǎng)關(guān)系表達(dá)式求解

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

圖1 太陽(yáng)照射物體影子示意圖

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

(1)

式中，α為太陽(yáng)高度角，L為直竿長(zhǎng)度，l為陰影長(zhǎng)度。

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

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

(2)

式中,ψ為緯度，δ為太陽(yáng)赤緯角,ω為時(shí)角。

圖2 太陽(yáng)高度角相關(guān)因素球面圖(注：圖中ω0為日出時(shí)角)

查閱資料知太陽(yáng)赤緯角[3](單位：度)，即太陽(yáng)直射點(diǎn)緯度與日角的關(guān)系如下[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部分構(gòu)成：

T=N+ΔN+N0

(5)

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

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

(6)

(7)

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

(8)

(9)

式中，Tc為觀測(cè)時(shí)刻(單位：h)，λ為測(cè)試點(diǎn)經(jīng)度。

文獻(xiàn)[3]中同時(shí)可以查到時(shí)角(單位：rad)的表達(dá)式為：

(10)

(11)

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

7.092 4cosθ-0.682 2cos2θ

(12)

式中，Trs為真太陽(yáng)時(shí)，Eq為平太陽(yáng)時(shí)和真太陽(yáng)時(shí)時(shí)差。

2 公式項(xiàng)的誤差分析

2.1 天安門(mén)廣場(chǎng)的影長(zhǎng)變化

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

圖3 天安門(mén)廣場(chǎng)影長(zhǎng)變化

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

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

利用專業(yè)計(jì)算太陽(yáng)高度角的軟件“云算筆記”將對(duì)應(yīng)時(shí)刻的天安門(mén)影長(zhǎng)計(jì)算出來(lái)，與前文公式計(jì)算結(jié)果對(duì)比如圖4所示。

圖4 陰影長(zhǎng)度對(duì)比

本文定義平均差異比率為每個(gè)時(shí)間點(diǎn)差異的比率均值，即：

(13)

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

η值相對(duì)較小，表明模型較為可信。

3 已知影長(zhǎng)變化反演特定時(shí)空

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

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

表1 某地直桿影長(zhǎng)變化

3.1 已知日期反演空間位置

3.1.1 定值匹配搜索方法

1)確定定值及匹配度

定值匹配可以作為遍歷法的一個(gè)評(píng)判標(biāo)準(zhǔn)，其意為對(duì)于每一個(gè)搜索點(diǎn)建立統(tǒng)一的標(biāo)準(zhǔn)來(lái)判斷該點(diǎn)是否符合要求。以二維空間為例，設(shè)在遍歷法中需要搜尋的點(diǎn)為(xi,yj)，在遍歷的過(guò)程中可以對(duì)每個(gè)點(diǎn)(xi,yj)求出其匹配程度,從中尋找出匹配度最大的點(diǎn)即為目標(biāo)點(diǎn)。

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

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

(14)

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

2)定值匹配方法計(jì)算結(jié)果

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

圖5 定值匹配搜索匹配度

可以看出，在(109°E，19°N) 附近匹配程度最高。考慮到測(cè)量的誤差及精準(zhǔn)度問(wèn)題，取1°的誤差限進(jìn)行細(xì)精度搜索，即初步定位到：

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

3)增大精度小范圍搜索

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

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

如果此時(shí)將誤差限設(shè)為0.1°,則空間復(fù)雜度為400；如果將誤差限設(shè)為0.01°,則空間復(fù)雜度為40 000。相對(duì)于現(xiàn)代最普通的計(jì)算機(jī)來(lái)說(shuō)，這種級(jí)別的空間復(fù)雜度也可以忽略不計(jì)。

4)細(xì)精度定位結(jié)果及誤差分析

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

實(shí)驗(yàn)測(cè)試地點(diǎn)經(jīng)緯度為(109.5°E，18.3°N)，其經(jīng)度誤差為：

(16)

緯度誤差為：

(17)

其誤差相對(duì)來(lái)說(shuō)比較小，認(rèn)為模型較為可信。

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

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

1)最小二乘擬合使用

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

(18)

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

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

l=f(L,ψ,λ)

(19)

2)最小二乘方法結(jié)果

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

3.2 未知日期反演空間位置和時(shí)間

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

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

(20)

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

擬合結(jié)果為如表所示。

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

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

圖6 算法流程圖

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

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

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

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

(20)

圖7 大氣折射太陽(yáng)光路圖

所以實(shí)際太陽(yáng)高度角為：

(21)

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

6 結(jié)束語(yǔ)

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

目前，國(guó)內(nèi)和國(guó)際對(duì)太陽(yáng)影子的精確定位還沒(méi)有很多研究，此類研究對(duì)于衛(wèi)星信號(hào)相對(duì)較差的極端天氣或惡劣地帶的定位有實(shí)際意義。

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

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

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

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

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

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

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

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

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

[10]李小文,高峰,王錦地.單一太陽(yáng)角BRDF數(shù)據(jù)反演過(guò)程中誤差傳播的估計(jì)[J].中國(guó)科學(xué):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

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

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

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

O151.2

A

10.3969/j.issn.1672-4550.2017.01.011