太陽影子定位模型

佘梓航，馬曉娜

（韓山師范學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，廣東潮州521041）

太陽影子定位模型

佘梓航，馬曉娜

（韓山師范學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，廣東潮州521041）

摘要：以2015年全國大學(xué)生數(shù)學(xué)建模競(jìng)賽A題的附件1和附件4為對(duì)象.通過建立太陽影子定位的模型，分析確定影子對(duì)應(yīng)物體所在的經(jīng)緯度以及拍攝時(shí)間.在問題1中，利用直角三角形的性質(zhì)，建立物體影子長(zhǎng)度關(guān)于太陽赤緯、觀測(cè)點(diǎn)緯度、真太陽時(shí)和物體高度變化模型，在MATLAB環(huán)境下繪制影子長(zhǎng)度與各參數(shù)的三維相關(guān)圖，得到其變化規(guī)律；在問題2中,基于問題1中建立的影子長(zhǎng)度變化模型，結(jié)合真太陽時(shí)與觀測(cè)點(diǎn)經(jīng)度的關(guān)系，建立影子長(zhǎng)度關(guān)于觀測(cè)點(diǎn)緯度和觀測(cè)點(diǎn)經(jīng)度關(guān)系模型，通過非線性擬合估計(jì)觀測(cè)點(diǎn)經(jīng)緯度；在問題3中，將視頻數(shù)據(jù)進(jìn)行轉(zhuǎn)化，獲取太陽影子頂點(diǎn)坐標(biāo)數(shù)據(jù),基于問題2中建立的影子長(zhǎng)度變化模型，結(jié)合日期與太陽赤緯的關(guān)系，通過非線性擬合估計(jì)視頻拍攝的可能地點(diǎn).

關(guān)鍵詞：控制變量法；數(shù)據(jù)可視化；非線性最小二乘；非線性擬合；真太陽時(shí)

太陽影子定位技術(shù)就是通過分析物體的太陽影子變化，確定視頻拍攝地點(diǎn)經(jīng)緯度和日期的一種方法.如何確定拍攝地點(diǎn)和拍攝日期是數(shù)據(jù)分析的重要方面，本文采用2015年全國大學(xué)生數(shù)學(xué)建模競(jìng)賽A題附件1和附件4的數(shù)據(jù)，嘗試解決以下問題.

問題1：建立影子長(zhǎng)度變化的數(shù)學(xué)模型，分析影子長(zhǎng)度關(guān)于各個(gè)參數(shù)的變化規(guī)律，并畫出2015年10 月22日北京時(shí)間9:00-15:00之間天安門廣場(chǎng)（北緯39°54'26"，東經(jīng)116°23'29"）3 m高的直桿的太陽影子長(zhǎng)度的變化曲線.

問題2：根據(jù)附件1中某固定直桿在水平地面上的太陽影子頂點(diǎn)坐標(biāo)數(shù)據(jù)，建立數(shù)學(xué)模型確定直桿所處的地點(diǎn).

問題3：根據(jù)附件4中一根2 m高直桿在太陽下的影子變化的視頻，建立數(shù)學(xué)模型確定若干個(gè)可能的視頻拍攝地點(diǎn).

注：題目及數(shù)據(jù)附件都可到全國大學(xué)生數(shù)學(xué)建模競(jìng)賽網(wǎng)站（www.mcm.edu.cn）下載.

1模型的假設(shè)

（1）假設(shè)物體較規(guī)則且物體直立；

（2）拍攝當(dāng)天天氣狀況良好，受太陽光持續(xù)照射；

（3）物體四周無其余物體遮擋，保證影子軌跡不受外界影響；（4）不考慮空氣折射和太陽視半徑.

2模型的建立與求解

2.1問題1分析與建模求解

問題1要求建立影子長(zhǎng)度變化的數(shù)學(xué)模型，分析影子長(zhǎng)度關(guān)于各個(gè)參數(shù)的變化規(guī)律，并繪制太陽影子長(zhǎng)度的變化曲線.本文利用控制變量法及數(shù)據(jù)可視化處理分析太陽影子長(zhǎng)度與各相關(guān)參數(shù)的變化關(guān)系建立模型1，繪制2015年10月22日北京時(shí)間9:00-15:00之間天安門廣場(chǎng)3 m高的直桿的太陽影子長(zhǎng)度的變化曲線.

