田新宇，苗 楓，劉卓瀾(.中南大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院，湖南 長沙 40083；. 中南大學(xué) 機電工程學(xué)院，湖南 長沙 40083)

2015年“高教社杯”全國大學(xué)生數(shù)學(xué)建模競賽A題“太陽影子定位”要求通過分析提供的影長測量數(shù)據(jù)及視頻中物體太陽影子的變化，確定視頻拍攝的地點和日期.而后，國內(nèi)學(xué)者先后提出了利用影長與位置關(guān)系進行反演定位的模型[1-3]，但是針對定位精度的研究較少，對精度的提高僅針對于模型優(yōu)化的算法研究.例如，于鵬等[4]提出了利用并列選擇遺傳算法提高定位精度,毛廣瑋等[5]利用模擬退火算法提高定位精度.然而，針對模型本身的構(gòu)建誤差研究尚無相關(guān)文獻，而基于正確的模型所做的后續(xù)優(yōu)化算法等研究才是有意義的. 本文基于模型本身的構(gòu)建，進行誤差分析，以減小模型的定位誤差.

在利用視頻影長數(shù)據(jù)進行定位的研究上，何偉梁等[6]利用線性相機模型實現(xiàn)了視頻數(shù)據(jù)定位.但是該論文僅著眼于視頻數(shù)據(jù)的還原，對文中120幀圖片數(shù)據(jù)的獲取方法未有提及.在影長定位的實際應(yīng)用中，手動測量視頻數(shù)據(jù)是極不方便的，本文提出了影長數(shù)據(jù)的自動提取算法以及還原的模型，提高影長定位模型的實用性.

1 影長定位的反演模型

影子長度與角度的變化包含著許多地理位置和時間的信息，這些信息為我們確定出物體所處的位置和時間提供了依據(jù)[7].

在經(jīng)緯度給定時，影子方位角以及高度角的變化公式為：

sinθ=sinφsinδ+cosφcosδcosα

(1)

(2)

式中θ為太陽高度角，ω為太陽方位角，δ為太陽赤緯角，φ為緯度，α為太陽時角.α和δ的計算中包含經(jīng)度和日期值，其詳細計算可見文獻[8].

在桿長H、日期已知的條件下，可以計算出桿的影長.

(3)

在已知影長數(shù)據(jù)及桿長、日期時，可以通過公式(1)、(2)進行反演求取經(jīng)緯度，利用影長誤差平方和最小建立如下的優(yōu)化模型.

min∑(l′-l)2

(4)

s.t. |φ|≤90,|ψ|≤180,θi>0

(5)

式中l(wèi)′為通過公式計算的理論影長，l為實際測量的影長，φ和ψ分別為緯度和經(jīng)度.

該模型是非線性優(yōu)化問題，可以用Matlab或lingo中帶約束條件的優(yōu)化命令求解，也可以用牛頓法、擬牛頓法求解，還可用遺傳算法、粒子群算法等智能算法求解.本文借助R軟件中的全局優(yōu)化GenSA包利用模擬退火算法進行求解[9].若桿長未知，則可以利用相鄰影長變化的相對值代替絕對值，改變優(yōu)化目標(biāo)，即

(6)

若日期未知，將日期做為離散值，分別得到每一天的最優(yōu)解及相應(yīng)的優(yōu)化變量值，然后在得到的365個最優(yōu)值中選取最優(yōu)解，或者增加日期作為連續(xù)的優(yōu)化參數(shù)，若日期求解結(jié)果不為整數(shù)，則分別對其進行向下取整和向上取整，選取最優(yōu)值.

2 影長定位的誤差分析

在反演模型中，影響定位精度的誤差來源主要有兩方面，一是物體影長的測量誤差，二是優(yōu)化模型的模型誤差.前者主要為測量數(shù)據(jù)中的隨機誤差，后者則包括計算精度、優(yōu)化算法的系統(tǒng)誤差以及模型建立不當(dāng)導(dǎo)致的人為誤差.本文中主要研究測量的隨機誤差以及模型建立不當(dāng)對定位精度造成的影響.

