孫慧宇，賀小龍，吳 飛

(1.安徽財(cái)經(jīng)大學(xué)金融學(xué)院，安徽 蚌埠 233030；2.安徽財(cái)經(jīng)大學(xué)管理科學(xué)與工程學(xué)院，安徽 蚌埠 233030； 3.安徽財(cái)經(jīng)大學(xué)統(tǒng)計(jì)與應(yīng)用數(shù)學(xué)學(xué)院，安徽 蚌埠 233030)

基于全局搜索算法的太陽影子定位研究

孫慧宇1，賀小龍2，吳 飛3

(1.安徽財(cái)經(jīng)大學(xué)金融學(xué)院，安徽 蚌埠 233030；2.安徽財(cái)經(jīng)大學(xué)管理科學(xué)與工程學(xué)院，安徽 蚌埠 233030； 3.安徽財(cái)經(jīng)大學(xué)統(tǒng)計(jì)與應(yīng)用數(shù)學(xué)學(xué)院，安徽 蚌埠 233030)

目的 為更好地解決已知某物體影長(zhǎng)、測(cè)量的時(shí)間和日期的情況下確定物體的測(cè)量地點(diǎn)，或已知物體影長(zhǎng)、測(cè)量日期及地點(diǎn)的情況下確定測(cè)量時(shí)間的問題，進(jìn)而解決野外探險(xiǎn)時(shí)地理位置的確定和根據(jù)某物體的視頻或照片等影像資料而確定其拍攝地等實(shí)際問題。方法 針對(duì)太陽影子定位，使用基于窮舉法的全局搜索式算法、最小二乘法的參數(shù)估計(jì)法、比例消參法等方法，分別構(gòu)建全局搜索網(wǎng)絡(luò)模型、參數(shù)方程模型等模型，結(jié)合MATLAB軟件進(jìn)行算法的編寫、圖形的繪制、數(shù)據(jù)的擬合與預(yù)測(cè)等處理。結(jié)果 結(jié)合第一組數(shù)據(jù)，所推算的測(cè)量地位于廣西，結(jié)合第二組數(shù)據(jù)，所推算的測(cè)量地在蒙古國(guó)。結(jié)論 利用全局搜索式算法及參數(shù)方程模型，結(jié)合MATLAB軟件，可以解決太陽影長(zhǎng)、直桿長(zhǎng)度、測(cè)量地點(diǎn)、測(cè)量日期、測(cè)量時(shí)間等5個(gè)變量中知四求一、知三求二等問題。

太陽影子定位；全局網(wǎng)絡(luò)搜索；測(cè)量；MATLAB

0 引 言

物體在燈光或自然光線的照射下會(huì)產(chǎn)生影子并且影子的長(zhǎng)度會(huì)因?yàn)槟承┮蛩匕l(fā)生改變，但是如何確定太陽影子的作用機(jī)制并根據(jù)其中多個(gè)未知量進(jìn)行計(jì)算而求得其他幾個(gè)未知量并沒有深入的研究，而太陽影子定位可廣泛應(yīng)用于野外求生、刑事偵查、軍事偵察以及科研計(jì)算中，因此，對(duì)太陽影子及其影響機(jī)制以及相互遞推關(guān)系進(jìn)行探討顯得尤為重要。

1 太陽影子的影響機(jī)制

為分析不同因素對(duì)于物體影長(zhǎng)的影響機(jī)制并建立相應(yīng)的數(shù)學(xué)模型來反映這種影響機(jī)制，同時(shí)利用所建立的模型繪制2015年10月22日北京時(shí)間9∶00～15∶00北京天安門廣場(chǎng)3 m高的直桿的太陽影子變化曲線,并以此為例。明確太陽高度角α、測(cè)點(diǎn)緯度φ、太陽赤緯σ以及太陽時(shí)角t等參數(shù)，參考三角函數(shù)表，對(duì)太陽影子長(zhǎng)度的變化機(jī)制做出評(píng)價(jià)并繪制相應(yīng)的曲線[1]。

1.1 數(shù)據(jù)處理

(1)太陽赤緯σ：表示太陽光線與地球赤道面的夾角，一年四季每天都在變動(dòng)著，冬至日σ=-23 °27′，春分日和秋分日σ=0 °，夏至日σ=23 °27′。若定義每年的1月1日日期序數(shù)為1，且之后的每一天日期序數(shù)依次遞增，則結(jié)合相關(guān)資料可得太陽赤緯計(jì)算公式：

例如太陽赤緯為-15 °13′，即太陽直射15 °13′S。