2.1.1模型1的建立

太陽高度角的正弦值［1］計(jì)算如下：

根據(jù)物體影子長(zhǎng)度與太陽高度角的余切關(guān)系，影子長(zhǎng)度公式為：

根據(jù)（2）式，計(jì)算各個(gè)參數(shù)，最終得到影子長(zhǎng)度變化的余切模型：

其中L為物體影子長(zhǎng)度，H為物體的高度或長(zhǎng)度，φ表示當(dāng)?shù)鼐暥?，δ為太陽赤緯角，Ω為太陽時(shí)角.

2.1.2模型1的求解

保持太陽赤緯角δ不變，取3 m高的較規(guī)則物體，并假設(shè)其直立.結(jié)合（3）式，利用MATLAB編程，繪制該物體的影子長(zhǎng)度在太陽赤緯角δ為0度時(shí)，隨太陽時(shí)角和觀測(cè)點(diǎn)緯度變化的三維圖，如圖1所示.

由圖1觀測(cè)到：影子長(zhǎng)度在日出到日落的變化趨勢(shì)先是由長(zhǎng)變短再變長(zhǎng)，以太陽時(shí)角為0度即真太陽時(shí)為12點(diǎn)作為分界點(diǎn)，長(zhǎng)度最短.利用（3）式，將題目中所給數(shù)據(jù)代入對(duì)應(yīng)公式，得到天安門廣場(chǎng)2015年10月22日北京時(shí)間9:00-15:00之間天安門廣場(chǎng)（北緯39°54'26"，東經(jīng)116°23'29"）3 m高的直桿的太陽影子長(zhǎng)度，繪制太陽影子長(zhǎng)度與真太陽時(shí)的曲線圖.

圖1　影子長(zhǎng)度隨太陽時(shí)角和觀測(cè)點(diǎn)緯度變化的三維圖

在（3）式中，求影子長(zhǎng)度需要得到太陽赤緯角δ和太陽時(shí)角Ω，其計(jì)算步驟如下.

（1）赤緯角的計(jì)算

根據(jù)有關(guān)文獻(xiàn)，太陽赤緯角的計(jì)算公式［2］為：

其中，θ表示日角，計(jì)算公式為：

此處t由兩部分組成，即：

式中N為積日，積日表示日期在年內(nèi)的順序號(hào).且：

INT（X）為取不大于X的最大整數(shù)，Y為年份.

（2）太陽時(shí)角的計(jì)算

規(guī)定正午太陽時(shí)角Ω為0，上午Ω為負(fù)值，下午Ω為正值.由于地球自轉(zhuǎn)一周360°需要24小時(shí)，每小時(shí)對(duì)應(yīng)時(shí)角為15°.太陽時(shí)角的計(jì)算公式為：

其中tO為真太陽時(shí)，以24小時(shí)計(jì).在某地要取得正確的當(dāng)?shù)卣嫣枙r(shí)tO，由北京時(shí)間來換算需要修正兩種差.

由于各地經(jīng)度不同而存在時(shí)間差.用Bt表示北京時(shí)間，用Et表示時(shí)差，那么平太陽時(shí)tm為：

其中時(shí)差與當(dāng)?shù)亟?jīng)度有關(guān)，其計(jì)算公式為：

其中λ表示當(dāng)?shù)亟?jīng)度，時(shí)差Et的單位為小時(shí).

平太陽時(shí)tm假設(shè)地球繞太陽是標(biāo)準(zhǔn)的圓形，即一年中每天都是24小時(shí)平均的.但地球的公轉(zhuǎn)軌道是一個(gè)橢圓，也就是說地球的公轉(zhuǎn)并不均勻，并非每天都為24小時(shí).因此要將平太陽時(shí)調(diào)整為真太陽時(shí)，需要修正真太陽時(shí)與平太陽時(shí)之差η.時(shí)差只與日期相關(guān)，與地區(qū)所處位置無關(guān).即：

