黃亞群，李星宇，任瑩瑩，田 潤，張懷雄，蘇 茜

(1.云南大學(xué) 信息學(xué)院，云南 昆明 650091；2.中國科學(xué)院大學(xué) 計算技術(shù)研究所，北京 石景山 100190；3.中南大學(xué) 信息科學(xué)與工程學(xué)院，湖南 長沙 410083)

太陽影子定位技術(shù)是通過分析視頻中物體的太陽影子變化，確定視頻拍攝的地點和日期的一種方法[1]。為保證根據(jù)物體的太陽影子變化精確地確定拍攝地點和拍攝日期，關(guān)鍵在于分析影子長度關(guān)于各個參數(shù)的變化規(guī)律，建立太陽影子定位的數(shù)學(xué)模型。將直桿在陽光下產(chǎn)生的影子軌跡近似代替太陽運動的軌跡，根據(jù)相對運動的原理，將地球繞太陽自轉(zhuǎn)和公轉(zhuǎn)的運動簡化為地球不動太陽繞地球轉(zhuǎn)動，并將太陽繞地球的運動看作在近圓形的軌道上運行。首先，建立影子長度隨日期、當(dāng)?shù)貥?biāo)準(zhǔn)時間、目標(biāo)所在經(jīng)度和緯度及直桿長度變化的動態(tài)數(shù)學(xué)模型，從而計算并繪制太陽影子長度的變化曲線；其次，考慮影子長度隨當(dāng)?shù)貥?biāo)準(zhǔn)時間(北京時間)的變化曲線，以直桿所在經(jīng)度和緯度及直桿長度為決策變量，以實際影子長和計算影子長之差的平方和最小為目標(biāo)函數(shù)，建立太陽影子定位的非線性規(guī)劃模型，結(jié)合曲線擬合選取初值，對視頻信息進(jìn)行預(yù)處理后提取相關(guān)數(shù)據(jù)，從而得到視頻拍攝地點和日期。

1 影子長度變化模型

以直桿在陽光下產(chǎn)生的影子軌跡近似代替太陽運動的軌跡，并根據(jù)相對運動的原理，將地球繞太陽自轉(zhuǎn)和公轉(zhuǎn)的運動簡化為地球不動，太陽繞地球轉(zhuǎn)動，于是太陽繞地球在近圓形的軌道上運行。以直桿所處地面為參考平面，其底端為原點，以南為x軸，以東為y軸建立直角坐標(biāo)，如圖1所示。

設(shè)當(dāng)?shù)鼐暥葹閣d，太陽赤緯角[2]為δ，時角為Ω，太陽高度角[3]為h，太陽方位角[4]為A，則計算公式[5－6]分別為:

圖1 直桿投影坐標(biāo)系

式中，n為日期序列，即某日期在一年內(nèi)的順序號，如n＝1時表示1月1日。

設(shè)直桿的高度為H，影子長度為L，則影子長度動態(tài)變化模型為:

若直桿影子頂點的坐標(biāo)為P(x0，y0)，則影子頂點的橫縱坐標(biāo)分別為:

于是影子長度L關(guān)于各參數(shù)的變化規(guī)律為:

把文獻(xiàn)[1]的附件1中數(shù)據(jù)帶入式(7)，利用MATLAB編程畫出直桿影子隨時間的變化曲線如圖2所示，影子軌跡曲線如圖3所示。

圖2 影子長度時間變化曲線

圖3 影子軌跡變化曲線

2 太陽影子定位模型

確定了直桿影子長度變化規(guī)律后，根據(jù)直桿在水平地面上的影子頂點坐標(biāo)數(shù)據(jù)，分別從已知和缺失拍攝的日期信息兩種情況，建立非線性規(guī)劃模型確定直桿所處的地點和日期，并給出若干個可能的拍攝地點和日期。

2.1 拍攝地點的優(yōu)化模型

式(7)表明，影子長度L的計算值Lc由直桿所在位置的經(jīng)度jd、緯度wd及直桿高度H共同決定，而第i(i＝1，2，…，21)個時刻影子頂點為(xi，yi)時的影子長度的計算值為1，2，…，21)。

在已知測量的準(zhǔn)確日期、當(dāng)?shù)貥?biāo)準(zhǔn)時間及相應(yīng)時刻的影子頂點坐標(biāo)等數(shù)據(jù)的情況下，建立拍攝地點的優(yōu)化模型，利用影子軌跡確定直桿可能的拍攝位置的經(jīng)緯度。

以直桿所在的經(jīng)度、緯度及直桿的長度作為決策變量，以第i(i＝1，2，…，21)個時刻的測量影長Li和計算影長Lci間的誤差最小為目標(biāo)函數(shù)，將描述太陽位置的參數(shù)作為約束條件，建立太陽影子定位的非線性規(guī)劃模型[7－8]:

式中，δ是太陽赤緯角，e為時差，Tt為真太陽時[9]，T為當(dāng)?shù)貢r間，n為日期序號，h為太陽高度角，Ω為時角，wd為當(dāng)?shù)氐木暥?，H為直桿長度，L為直桿的影子長度。

