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

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

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

1 影長(zhǎng)定位的反演模型

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

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

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

(1)

(2)

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

在桿長(zhǎng)H、日期已知的條件下，可以計(jì)算出桿的影長(zhǎng).

(3)

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

min∑(l′-l)2

(4)

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

(5)

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

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

(6)

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

2 影長(zhǎng)定位的誤差分析

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

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

2.1 影長(zhǎng)測(cè)量誤差對(duì)定位誤差的影響

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

(7)

(8)

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

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

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

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

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

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

(9)

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

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

變量均值估計(jì)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%置信水平下顯著.

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

2.3 地球橢球體模型下的影長(zhǎng)定位

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

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

太陽(yáng)入射方向的向量為

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

則橢球體模型中太陽(yáng)高度角的變化公式為

(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ù)進(jìn)行擬合，并對(duì)原模型結(jié)果與橢球體模型結(jié)果進(jìn)行T檢驗(yàn)，檢驗(yàn)結(jié)果見(jiàn)表2.

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

變量均值估計(jì)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未考慮橢球體求解的緯度；φ實(shí)際的緯度；ψ1橢球模型求解的經(jīng)度；ψ0:未考慮橢球體求解的經(jīng)度；ψ實(shí)際的經(jīng)度;*95%置信水平下顯著.

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

2.4 誤差分析結(jié)論

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

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

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

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

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

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

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

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

(11)

接著對(duì)差分0-1矩陣做掩膜運(yùn)算

(12)

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

圖3 影長(zhǎng)圖像提取Fig.3 Shadow image extraction

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

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

經(jīng)過(guò)實(shí)驗(yàn)，該自動(dòng)提取算法具有較好的精度，其測(cè)量的影子長(zhǎng)度與手動(dòng)測(cè)量的結(jié)果相對(duì)誤差小于0.25%.

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

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

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

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

L點(diǎn)坐標(biāo)為(0,0,0)，桿長(zhǎng)為H，則M(0,0,H),O點(diǎn)坐標(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’三點(diǎn)共線，其直線方程為

(16)

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

(17)

(18)

(19)

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

(20)

因此像平面中影長(zhǎng)為|P′M′|

(21)

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

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

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

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

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

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

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

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

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

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