其中：

式中，t表示每年自1月1日后的天數(shù)，其中不包括1月1日.η單位為分鐘.

根據(jù)上述公式，且真太陽時(shí)為平太陽時(shí)與真平太陽時(shí)差相加之和［3］，得到真太陽時(shí)為：

結(jié)合（3）式，應(yīng)用以上的求解步驟，利用MATLAB編程實(shí)現(xiàn)，繪制2015年10月22日北京時(shí)間9:00-15:00之間天安門廣場(chǎng)3 m高的直桿的太陽影子長(zhǎng)度的變化曲線.以3 min為時(shí)間間隔取點(diǎn)，計(jì)算120個(gè)時(shí)間點(diǎn)各個(gè)參數(shù)的值，最終得到直桿影子長(zhǎng)度隨時(shí)間變化的圖像如圖2所示.

由圖2曲線可知，天安門廣場(chǎng)3 m高的直桿的太陽影子在指定的時(shí)間段北京時(shí)間9:00-15:00內(nèi)，隨著時(shí)間變化，影子長(zhǎng)度總體呈現(xiàn)出先變短后變長(zhǎng)的趨勢(shì)，在9:00時(shí)，影子長(zhǎng)度值達(dá)到最高點(diǎn)為6.438 6 m；接著長(zhǎng)度值逐漸減少，約在12:00處到達(dá)最低點(diǎn)；12:00過后長(zhǎng)度值曲線回升，最后在15:00處.

2.2問題2分析與建模求解

問題2要求根據(jù)某固定直桿在水平地面上的太陽影子頂點(diǎn)坐標(biāo)數(shù)據(jù)（附件1），建立數(shù)學(xué)模型確定直桿所處地的緯經(jīng)度.由（3）式可知影子長(zhǎng)度與參數(shù)桿長(zhǎng)H，緯度φ，太陽赤緯角δ，真太陽時(shí)tO存在確定的非線性關(guān)系且真太陽時(shí)與觀測(cè)點(diǎn)經(jīng)度存在確定關(guān)系，此問題可以轉(zhuǎn)化為非線性最小二乘估計(jì)問題.已知數(shù)據(jù)的測(cè)量時(shí)間為2015年4月18日，可以算出太陽赤緯.利用非線性最小二乘法［4］擬合估計(jì)出觀測(cè)點(diǎn)的經(jīng)度和緯度，再根據(jù)經(jīng)緯度在地圖上查詢直桿的具體位置.

2.2.1模型2的建立

真太陽時(shí)與經(jīng)度關(guān)系為：

圖2　直桿影子長(zhǎng)度曲線變化圖

其中Bt表示北京時(shí)間，η為真平太陽時(shí)時(shí)差，且：

太陽時(shí)角與真太陽時(shí)關(guān)系：

將（14）和（15）式代入（16）式可得太陽時(shí)角與經(jīng)度的關(guān)系式：

要使擬合殘差平方之和最小，即目標(biāo)函數(shù)為：

其中n=2，Li為直桿真實(shí)影長(zhǎng)，為每一次估計(jì)的直桿擬合影長(zhǎng)，E表示殘差平方之和計(jì)算如下：

綜合以上結(jié)果，可以得到使得殘差平方之和最小的最優(yōu)模型：

2.2.2模型2的求解

Step2：利用MATLAB編程調(diào)用非線性擬合函數(shù)lsqnonlin擬合參數(shù)桿長(zhǎng)L，緯度φ，觀測(cè)點(diǎn)經(jīng)度λ的估計(jì)值，且設(shè)置的參數(shù)初始值，可得到多個(gè)擬合結(jié)果.

Step3：繪制擬合圖，利用殘差檢驗(yàn)擬合優(yōu)度.若擬合效果差，應(yīng)返回Step2，重新擬合；擬合效果好，通過檢驗(yàn).仍返回Step2重新設(shè)置區(qū)間，得到多組可能的觀測(cè)點(diǎn)經(jīng)緯度，如表1所示.

表1　直桿可能的所在地