利用制約函數(shù)法，將求解非線性規(guī)劃的問題轉(zhuǎn)化為一系列無約極值問題來求解，再使用內(nèi)點法，即對企圖從內(nèi)部穿越可行域邊界的點在目標(biāo)函數(shù)中加入相應(yīng)的 “障礙約束”，從而保證迭代一直在可行域內(nèi)部進(jìn)行[5]。

利用文獻(xiàn)[1]的附件1中數(shù)據(jù)，通過Lingo程序迭代求解[10]，得到直桿的地理位置為 (108.644 8°E，19.273 5°N)。

將測量影長Li和計算影長Lci進(jìn)行對比，如圖4所示。由圖4可以看出，文獻(xiàn)[1]的附件1中給定時刻的測量影長Li和該時刻對應(yīng)的計算影長Lci幾乎完全一致，因此使用非線性規(guī)劃模型(8)能夠較為精確地計算出直桿的拍攝地點。

圖4 測量影長和計算影長對比

2.2 拍攝時間和地點的優(yōu)化模型

當(dāng)缺失拍攝時間信息時，赤緯角等描述太陽位置的參數(shù)無法確定，因此采用窮舉法添加日期信息，然后再建立優(yōu)化模型求出全局最優(yōu)解。由于非線性規(guī)劃的初始值選取非常重要，所以首先采用曲線擬合的方法粗略地計算出每個日期的初始迭代值，然后對每個日期利用模型(8)進(jìn)行非線性規(guī)劃，并記錄下每個日期的目標(biāo)函數(shù)值，用該日期計算的位置信息利用式(5)求出相應(yīng)時刻的影子頂點坐標(biāo)，最后通過最小二乘法分析出誤差最小的日期，其相應(yīng)的地點即為直桿可能所在地點。

1)非線性規(guī)劃初值選取——擬合粗略計算經(jīng)緯度。

利用文獻(xiàn)[1]的附件2和附件3中的時間－影子長度數(shù)據(jù)，用二次函數(shù)進(jìn)行擬合，得到影子最短時的長度Lemin和此時所對應(yīng)的當(dāng)?shù)貥?biāo)準(zhǔn)時間TS12。

附件2中數(shù)據(jù)得到的擬合函數(shù)為:顯然當(dāng)t＝－b/2a時，函數(shù)有最小值，即附件2中TS12＝15.207 5 h，Lemin＝0.622 2 m。附件3中數(shù)據(jù)得到的擬合函數(shù)為:于是TS12＝12.735 4 h，Lemin＝3.483 3 m，擬合結(jié)果如圖5和圖6所示。

圖5 附件2測量影長與擬合影長對比圖

根據(jù)真太陽時和太陽高度角公式，由當(dāng)?shù)貥?biāo)準(zhǔn)時間TS12、當(dāng)?shù)貢r間T(T＝12∶00)及直桿高度H、赤緯角δ、經(jīng)度jd和緯度wd計算公式分別為:

圖6 附件3測量影長與擬合影長對比圖

綜上所述，根據(jù)式(12)和式(13)選取非線性規(guī)劃中經(jīng)緯度的初始值，能夠得到比較優(yōu)化的解。

2)基于窮舉法求解非線性規(guī)劃模型。

缺少日期信息的情況下，使用窮舉法以每天為步長對一年中所有日期用非線性規(guī)劃模型式(4)進(jìn)行求解，記錄下全局最優(yōu)解，即該日期的經(jīng)緯度的值及所有的誤差，也就是目標(biāo)函數(shù)的值。依據(jù)誤差的大小將誤差和日期按升序排列，便得到直桿位置的經(jīng)緯度。

3)利用最小二乘法誤差分析確定可能的日期。

將文獻(xiàn)[1]的附件2及附件3中的數(shù)據(jù)利用非線性規(guī)劃模型式(8)進(jìn)行求解，進(jìn)行窮舉誤差排序。在已升序排列的誤差中，選取誤差在同一量級的日期作為候選日期序號，并取出誤差最小的前6個，如表1和表2所示。

表1 附件2計算出的前6個候選日期序號及相應(yīng)經(jīng)緯度

表2 附件3計算出的前6個候選日期序號及相應(yīng)經(jīng)緯度

利用非線性規(guī)劃模型式(8)計算的影長和測量影長如圖7和圖8所示。

圖7 附件2中n＝200測量影長和規(guī)劃影長對比圖

圖8 附件3中n＝307測量影長和規(guī)劃影長對比圖

由表1和表2可知，附件2中直桿的可能地點為 (78.907 6°E，39.871 9°N)，日期為5月24日或 (81.239 2°E，39.939 2°N)，日期為 7 月 19日。附件3中直桿的可能地點為 (106.149 5°E，32.844 4°N)，日期為11月3日或 (113.829 5°E，32.904 7°N)，日期為2月7日。

3 監(jiān)控視頻定位模型

在實際應(yīng)用中，需要根據(jù)直桿在太陽下的影子變化的視頻，確定視頻的拍攝地點和時間。

