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

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

在利用視頻影長數(shù)據(jù)進行定位的研究上，何偉梁等[6]利用線性相機模型實現(xiàn)了視頻數(shù)據(jù)定位.但是該論文僅著眼于視頻數(shù)據(jù)的還原，對文中120幀圖片數(shù)據(jù)的獲取方法未有提及.在影長定位的實際應用中，手動測量視頻數(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)化目標，即

(6)

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

2 影長定位的誤差分析

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

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

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

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

(7)

(8)

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

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

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

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

2.2 大氣折射對定位結果的影響

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

(9)

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

表1 大氣折射對定位影響的T檢驗結果
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%置信水平下顯著.

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

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

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

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

太陽入射方向的向量為

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

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

(10)

圖2 地球橢球坐標系Fig.2 Geographic coordinate system with Earth as an ellipsoid

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

表2 地球橢球模型對定位影響的T檢驗結果
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%置信水平下顯著.

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

2.4 誤差分析結論

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

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

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

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

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

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

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

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

(11)

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

(12)

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

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

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

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

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

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

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

圖4 三維透視坐標系Fig.4 Three- dimensional perspective coordinate system

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

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

(13)

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

(14)

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

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

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

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

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

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

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

[7]蔡志杰. 太陽影子定位[J]. 數(shù)學建模及其應用,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.