根據(jù)本文計(jì)算的數(shù)據(jù)所得的擬合圖如圖3所示.由圖3可知，影長(zhǎng)真實(shí)值和估計(jì)值基本吻合，擬合優(yōu)度較好.

2.3問題3分析與建模求解

本題要求根據(jù)一根2 m高的直桿在水平地面上的太陽影子變化的視頻（附件4），建立數(shù)學(xué)模型確定拍攝地點(diǎn)經(jīng)緯度.本題給的是視頻數(shù)據(jù)，則需要將視頻的數(shù)據(jù)進(jìn)行轉(zhuǎn)化，獲取太陽影子頂點(diǎn)坐標(biāo)數(shù)據(jù)，建立太陽高難度角相對(duì)于觀測(cè)點(diǎn)緯度和經(jīng)度的變化模型，其中3個(gè)參數(shù)分別為：太陽高度角h，觀測(cè)點(diǎn)緯度φ，觀測(cè)點(diǎn)經(jīng)度λ，可以利用非線性最小二乘法擬合估計(jì)，得到經(jīng)度和緯度數(shù)據(jù)，從而直桿所在地得以確定.

圖3　影子長(zhǎng)度與真實(shí)值的擬合圖

2.3.1模型3的建立

由于桿長(zhǎng)和時(shí)間是已知的，太陽赤緯可以根據(jù)時(shí)間代入公式（4）～（7）求出；通過讀取視頻處理得到的坐標(biāo)數(shù)據(jù)和已知的桿長(zhǎng)，可以得到太陽高度角，結(jié)合（1）式和（17）式得到太陽高度關(guān)于觀測(cè)點(diǎn)緯度與經(jīng)度的關(guān)系式為：

此問中利用太陽高度角和觀測(cè)點(diǎn)緯度的非線性最小二乘估計(jì)，要使擬合殘差之和最小，則有基本模型：

2.3.2模型3的求解

Step1：利用MATLAB軟件中的VideoReader函數(shù)，將視頻轉(zhuǎn)化為圖片，由于視頻的時(shí)長(zhǎng)約為40 min，且附件1～3中均以21個(gè)數(shù)據(jù)點(diǎn)作為分析數(shù)據(jù).不失一般性，每間隔2 min取圖片，從而獲取21張圖片.

Step2：在MATLAB環(huán)境下，調(diào)用rgb2gray編程實(shí)現(xiàn)將Step1獲取的21張真彩色圖像轉(zhuǎn)換為灰度圖像，記錄下每一時(shí)間間隔桿頂影子的像素坐標(biāo)a（x1，y1），桿頂?shù)南袼刈鴺?biāo)b（x2，y2），桿底的像素坐標(biāo)c（x3，y3）.

Step3：將Step2得到的桿頂影子像素坐標(biāo)a（x1，y1）轉(zhuǎn)化為以桿底像素坐標(biāo)為原點(diǎn)的像素坐標(biāo)（x1-y3，y1-y3），對(duì)附件4中的數(shù)據(jù)進(jìn)行處理，結(jié)果如表2所示.

表2　桿頂點(diǎn)陰影的像素坐標(biāo)數(shù)據(jù)　（像素）

Step4：根據(jù)桿頂點(diǎn)與底端的像素坐標(biāo)，其縱坐標(biāo)之差為d=y3-y2，結(jié)合題目中給出桿的高度為2 m，則可將桿頂影子像素坐標(biāo)轉(zhuǎn)化為以米為單位的坐標(biāo)，即可得桿影的頂點(diǎn)坐標(biāo)數(shù)據(jù)，將附件4中的數(shù)據(jù)代入，以桿底為原點(diǎn)坐標(biāo)，結(jié)果如表3所示.

表3　桿影頂點(diǎn)坐標(biāo)數(shù)據(jù) （m）

類比模型2，對(duì)附件4中視頻數(shù)據(jù)轉(zhuǎn)化后，結(jié)合 （21）式，通過MATLAB編程實(shí)現(xiàn)，給出直桿可能的地點(diǎn)，見表4.

表4　拍攝地點(diǎn)經(jīng)緯度

2.3.3誤差分析