在影長定位的誤差分析中，所采用的數(shù)據(jù)皆為仿真數(shù)據(jù)，即利用影長公式結(jié)合研究目的進行系統(tǒng)誤差的設(shè)計得到模擬數(shù)據(jù)，利用模擬數(shù)據(jù)反演經(jīng)緯度與實際經(jīng)緯度進行比較分析，分析的方法主要為統(tǒng)計檢驗以及可視化比較.

2.1 影長測量誤差對定位誤差的影響

影響測量精度的因素有很多，由中心極限定理，測量的隨機誤差服從均值為0的正態(tài)分布，測量誤差的偏倚由隨機誤差的方差反映，在測量絕對誤差限確定時，可由公式(7)、(8)得到隨機誤差項方差的估計.

(7)

(8)

式中d為絕對誤差限，Zα/2為標(biāo)準(zhǔn)正態(tài)分布的α/2分位數(shù).

進而利用不同方差的正態(tài)分布隨機數(shù)構(gòu)造不同測量誤差的仿真數(shù)據(jù)，并通過求解模型(4)式得到對經(jīng)緯度的擬合結(jié)果.在同一地點，不同測量方差對物體影長定位的誤差影響如圖1所示.

圖1 定位誤差與測量誤差關(guān)系圖Fig.1 Quadratic relationship between positioning error and measurement error

從圖中擬合二次曲線的系數(shù)，可以看出經(jīng)度誤差平方隨測量方差增大幅度約為緯度誤差平方的兩倍.測量絕對誤差限增長到10cm時，緯度偏誤將達到0.28度，經(jīng)度偏誤達到0.39度.

2.2 大氣折射對定位結(jié)果的影響

大氣折射會造成影長測量數(shù)據(jù)與理論數(shù)據(jù)間存在定向的偏差.設(shè)太陽的入射角為π/2-θ，大氣的折射率為n，折射角為π/2-θ′，由折射公式可得太陽高度角的余弦值為

(9)

在大氣折射率已知的情況下，可由折射率公式計算出大氣折射下的影長模擬值，進而用模擬值反演經(jīng)緯度，并與未考慮大氣折射的模型進行比較.由于大氣折射率受溫度濕度影響，其數(shù)值并非不變，在大氣折射率未知的情況下，依然可仿真數(shù)據(jù)進行三參數(shù)(經(jīng)緯度、大氣折射率)的擬合與誤差分析. 對一系列位置數(shù)據(jù)，分別用未考慮大氣折射、考慮到大氣折射(大氣折射率未知)以及考慮到大氣折射(大氣折射率已知)的模型進行擬合，并對其結(jié)果進行成對樣本T檢驗，其檢驗結(jié)果見表1.

表1 大氣折射對定位影響的T檢驗結(jié)果
Tab.1 The result of T-test on atmospheric refraction′s impact on positioning

變量均值估計95%置信區(qū)間TP?valueφ1-φ0-0．1654[-0．2511，-0．0797]-4．72470．0032ψ1-ψ0-0．2129[-0．2474，-0．1784]-15．125．2786×10-6φ1-φ2-0．0175[-0．0413，0．0064]-1．79080．1235ψ1-ψ2-0．0247[-0．0566，0．0072]-1．89690．1066

注:φ1考慮大氣折射且折射率已知求解的緯度；φ0未考慮大氣折射求解的緯度；φ2考慮大氣折射且折射率未知求解的緯度；ψ1考慮大氣折射且折射率已知求解的經(jīng)度；ψ0未考慮大氣折射求解的經(jīng)度；ψ2考慮大氣折射且折射率未知求解的經(jīng)度;*95%置信水平下顯著.

由檢驗結(jié)果可以發(fā)現(xiàn)，未考慮大氣折射的模型與實際考慮大氣折射的模型求解結(jié)果在緯度和經(jīng)度上都有顯著差異，而在大氣折射未知時增加參數(shù)進行擬合對求解結(jié)果影響并不顯著.因此，在實際進行反演定位時，應(yīng)當(dāng)考慮大氣折射的存在，確定其值或設(shè)為未知參數(shù)代入影長定位模型.