(2)太陽時(shí)角t：t為太陽時(shí)角，以當(dāng)?shù)卣鐬? °，上午為負(fù)，每后退1小時(shí)減15 °，下午為正，每前進(jìn)1小時(shí)加15 °，僅以t表示太陽時(shí)。例如北京時(shí)間9∶00～15∶00的太陽時(shí)角為-45～45 °，每分鐘變化0.25 °。

(3)太陽高度角α：太陽高度角即為太陽入射光線和地面切線的夾角(見圖1)，其中α為太陽高度角，I1為太陽入射平行光線，I2地面切線，通過查閱相關(guān)資料[2]可得太陽高度角計(jì)算公式：

sinα=sinσsinφ+cosσcosφcost

(1)

其中α為太陽高度角、φ為測(cè)點(diǎn)緯度、σ為太陽赤緯，t為太陽時(shí)角。

1.2 確定參數(shù)變量關(guān)系

將太陽高度角、桿長(zhǎng)以及影長(zhǎng)的互相影響抽象為幾何問題，可得影長(zhǎng)示意圖(圖2)。[1]

圖1 太陽高度角示意圖 圖2 影長(zhǎng)示意圖

據(jù)此可得桿長(zhǎng)、影長(zhǎng)以及太陽高度角之間的關(guān)系：

(2)

其中l(wèi)為影長(zhǎng)，h為桿長(zhǎng)，α為太陽高度角。

結(jié)合公式(1)和公式(2)可得影響物體影長(zhǎng)與太陽時(shí)角、測(cè)點(diǎn)緯度、太陽赤緯以及桿長(zhǎng)的參數(shù)函數(shù)：

(3)

其中:太陽赤緯、測(cè)地緯度、太陽時(shí)角保持不變時(shí)，則影長(zhǎng)隨桿長(zhǎng)的增加而增加；

桿長(zhǎng)、測(cè)地緯度、太陽時(shí)角保持不變時(shí)，則太陽直射點(diǎn)離測(cè)量地越遠(yuǎn)，影長(zhǎng)越短；

桿長(zhǎng)、太陽赤緯、太陽時(shí)角均保持不變時(shí)，則測(cè)量地離太陽直射點(diǎn)越遠(yuǎn)，影長(zhǎng)越短；

桿長(zhǎng)、太陽赤緯、測(cè)地緯度均保持不變時(shí)，則當(dāng)太陽時(shí)角小于零時(shí)，影長(zhǎng)隨太陽時(shí)角的增大而減小，當(dāng)太陽時(shí)角大于零時(shí)，影長(zhǎng)隨太陽時(shí)角的增大而增大。

圖3 直桿的影長(zhǎng)隨時(shí)間變化圖

以2015年10月22日北京時(shí)間9∶00～15∶00時(shí)北京天安門廣場(chǎng)3m高的直桿為例，利用太陽影子圖像的繪制，結(jié)合影子定點(diǎn)坐標(biāo)，利用MATLAB軟件可繪制其圖像(見圖3)。

2 太陽影子影響因素中單因素未知求解方法

根據(jù)平面上影子的頂點(diǎn)坐標(biāo)求解出影子的長(zhǎng)度，從而針對(duì)時(shí)間以及影長(zhǎng)進(jìn)行擬合[3]，進(jìn)而預(yù)測(cè)北京時(shí)間14∶42時(shí)之前的時(shí)點(diǎn)在觀測(cè)地點(diǎn)直桿的影子長(zhǎng)度，因?yàn)檎鐣r(shí)刻影子是一天中最短的時(shí)刻，所以通過所擬合的函數(shù)頂點(diǎn)來推測(cè)測(cè)量地的正午時(shí)刻所對(duì)應(yīng)的北京時(shí)間，從而通過時(shí)間差來計(jì)算兩地經(jīng)度差，進(jìn)一步確定觀測(cè)地點(diǎn)的經(jīng)度；針對(duì)觀測(cè)點(diǎn)緯度，結(jié)合公式(3)及附件(1)[4]中的數(shù)據(jù)通過計(jì)算不同時(shí)刻的影長(zhǎng)并構(gòu)建比例消元模型以消除桿長(zhǎng)對(duì)計(jì)算的影響，進(jìn)而對(duì)以測(cè)地緯度為未知量一元線性方程進(jìn)行求解，最終得到觀測(cè)點(diǎn)的經(jīng)緯度數(shù)據(jù)。

2.1 數(shù)據(jù)的處理

(1)相關(guān)數(shù)據(jù)的獲取

通過觀察和測(cè)量，得到太陽影子的定點(diǎn)坐標(biāo)(見表1)。

