黃亞群，李星宇，任瑩瑩，田 潤(rùn)，張懷雄，蘇 茜

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

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

1 影子長(zhǎng)度變化模型

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

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

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

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

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

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

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

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

圖2 影子長(zhǎng)度時(shí)間變化曲線

圖3 影子軌跡變化曲線

2 太陽(yáng)影子定位模型

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

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

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

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

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

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

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

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

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

圖4 測(cè)量影長(zhǎng)和計(jì)算影長(zhǎng)對(duì)比

2.2 拍攝時(shí)間和地點(diǎn)的優(yōu)化模型

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

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

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

附件2中數(shù)據(jù)得到的擬合函數(shù)為:顯然當(dāng)t＝－b/2a時(shí)，函數(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測(cè)量影長(zhǎng)與擬合影長(zhǎng)對(duì)比圖

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

圖6 附件3測(cè)量影長(zhǎng)與擬合影長(zhǎng)對(duì)比圖

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

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

缺少日期信息的情況下，使用窮舉法以每天為步長(zhǎng)對(duì)一年中所有日期用非線性規(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)行窮舉誤差排序。在已升序排列的誤差中，選取誤差在同一量級(jí)的日期作為候選日期序號(hào)，并取出誤差最小的前6個(gè)，如表1和表2所示。

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

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

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

圖7 附件2中n＝200測(cè)量影長(zhǎng)和規(guī)劃影長(zhǎng)對(duì)比圖

圖8 附件3中n＝307測(cè)量影長(zhǎng)和規(guī)劃影長(zhǎng)對(duì)比圖

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

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

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

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

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

圖9 灰度化處理圖

圖10 二值化處理圖

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

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

其中紅色矩形框?yàn)橐獧z測(cè)的影子區(qū)域。取直桿的腳點(diǎn)為原點(diǎn)，從原點(diǎn)開(kāi)始從左向右進(jìn)行豎直方向的掃描線檢測(cè)，找到影子的頂點(diǎn)。

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

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

由模型計(jì)算的影子長(zhǎng)度和實(shí)際測(cè)量影子長(zhǎng)度高度吻合，模擬效果良好，如圖12所示。

圖12 視頻提取影長(zhǎng)與計(jì)算影長(zhǎng)對(duì)比圖

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

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

[1]全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽組委會(huì).2015年高教社杯全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽[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]百度百科.太陽(yáng)高度角 [EB/OL].[2015－12－15].http://baike.baidu.com/view/86609.htm.

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

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

[6]單黎明.太陽(yáng)跟蹤定位技術(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].北京:國(guó)防工業(yè)出版社，2011.

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

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

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

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