2.3 地球橢球體模型下的影長定位

影長公式(1)、(2)的建立是基于地球正球體假設(shè)，而實際上地球是個不規(guī)則的橢球體.以地球中心為原點，赤道平面為x-y平面，指向北極的方向為z軸正向，固定某地P緯度為φ，正午時刻該地所在大圓的平面為x-z平面，建立如圖2所示的坐標(biāo)系，則P點在時角為α?xí)r的坐標(biāo)為

R為赤道橢圓面的長軸半徑，e為第一偏心率.則其切平面的法向量為

太陽入射方向的向量為

m=(-cosδ,0,-sinδ)

則橢球體模型中太陽高度角的變化公式為

(10)

圖2 地球橢球坐標(biāo)系Fig.2 Geographic coordinate system with Earth as an ellipsoid

在WGS-84坐標(biāo)系中e2=0.006694379995，代入優(yōu)化模型，利用修正后的橢球模型的模擬數(shù)據(jù)進行擬合，并對原模型結(jié)果與橢球體模型結(jié)果進行T檢驗，檢驗結(jié)果見表2.

表2 地球橢球模型對定位影響的T檢驗結(jié)果
Tab.2 The result of T-test on the impact of Earth shape on positioning

變量均值估計95%置信區(qū)間TP-valueφ1-φ00．1854[0．1780，0．1915]67．19407．308×10-10ψ1-ψ08．418×10-10[-8．3232×10-10，2．5159×10-9]1．23040．2646φ1-φ-0．0043[-0．0175，0．0088]-0．81010．4488ψ1-ψ-0．0018[-0．0109，0．0073]-0．47670．6504

注:φ1橢球模型求解的緯度；φ0未考慮橢球體求解的緯度；φ實際的緯度；ψ1橢球模型求解的經(jīng)度；ψ0:未考慮橢球體求解的經(jīng)度；ψ實際的經(jīng)度;*95%置信水平下顯著.

對橢球體考慮與否的模型求解結(jié)果在緯度擬合上有顯著差異，故地球橢球體模型在進行定位反演時不可忽略.

2.4 誤差分析結(jié)論

通過控制測量誤差可以有效地提高定位的精度，而不可控的大氣折射與地球形狀引起的誤差可通過模型修正規(guī)避.因此利用影長定位，在測量條件控制得當(dāng)?shù)臈l件下，可以實現(xiàn)較高精度，具有實用性.

而利用視頻數(shù)據(jù)還原影長并定位，可以得到足夠大的樣本以減小測量方差，因此可以通過提高數(shù)據(jù)還原精度，控制影長定位的誤差.

3 視頻數(shù)據(jù)的反演模型

3.1 視頻數(shù)據(jù)的提取

在對運動目標(biāo)進行檢測時，幀間差分法是常用的方法，但由于定位模型中使用的數(shù)據(jù)多為有一定時間間隔的離散序列數(shù)據(jù)，因此使用逐幀差分或三幀差分，都會造成巨大的計算量，且在只有序列圖像數(shù)據(jù)時并不可用.因此結(jié)合三幀差分，提出了離散幀的三幀差分.