表1 直桿影子定點(diǎn)坐標(biāo)表

(2)觀測(cè)地經(jīng)度的數(shù)據(jù)準(zhǔn)備

圖4 直桿影長(zhǎng)示意圖

根據(jù)文獻(xiàn)[4]中的數(shù)據(jù)利用Excel進(jìn)行二次擬合,可得到北京時(shí)間14∶42～15∶42時(shí)之間的影長(zhǎng)與時(shí)間的函數(shù)關(guān)系：

y=0.0004x2+0.0305x+1.12083

(5)

其中定義北京時(shí)間14∶42時(shí)刻為x=1，此后時(shí)間每增加3min，x的值便增加1；同理，時(shí)間每減少3min，x的值便減少1。通過MATLAB繪制函數(shù)曲線進(jìn)行預(yù)測(cè)可得到直桿影長(zhǎng)示意圖(見圖4)。

利用MATLAB數(shù)據(jù)游標(biāo)工具不難得出該曲線的最低點(diǎn)為(-38,0.538 9)，即北京時(shí)間為12∶45時(shí)影長(zhǎng)最短，為0.538 9m。由影長(zhǎng)變化規(guī)律[5]可知，北京時(shí)間為12∶45時(shí)測(cè)量地的時(shí)間為正午12∶00時(shí)，因此測(cè)量地與北京存在45min的時(shí)差，即可求得測(cè)量地的經(jīng)度約為東經(jīng)109°10′。

2.2 求解方法

根據(jù)表1提供的21組不同時(shí)間對(duì)應(yīng)的影子坐標(biāo)，計(jì)算出相應(yīng)的影子長(zhǎng)度h，再根據(jù)公式(3)，在21組數(shù)據(jù)中任意選取2組數(shù)據(jù)進(jìn)行比例消元，可以得到如下公式：

其中每天求的太陽赤緯[6]是一樣的，并且知道日期，可以根據(jù)太陽赤緯公式y(tǒng)=23.45×sin(360×(284+n)/365)求出每天相應(yīng)的赤緯角度。其時(shí)角也可以根據(jù)不同時(shí)間對(duì)應(yīng)求出，最后模型中只剩下測(cè)量地的緯度是未知的。

2.3 結(jié)果及分析

為了消除桿長(zhǎng)未知對(duì)計(jì)算結(jié)果的影響[7]，選取文獻(xiàn)[4]14∶45時(shí)與15∶15時(shí)2個(gè)時(shí)點(diǎn)為例，從而構(gòu)建比例式以求解測(cè)地緯度。計(jì)算結(jié)果：

(6)

圖5 影長(zhǎng)比例圖

利用MATLAB繪制圖像(見圖5)。

利用MATALB數(shù)據(jù)游標(biāo)工具可知當(dāng)x分別取-1.244、0.018 14、0.337 4時(shí)公式(6)成立，此時(shí)的測(cè)地緯度為南緯71.312 °、南緯1.040 °、北緯19.341 °。通過計(jì)算可得到剩余9組數(shù)據(jù)，通過查閱相關(guān)經(jīng)緯度即可確定觀測(cè)地，經(jīng)整理可得觀測(cè)地時(shí)間、零點(diǎn)x值、影長(zhǎng)比、緯度匯總表(見表2)。

表2 觀測(cè)地對(duì)應(yīng)數(shù)據(jù)匯總表

根據(jù)附件(1)[4]給出的影子坐標(biāo)，結(jié)合相關(guān)數(shù)據(jù)求得可能的觀測(cè)點(diǎn)為廣西、海南、貴州。其中廣西省的頻數(shù)為8，因此推斷最可能的觀測(cè)地點(diǎn)是廣西省。

3 太陽影子影響因素中雙因素未知求解方法

根據(jù)平面上影子的頂點(diǎn)坐標(biāo)求解出影子的長(zhǎng)度，建立基于“最小二乘法”[8]的數(shù)據(jù)擬合模型。根據(jù)時(shí)間關(guān)系對(duì)影長(zhǎng)進(jìn)行擬合，進(jìn)而求出影子最短時(shí)的時(shí)間。再根據(jù)正午時(shí)刻影子最短原理[9]，通過擬合函數(shù)推測(cè)出測(cè)量地的正午時(shí)刻對(duì)應(yīng)的北京時(shí)間，從而根據(jù)時(shí)間差來計(jì)算兩地經(jīng)度差。針對(duì)觀測(cè)點(diǎn)緯度，建立基于“窮舉法”的全局搜索算法模型，利用窮舉法分別對(duì)桿長(zhǎng)、日期、緯度進(jìn)行三層循環(huán)遍歷[10]，利用公式(3)求出相對(duì)應(yīng)的影長(zhǎng)，將求出的影長(zhǎng)于實(shí)際影長(zhǎng)的差的絕對(duì)值進(jìn)行累加，最后比較累加值[11]，最小的累加值相對(duì)應(yīng)的緯度和日期即為所需要求的緯度和日期。

