孫慧宇，賀小龍，吳 飛

(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)

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

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

0 引 言

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

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

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

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

(1)太陽赤緯σ：表示太陽光線與地球赤道面的夾角，一年四季每天都在變動著，冬至日σ=-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)

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

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

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

圖1 太陽高度角示意圖 圖2 影長示意圖

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

(2)

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

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

(3)

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

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

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

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

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

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

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

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

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

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

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

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

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

圖4 直桿影長示意圖

根據(jù)文獻(xiàn)[4]中的數(shù)據(jù)利用Excel進(jìn)行二次擬合,可得到北京時(shí)間14∶42～15∶42時(shí)之間的影長與時(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ù)測可得到直桿影長示意圖(見圖4)。

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

2.2 求解方法

根據(jù)表1提供的21組不同時(shí)間對應(yīng)的影子坐標(biāo)，計(jì)算出相應(yī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í)間對應(yīng)求出，最后模型中只剩下測量地的緯度是未知的。

2.3 結(jié)果及分析

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

(6)

圖5 影長比例圖

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2 求解方法

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

3.3 結(jié)果分析

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

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

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

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

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

4 結(jié) 論

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

并以此模型為藍(lán)本，通過比例消元、參數(shù)估計(jì)以及基于窮舉法的全局搜索算法等方法，針對不同的觀測數(shù)據(jù)可以做到知四求一甚至知三求二。即測量地、測量時(shí)間、測量日期、物體影長和直桿長度中3個或3個以上為已知時(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]中國工業(yè)與應(yīng)用數(shù)學(xué)學(xué)會，2015全國大學(xué)生數(shù)學(xué)建模競賽A題[EB/OL].http://www.mcm.edu.cn/html.cn/node/ac8b96613522ef62c019dlcd45a125e3.html.2015-09-14.

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

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

[7]陳佳.周濤.曾祥平.城市異常路段交通擁擠——消散分析[J].中國水運(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].長春：吉林大學(xué)出版社，2005:4-59.

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

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

[12]肖智勇,劉宇翔.一種新的緯度測量方法[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

教育部人文社會科學(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