對于序列圖像1，…，n,取第i(i

首先讀取灰度圖，得到像素矩陣M1,M2,M3；

然后分別用M2,M3對M1做差分，設(shè)定閾值T，進行二值化處理

(11)

接著對差分0-1矩陣做掩膜運算

(12)

同樣的，對于i>n/2，取第i(i>n/2)張圖片為pic1，第i-n/4張圖片為pic2,第i-n/2張圖片為pic3.使用相同的方法進行計算取得圖像矩陣.計算得到的矩陣B，即為影子圖像矩陣.應(yīng)用該方法得到的提取效果如圖3所示.

圖3 影長圖像提取Fig.3 Shadow image extraction

在圖像數(shù)據(jù)時間間隔較小，或者影子變化不顯著時，可以適當(dāng)增大取幀步長，或是減小pic2的閾值.

提取出影長圖像后，可以先進行閉運算再進行開運算以剔除數(shù)值為1的離群噪聲點，然后求取圖像中數(shù)值為1的點到桿底像素點的距離的最大值，即為圖像中的影長.

經(jīng)過實驗，該自動提取算法具有較好的精度，其測量的影子長度與手動測量的結(jié)果相對誤差小于0.25%.

3.2 視頻數(shù)據(jù)的反演

視頻數(shù)據(jù)是經(jīng)過透視變換的，不能僅利用大小固定的比例進行還原，需建立三維空間坐標(biāo)系(圖4)進行透視還原.

圖4 三維透視坐標(biāo)系Fig.4 Three- dimensional perspective coordinate system

以直桿底端L(0,0,0)為原點，地平面為x-y平面，向東方向為x軸正向，向北方向為y軸正向，垂直向上方向為z軸正向建立直角坐標(biāo)系.M(0,0,H)為桿子頂部，O(g,h,i)為攝像機光心,N(0,0,i)為桿上與光心等高的點,P(x,y,z)為影子的端點，P′、M′ 、N′、L′ 等分別為對應(yīng)的成像后的點.

L點坐標(biāo)為(0,0,0)，桿長為H，則M(0,0,H),O點坐標(biāo)已知為(g,h,i)，則N(0,0,i).

(13)

故N′坐標(biāo)為(a,b,i).

(14)

同理可得M′與L′坐標(biāo)分別為(a,b,(1/k+1)i)、(a,b,i+(i-H)/k).又有NO與像平面垂直，可得像平面方程為

(15)

而P、O、P’三點共線，其直線方程為

(16)

聯(lián)立上述兩個方程，可解得P′(x′,y′,z′)為

(17)

(18)

(19)

若P與L處于同一水平面，則z=0(若P處于一傾斜平面，可聯(lián)立PM直線與所在平面方程求解),

(20)

因此像平面中影長為|P′M′|

(21)

將理論值的l′與實際測得的圖上影長構(gòu)造誤差平方和作為優(yōu)化目標(biāo)，建立形同式(4)-(5)的優(yōu)化模型.若k值未知，可以利用式(6)消去k.除了使用計算圖上影長的方法，也可以利用圖像中影子端點坐標(biāo)對實際影長進行還原，構(gòu)建優(yōu)化模型.經(jīng)實驗，在所需透視模型參數(shù)數(shù)據(jù)已知的情況下，視頻數(shù)據(jù)定位具有較高的精度.

[1]馬玉雪,令狐林玉,楊東海,等. 基于二次函數(shù)擬合的太陽影子定位技術(shù)研究[J]. 運城學(xué)院學(xué)報,2015(06):24-28.

[2]於煒力,竇如鳳. 基于數(shù)學(xué)建模的太陽影子定位技術(shù)的研究[J]. 通訊世界,2016(18):270-271.

[3]胡毅華,楊旭龍,劉媛萍. 太陽影子定位模型的構(gòu)建[J]. 洛陽師范學(xué)院學(xué)報,2015(11):13-18.

[4]于鵬,劉澤鋒,郭改慧,等. 基于并列選擇遺傳算法的太陽影子定位方法[J]. 陜西科技大學(xué)學(xué)報(自然科學(xué)版),2017(1):193-197.

[5]毛廣瑋,袁洋,吳涵旭,文學(xué). 基于模擬退火算法利用太陽影長估計經(jīng)緯度的方法[J]. 科技展望,2016(15):139.

[6]何偉梁,凡皓,常欽皓. 基于極值優(yōu)化的日影定位技術(shù)研究[J]. 南通職業(yè)大學(xué)學(xué)報,2016(4):63-66.

[7]蔡志杰. 太陽影子定位[J]. 數(shù)學(xué)建模及其應(yīng)用,2015(4):25-33.

[8]于賀軍. 氣象用太陽赤緯和時差計算方法研究[J]. 氣象水文海洋儀器,2006(3):50-53.

[9] MULLEN K M. Continuous global optimization in R [J]. Journal of Statistical Software, 2014, 60(6): 1-45.