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

孫慧宇，賀小龍，吳 飛

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

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

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

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

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

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

0 引 言

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

1 太陽影子的影響機制

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

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

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

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

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

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

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

(1)

其中α為太陽高度角、φ為測點緯度、σ為太陽赤緯，t為太陽時角。

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

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

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

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

(2)

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

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

(3)

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

桿長、測地緯度、太陽時角保持不變時，則太陽直射點離測量地越遠，影長越短；

桿長、太陽赤緯、太陽時角均保持不變時，則測量地離太陽直射點越遠，影長越短；

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

圖3 直桿的影長隨時間變化圖

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

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

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

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

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

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

表1 直桿影子定點坐標表

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

圖4 直桿影長示意圖

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

y=0.0004x2+0.0305x+1.12083

(5)

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

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

2.2 求解方法

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

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

2.3 結(jié)果及分析

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

(6)

圖5 影長比例圖

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

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

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

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

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

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

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

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

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

表3 觀測地坐標及對應影子長度匯總表

注:定義x=1對應文獻[4]中的起始觀測時刻，時間每增加3min，x的值增加1。

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

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

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

推算出相應的赤緯角度y。

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

3.2 求解方法

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

3.3 結(jié)果分析

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

圖6 影子長度的散點圖 圖7 影子長度的擬合曲線圖

圖8 影子長度的散點圖 圖9 影子長度的擬合曲線圖

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

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

4 結(jié) 論

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

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

[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è)與應用數(shù)學學會，2015全國大學生數(shù)學建模競賽A題[EB/OL].http://www.mcm.edu.cn/html.cn/node/ac8b96613522ef62c019dlcd45a125e3.html.2015-09-14.

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

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

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

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

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

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

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

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

[責任編輯：關(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

教育部人文社會科學青年科研資助項目(10YJC630143)

孫慧宇(1994-)，女，江蘇南京人，安徽財經(jīng)大學金融學院在讀學生，研究方向為應用數(shù)學。

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

TP 242.6

A

10.3969/j.issn.1673-1492.2017.01.004

來稿日期：2016.04.05