拍攝日期和時間可以從視頻中直接得到，所以該問題的關(guān)鍵是如何從視頻中提取直桿和影子的位置及影子頂點坐標(biāo)等信息，從而將該問題轉(zhuǎn)化為第2節(jié)中討論的問題，然后利用已建立的模型進(jìn)行求解，得到視頻拍攝的可能地點。

首先，對視頻進(jìn)行圖像二值化處理。利用MATLAB按固定頻率取各幀圖像，將圖像進(jìn)行灰度、歸一化及二值化處理[11]，過濾掉不需要的背景部分，使包含影子的目標(biāo)區(qū)域更為明顯，如圖9和圖10所示。

圖9 灰度化處理圖

圖10 二值化處理圖

其次，通過邊緣檢測方法[12]提取影子。為處理方便，裁剪掉上方復(fù)雜部分，保留下方1/3部分，如圖11所示。

圖11 影子區(qū)域水平線檢測示意圖

其中紅色矩形框為要檢測的影子區(qū)域。取直桿的腳點為原點，從原點開始從左向右進(jìn)行豎直方向的掃描線檢測，找到影子的頂點。

再次計算影子長度，得到影子頂點坐標(biāo)。由于所給視頻的影子是經(jīng)過相機的中心投影形成的，而且是從三維的世界坐標(biāo)系映射到了二維的設(shè)備坐標(biāo)系，所以要想準(zhǔn)確計算影子長度與桿的比例，必須進(jìn)行透視處理，做坐標(biāo)轉(zhuǎn)換，求出相機的校準(zhǔn)矩陣[5，11]。

最后，利用第2節(jié)中建立的非線性規(guī)劃模型式(8)，在已知拍攝時間的情況下，得到拍攝地點為(115.248 7°E，40.142 5°N)，在缺失拍攝時間的情況下，得到拍攝地點為 (110.147 2°E，41.687 6°N)。

由模型計算的影子長度和實際測量影子長度高度吻合，模擬效果良好，如圖12所示。

圖12 視頻提取影長與計算影長對比圖

4 結(jié)束語

本文將直桿在陽光下產(chǎn)生的影子軌跡近似代替太陽運動的軌跡，將太陽繞地球的旋轉(zhuǎn)軌道近似為圓形，確定了太陽影子長度隨日期、時間、直桿位置及長度等參數(shù)的變化規(guī)律；然后利用非線性規(guī)劃模型，求出了拍攝時間和拍攝地點。采用曲線擬合的方法確定非線性規(guī)劃的初值，提高了精度，并通過實際數(shù)據(jù)的檢驗，得到較合理的結(jié)果。在模型求解過程中如果利用遺傳算法、粒子群算法等智能算法將會進(jìn)一步減少誤差，這是今后研究的內(nèi)容。

[1]全國大學(xué)生數(shù)學(xué)建模競賽組委會.2015年高教社杯全國大學(xué)生數(shù)學(xué)建模競賽[EB/OL]. [2015－12－15].http://www.mcm.edu.cn/html＿cn/block/c61dfec317 d7a5bd9b2b8efed81c8af3.html.

[2]百度百科.赤緯角[EB/OL].[2015－12－15].http://baike.baidu.com/link？ url＝ xAO9Hmxqeo8I4－1gkHkrc NfobhF2PnLeVa5SkOxldcIFJ＿Nsz－giDraluaghm Gamt E0NSExTbxrL2＿pO162B1K＆qq－pf－to ＝pcqq.c2c.

[3]百度百科.太陽高度角 [EB/OL].[2015－12－15].http://baike.baidu.com/view/86609.htm.

[4]百度百科.太陽方位角 [EB/OL].[2015－12－15].http://baike.baidu.com/view/862818.htm.

[5]武琳.基于太陽陰影軌跡的經(jīng)緯度估計技術(shù)研究[D].天津:天津大學(xué)，2010.

[6]單黎明.太陽跟蹤定位技術(shù)及其應(yīng)用研究[J].空間控制技術(shù)與應(yīng)用，2012，38(3):58－62.

[7]韓中庚.數(shù)學(xué)建模方法及其應(yīng)用[M].北京:高等教育出版社，2009.

[8]司守奎，孫璽菁.數(shù)學(xué)建模算法與應(yīng)用[M].北京:國防工業(yè)出版社，2011.

[9] 維基百科.真太陽時[EB/OL].[2015－12－15].https://en.wikipedia.org/wiki/Solar＿time＃Apparent＿ solar＿time.

[10]汪曉銀，鄒庭榮，周保平.數(shù)學(xué)軟件與數(shù)學(xué)實驗[M].北京:科學(xué)出版社，2010.

[11]RAFAEL C，GONZALE Z，RICHARD E.數(shù)字圖像處理[M].北京:電子工業(yè)出版社，2002.

[12]黃楠.改進(jìn)的數(shù)學(xué)形態(tài)學(xué)圖像邊緣提取算法研究[J].計算機仿真，2012，29(3):288－291.