3.1 數(shù)據(jù)處理

(1)經(jīng)度求解的數(shù)據(jù)準(zhǔn)備

通過觀測(cè)可獲得一組太陽影子頂點(diǎn)坐標(biāo)并對(duì)其進(jìn)行整合后計(jì)算影子長(zhǎng)度(見表3)。

表3 觀測(cè)地坐標(biāo)及對(duì)應(yīng)影子長(zhǎng)度匯總表

注:定義x=1對(duì)應(yīng)文獻(xiàn)[4]中的起始觀測(cè)時(shí)刻，時(shí)間每增加3min，x的值增加1。

(2)緯度求解的數(shù)據(jù)準(zhǔn)備

a.桿長(zhǎng)h:假設(shè)桿長(zhǎng)的范圍為0.1～10 m。

b.太陽赤緯y:定義每年的1月1日的日期序數(shù)為1進(jìn)行編號(hào)(忽略閏年的影響)，其最大值為365。通過求太陽赤緯的公式y(tǒng)=23.45×sin(360×284+n)/365)

推算出相應(yīng)的赤緯角度y。

d.太陽時(shí)角θ:定義太陽時(shí)角的計(jì)算公式為θ=15×(t-12)，其中t表示測(cè)量地的測(cè)量時(shí)間與零時(shí)的差值，單位以小時(shí)計(jì)。再結(jié)合文獻(xiàn)[4]中給出的時(shí)間和前述結(jié)論中已求得的經(jīng)度推算出測(cè)量地的即時(shí)時(shí)間，就可以準(zhǔn)確算出太陽時(shí)角[12]，其中上午的太陽時(shí)角為正，下午的太陽時(shí)角為負(fù)。

3.2 求解方法

在經(jīng)度求解方面，可采用本文太陽影子影響因素中單因素未知求解中的方法求解。在緯度求解時(shí)，引入基于窮舉法的全局搜索式模型，將求桿長(zhǎng)所需要的4個(gè)參數(shù)，桿長(zhǎng)、太陽赤緯、測(cè)地緯度和太陽時(shí)角，根據(jù)實(shí)際經(jīng)驗(yàn)和理論上數(shù)據(jù)分別對(duì)其賦值遍歷，求出估計(jì)桿長(zhǎng)。再將求出的估計(jì)桿長(zhǎng)與文獻(xiàn)[4]的實(shí)際桿長(zhǎng)作差，將21組數(shù)據(jù)的差的絕對(duì)值累加。對(duì)每組數(shù)據(jù)的累加結(jié)果進(jìn)行比較，其中誤差最小即累加值最小，對(duì)應(yīng)的桿長(zhǎng)、太陽赤緯、測(cè)地緯度和太陽時(shí)角即符合條件。

3.3 結(jié)果分析

針對(duì)經(jīng)度求解可利用MAYTLAB繪制其擬合圖像并推測(cè)經(jīng)度(見圖6～圖9)。

圖6 影子長(zhǎng)度的散點(diǎn)圖 圖7 影子長(zhǎng)度的擬合曲線圖

圖8 影子長(zhǎng)度的散點(diǎn)圖 圖9 影子長(zhǎng)度的擬合曲線圖

針對(duì)2組數(shù)據(jù)分別得出如下結(jié)論：第一組數(shù)據(jù)中日期序列序號(hào)n=71，即3月12日，桿長(zhǎng)h=1.5 m，緯度為38度,結(jié)合API坐標(biāo)反查與所求得的經(jīng)度可知其測(cè)量地點(diǎn)大約在新疆維吾爾自治區(qū)。

第二組數(shù)據(jù)中日期序列序號(hào)n=294，即10月21日，桿長(zhǎng)h=1.5 m，緯度為46 °，結(jié)合API坐標(biāo)反查與所求的得經(jīng)度可確定其測(cè)量地點(diǎn)大約在蒙古國(guó)。

4 結(jié) 論

在太陽影子影響機(jī)制方面：建立的太陽影子影響機(jī)制模型：