通過分析得到的影子坐標(biāo)計(jì)算影子長(zhǎng)度L，則在其基礎(chǔ)上添加正態(tài)隨機(jī)數(shù)Z（0，0.005，1，21），即影子長(zhǎng)度"，再進(jìn)行非線性擬合求解，對(duì)擬合的結(jié)果進(jìn)行白噪聲檢驗(yàn)［5］，影子長(zhǎng)度添加正態(tài)隨機(jī)數(shù)后，R值（0.999）較為接近，且地點(diǎn)多為內(nèi)蒙古，則說明通過視頻讀取的數(shù)據(jù)具有合理性.

3結(jié)語

對(duì)于問題1，本文給出了影子長(zhǎng)度同各個(gè)參數(shù)之間的變化規(guī)律關(guān)系圖，為解答接下去的問題提供依據(jù).在畫出影子長(zhǎng)度變化圖時(shí)，能夠考慮到晚上沒有影子這個(gè)事實(shí)，畫出來的圖更具有代表性.同時(shí)，利用非線性最小二乘法給出的結(jié)果，接近數(shù)模組委會(huì)給出來的結(jié)果，正確率高.對(duì)于最后一個(gè)問題，由于視頻中給出的影子長(zhǎng)度受到很多不確定性因素的影響（如：拍攝角度，光的折射等），因此，最后一個(gè)問題不可能給出準(zhǔn)確的結(jié)果.即便如此，為了驗(yàn)證本文所給方法的正確性，本文通過對(duì)影子坐標(biāo)數(shù)據(jù)加入白噪聲，再進(jìn)行非線性擬合求解，得到同沒加白噪聲相近的擬合精度，說明了從視頻中讀取數(shù)據(jù)具有合理性.

參考文獻(xiàn)：

［1］鄭鵬飛，林大鈞，劉小羊，等.基于影子軌跡線反求采光效果的技術(shù)研究［J］.華東理工大學(xué)學(xué)報(bào)，2010，36（3）：459-460.

［2］陳曉勇，鄭科科.對(duì)建筑日照計(jì)算中太陽赤緯角公式的探討［J］.浙江建筑，2011，28（9）：7-8.

［3］孫吉山，洪薇.北京標(biāo)準(zhǔn)時(shí)間與各地真太陽時(shí)的換算法［J］.上海針灸雜志，1994，13（1）：25-26.

［4］薛毅.最優(yōu)化理論與算法［M］.北京：北京工業(yè)大學(xué)出版社，2004：201-202.

［5］王燕.應(yīng)用時(shí)間序列分析［M］.3版.北京：中國人民大學(xué)出版社，2012：27-31.

（責(zé)任編輯：邵曉軍）

中圖分類號(hào)：O29

文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1007-5348（2016）04-0006-07

［收稿日期］2016-03-12

［作者簡(jiǎn)介］佘梓航（1989-），男，廣東潮州人，韓山師范學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院助教，碩士；研究方向：高頻數(shù)據(jù)估計(jì)量構(gòu)建與數(shù)學(xué)建模.

Position Models of the Sun Shadow

SHE Zi-hang,MA Xiao-na
（DepartmentofMathematicsandStatistics,HanshanNormal University,Chaozhou521041,Guangdong,China）

Abstract:In this paper,a solution is presented for the ProblemA of Contemporary Undergraduate Mathematical Contest in Modeling in 2015.In question 1,by using the property of right triangle,it built the sun shadow model about the solar declination,latitude,apparent solar time and the length of object.It used MATLAB todrawthe 3D correlation diagramof the parameters and shadowlength.In question 2,based on the model of question 1 and the relationshipbetweentheapparent solar timeandthelongitudeof theobservationpoint,it built lengthof sunshadow model about the latitude and longitude by using the observation point.Therefore,it estimated the latitude and longitudeby NonlinearFittingMethod（NFM）.Inquestion3,itusedthevideodatatoobtainthevertex coordinatesof the sun shadow.Based on and the relationship between the date and the solar declination,it estimated the possible locationsofthevideobyNFM.

Key words:control variable method;data visualization;nonlinear least square method;nonlinear fitting method （NFM）;solardeclination