并以此模型為藍(lán)本，通過比例消元、參數(shù)估計(jì)以及基于窮舉法的全局搜索算法等方法，針對(duì)不同的觀測(cè)數(shù)據(jù)可以做到知四求一甚至知三求二。即測(cè)量地、測(cè)量時(shí)間、測(cè)量日期、物體影長(zhǎng)和直桿長(zhǎng)度中3個(gè)或3個(gè)以上為已知時(shí)，剩余變量可以求解得出。

[1]Wu L,Cao X C,Foroosh H.Camera calibration and geo-location estimation from two shadow trajectories[J].Comput Vis Image Underst,2010,114(08):915-927.

[2]Junejo I,Foroosh H.GPS coordinates estimation and camera calibration from solar shadows[J].Comput Vis Image Underst,2010,114(09):991-1003.

[3]屈名，王征兵，王德麾.基于交比不變性的太陽定位算法的研究[J].Silicon Valley,2013,(19):53-55.

[4]中國(guó)工業(yè)與應(yīng)用數(shù)學(xué)學(xué)會(huì)，2015全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽A題[EB/OL].http://www.mcm.edu.cn/html.cn/node/ac8b96613522ef62c019dlcd45a125e3.html.2015-09-14.

[5]林根石.利用太陽視坐標(biāo)的計(jì)算進(jìn)行物高測(cè)量與定位[J].南京林業(yè)大學(xué)學(xué)報(bào),1991,15(03):89-93.

[6]姚宋涵，王殿海，曲昭偉.基于二流理論的擁擠交通當(dāng)量排隊(duì)長(zhǎng)度模型[J].東南大學(xué)學(xué)報(bào)：自然科學(xué)版，2007,37(03)：521-526.

[7]陳佳.周濤.曾祥平.城市異常路段交通擁擠——消散分析[J].中國(guó)水運(yùn)，2007,7(02)：63-64.

[8]吳正.高速公路交通事故和干涉車流波的非線性數(shù)學(xué)模型[J].西安公路交通大學(xué)學(xué)報(bào)，2001,21(02)：77-80.

[9]洪文，馮守平，吳本中.利用LINGO建立最優(yōu)化模型[M].長(zhǎng)春：吉林大學(xué)出版社，2005:4-59.

[10]郭耀煌，鐘小鵬.動(dòng)態(tài)車輛路徑問題排隊(duì)模型分析[J].管理科學(xué)報(bào),2006，9(01)：34-37.

[11]何原標(biāo).用天文年歷進(jìn)行表差和經(jīng)度計(jì)算[J].測(cè)路通報(bào)1983,(06)：11-12.

[12]肖智勇,劉宇翔.一種新的緯度測(cè)量方法[J].大學(xué)物理，2010,29(09):51-54.

[責(zé)任編輯：關(guān)金玉 英文編輯：劉彥哲]

Orientation of Sun Shadow Based on Global Search Algorithm

SUN Hui-yu1,HE Xiao-long2,WU Fei3

(1.School of Finance,Anhui University of Finance and Economics,Bengbu,Anhui 233030,China;2.School of Management Science and Engineering,Anhui University of Finance and Economics,Bengbu,Anhui 233030,China； 3.School of Statistics and Applied Mathematics,Anhui University of Finance and Economics,Bengbu,Anhui 233030,China)

Objective To work out where the measurement locates with the known length of the shadows and the date of survey or the date of survey with the known length of the shadows and where the measurement locates,thus ensuring the position with a known wild adventure and the measurement address according to a photo or video.Methods Aiming at the orientation of sun shadow,based on global search algorithm,the return model parameter estimation by ordinary least squares and maximum likelihood,and the bias elimination,the global search network model,and parameter equation model were establish.By using MATLAB,the algorithm writing,graphing,data fitting and predicting were conducted.Results Combing with the data of group 1,the result was Guangxi Province;as for the question 2,the answer Mongolia.Conclusion By using the parameter equation model,global search algorithm and MATLAB,the length of the shadow,the length of the pole,the measuring location,the data of survey and the time of survey can be worked out.

sun shadow positioning;global web search;measurement;MATLAB

教育部人文社會(huì)科學(xué)青年科研資助項(xiàng)目(10YJC630143)

孫慧宇(1994-)，女，江蘇南京人，安徽財(cái)經(jīng)大學(xué)金融學(xué)院在讀學(xué)生，研究方向?yàn)閼?yīng)用數(shù)學(xué)。

朱家明(1973-)，安徽泗縣人，副教授，碩士，研究方向：應(yīng)用數(shù)學(xué)與建模。

TP 242.6

A

10.3969/j.issn.1673-1492.2017.01.004

來稿日期：2016.04.05