太陽影子定位,2015數(shù)學建模國賽A題

對太陽影子定位算法探究摘要本文是對2013年全國大學生數(shù)學建模競賽A題的解答.隨著人們對數(shù)據(jù)挖掘的深入，如何確定視頻的拍攝地點和拍攝日期已經(jīng)成為視頻數(shù)據(jù)分析的重要方面，太陽影子定位技術就是通過分析視頻中物體的太陽影子變化，確定視頻拍攝的地點和日期的一種方法，進而可以促進視頻分析定位技術發(fā)展。對于問題一，我們根據(jù)地球自轉公轉的自然規(guī)律，建立影子長度變化的數(shù)學模型，并且分析影子長度關于各個參數(shù)的變化規(guī)律?；趯栴}的分析以及理論的學習研究，畫出模擬概念圖，然后計算相關量（如太陽高度角、赤緯角等）的表達式，并按照相關地理知識建立起模型，得到桿子影長與時間函數(shù)表達式，再將題目所給的數(shù)據(jù)代入求解方程，并用MATLAB作出曲線圖，最后檢驗模型的準確性。對于問題二，我們以問題一所求出的表達式和資料作為基礎，繼而利用球面天文學求算太陽視坐標的簡化算法建立一模型直接求解出經(jīng)度，緯度的估算值。再代進數(shù)據(jù)并用利用多項式擬合出更長的時間序列曲線，用函數(shù)的特征值（最低點）加上時角，時區(qū)計算相關知識，再推算出經(jīng)度值。最后利用第一問模型，經(jīng)度，加上曲線獲得的幾組影長數(shù)據(jù)聯(lián)立求解出大致緯度，最后估算桿子所在的地區(qū)。對于問題三，結合問題一問題二所建立的模型，將附件2，附件3的數(shù)據(jù)先畫出散點圖并以多項式擬合出兩條相對完整的曲線，通過其曲線函數(shù)求得影長的最小值以及最小值所對應的時間求得經(jīng)度，緯度，將經(jīng)度和緯度代入赤緯角公式以及影長公式可求得相應的具體日期。對于問題四，首先將視頻數(shù)據(jù)利用MATLAB，并且編程處理視頻得到每分鐘一幀的圖片，再把相關圖片轉化為灰度圖矩陣，最后用語句轉化為二值圖（0為黑，1為白）。下一步把二值圖集分析并且分析出桿子影長的變化規(guī)律，求出視頻拍攝點經(jīng)度，利用模型一求出緯度，即是位置。關鍵字：影長位置MATLAB編程多項式擬合最小二乘法二值圖一．問題重述如何確定視頻的拍攝地點和拍攝日期是視頻數(shù)據(jù)分析的重要方面，太陽影子定位技術就是通過分析視頻中物體的太陽影子變化，確定視頻拍攝的地點和日期的一種方法。1.建立影子長度變化的數(shù)學模型，分析影子長度關于各個參數(shù)的變化規(guī)律，并應用你們建立的模型畫出2015年10月22日,北京時間9:00-15:00之間天安門廣場,3米高的直桿的太陽影子長度的變化曲線。2.根據(jù)某固定直桿在水平地面上的太陽影子頂點坐標數(shù)據(jù)，建立數(shù)學模型確定直桿所處的地點。將你們的模型應用于附件1的影子頂點坐標數(shù)據(jù)，給出若干個可能的地點。3.根據(jù)某固定直桿在水平地面上的太陽影子頂點坐標數(shù)據(jù)，建立數(shù)學模型確定直桿所處的地點和日期。將你們的模型分別應用于附件2和附件3的影子頂點坐標數(shù)據(jù)，給出若干個可能的地點與日期。4．附件4為一根直桿在太陽下的影子變化的視頻，并且已通過某種方式估計出直桿的高度為2米。請建立確定視頻拍攝地點的數(shù)學模型，并應用你們的模型給出若干個可能的拍攝地點。如果拍攝日期未知，你能否根據(jù)視頻確定出拍攝地點與日期？注：題目及數(shù)據(jù)附件都可到全國大學生數(shù)學建模競賽官方網(wǎng)站（）下載.二．模型的假設假設地球在自轉時角速度均勻；假設地球在公轉時的橢圓為規(guī)則橢圓；假設題目涉及的太陽時計算為真太陽時；假設在一天內(nèi)的太陽高度角不變；忽略大氣折射對陽光傳播的影響。注：1真太陽日分為24真太陽時。這個時間系統(tǒng)稱為真太陽時。真太陽時也稱為視太陽時，簡稱視時。真太陽時即真太陽視圓面中心的時角加12小時。即：真太陽時=平太陽時+真平太陽時差符號說明表1符號符號說明赤緯角方位角太陽高度角太陽時角(地球每個小時自轉15°稱之為時角)W地理緯度J經(jīng)度L影長D日期H桿子高度n某點的地方時間L桿影長度四．模型的建立與求解4.1問題一4.1.1問題1分析題目要求描建立影子長度變化的數(shù)學模型，分析影子長度關于各個參數(shù)的變化規(guī)律。首先通過查閱相關地理、物理文獻，找尋建模中可能遇到的概念、方法，再挖掘建立模型相關的方法，為之后的建模奠定堅實的基礎。模型建立第一步,理清自變量（四個），因變量，選取參考系（參考面還有參考點），建立空間直角坐標系。再做出太陽隨時間移動桿影的直觀視圖（日照桿影圖），最后，利用太陽高度角、桿高等相關量與桿影的幾何關系，構造對應函數(shù)，使之符合題意，同時根據(jù)相應的標準進行統(tǒng)計、分析和構建相關的數(shù)學模型。4.1.2模型1的建立通過對題目的探討，一根垂直地面的桿子，研究它的影子長度隨著日照光線的改變而改變，首先我們確定相應的參考系，以桿子所在地地平面α為參考面，以桿子的最低端為坐標原點，再得到以下的自變量：表2名稱符號時間n竿子高度（對地）H位置經(jīng)度J緯度W因變量：桿子在參考面上投影長度（可以從桿子頂點在參考面的投影求得）4.1.3步驟一：求出建立模型必須的物理，地理相關量值量在此基礎上我們分析地球運動的特性規(guī)律與研究桿影與時間、日期、經(jīng)度、緯度、桿長的相關關系，首先做出日照桿影分析圖1-1，如下:圖1-1設桿子頂點在地平面上投影點相對于原點坐標為Q(x0,y0),另外桿子高度為H，影子長度為L。太陽高度角可求：;從圖1-1中分析，太陽與桿子相夾的角近似為太陽高度角，太陽赤緯角：定義是又稱太陽赤緯，是地球赤道平面與太陽和地球中心的連線之間的夾角。由圖1-1幾何關系，查閱資料以及推理出計算公式如下：因為題中涉及到的是真太陽日所以太陽赤緯角的計算公式為①k為日期序號，如若3月22號,k就為81。時角計算公式：②太陽高度角θ計算公式：③根據(jù)三角函數(shù)知識以及圖中幾何關系得，太陽方位角計算公式：④4.1.4步驟二：利用上面所得的數(shù)據(jù)求出Q點的軌跡方程（即為影子軌跡）圖1-2對圖1-2進行分析，發(fā)現(xiàn)在某個地點某個時刻太陽照射桿子桿影的變化曲線圖的大致規(guī)律：一條類似于拋物線的曲線對桿影的因變量和自變量進行量化分析得影子長度,推并推算相關公式得：⑤⑥然后將上述式子聯(lián)立求解得Q點的x坐標的表達式：⑦代入數(shù)據(jù)得由再以時間t為自變量，可求得影子長度的變化曲線：如果數(shù)據(jù)充足，可以利用MATLAB做出Q點的軌跡，以及影子長度隨著自變量的變化而變化曲線圖，其中緯度被包含在公式中，經(jīng)度被時角所替代，地方時間n，日期也嵌套在赤緯角里面，這樣四個自變量就可以表示桿影的長度變化曲線。4.1.5步驟三：計算取值，作曲線計算畫出2015年10月22日北京時間9:00-15:00之間天安門廣場（北緯39度54分26秒,東經(jīng)116度23分29秒）3米高的直桿的太陽影子長度的變化曲線。由步驟二可知道,,.是主要公式，聯(lián)合步驟二的有關公式可以得出直桿的太陽影子長度隨時間的變化曲線為.⑩其中cosθ是時間的函數(shù)，因為H、W、J、n、日期已知，代入數(shù)據(jù)并用MATLAB求得桿子頂點的影子長度的變化曲線如曲線圖所示：圖1圖1-3檢驗準確性分析圖1-3，可見曲線直觀上是一條類似于拋物線的曲線，表明3米高的桿子在北京天安門廣場9:00到15:00的桿影長度曲線如圖1-3。圖1-4來自于·影端軌跡周年變化的實踐與分析_以北溫帶地區(qū)為例_吳濟廉再根據(jù)相關文獻中的依據(jù)，其中包含實驗驗證此現(xiàn)象的記載，表明結果的準確性程度較高。4.2問題二4.2.1問題二的分析題目要求根據(jù)某固定直桿在水平地面上的太陽影子頂點坐標數(shù)據(jù)，建立模型確定該桿所在的地理方位。附件一中有的條件是影長，時間，日期，不明確竿子高度H，求所測地點的經(jīng)緯度。4.2.2模型的建立與檢驗：設竿子長度為H（未知），以問題一所建立模型為基礎，利用球面天文學求算太陽視坐標的簡化算法對問題進行探究，建立此問題的數(shù)學模型，并且?guī)нM附件一中的數(shù)據(jù)，再利用數(shù)據(jù)擬合出一天之中的竿影長度隨時間的變化圖，得出最短陰影長度值，進而得到該點經(jīng)度值。估算出測驗地的地理未知，并檢驗。首先我們要計算某一時間的天頂距S（是指在天體方位圈上，天體與天頂之間的交角稱為天頂距來自百度百科），可以從天頂?shù)亩ㄎ蝗切?利用三角幾何關系求圖2-1百度百科上的天頂角示意圖天頂距S的代數(shù)關系式子：?時角的精確計算方式：?a0為視赤經(jīng)，f為真太陽時（忽略），r為地方時，為視赤緯；,太陽的視天頂距為利用最小二乘法進行擬合求觀測地點的經(jīng)緯度，桿子的高度與影長之間的關系:?當影長一定時，假設桿子高度H與經(jīng)度緯度的函數(shù)，即是當影長一定時,可把由影長而推算的物高H看作為緯度和經(jīng)度的函數(shù)?結合第一問所得的L-t圖，解微分方程求解經(jīng)度大約106°，再用多項式擬合求經(jīng)度與之比較。代入已知數(shù)據(jù)之后，用多項式擬合法擬合出出桿影隨時間的變化而變化圖像，如下:圖4-1圖4-1陰影部分為實際區(qū)，藍色部分為預測區(qū)最低點是：擬合度的檢驗，利用MATLAB進行regress擬合度命令的多次分析，擬合度b的數(shù)值的平均值為0.9998723，證明在區(qū)域內(nèi)，曲線的擬合效果優(yōu)良。接下來就是分析計算經(jīng)緯度值，由擬合出來的曲線可以計算桿子影子最短時的時間，擬合方程所代表式子為基礎條件，求經(jīng)度：測量最短桿影，地球自轉一周需要24小時,也就是1440分鐘,地球自轉360度,因此,可以算出每4分鐘,地球自轉一度。用上面的時間和正午12點之差,就可以求出你的經(jīng)度與當?shù)貢r區(qū)經(jīng)度的差。最短竿影下的時間約為12.60時即是12:36左右，查閱資料獲得東八區(qū)的經(jīng)度為120°,可知桿子所在地的經(jīng)度為.解得經(jīng)度為：120°±9°E即是111°E或者129°E.取最短影長為參考值，此時時間為12：36，日期為2015年4月8日，由第一問提取到的公式聯(lián)立有：.??綜合解得：緯度的可能取值是40.55°N和18.25°N。在問題一建立的模型的基礎上，通過對附件1所給的數(shù)據(jù)進行多項式擬合出時間序列的二次函數(shù)曲線，并通過驗證擬合度的方式保證該函數(shù)曲線的可靠性。接著通過分析二次函數(shù)最小值點以求得最小影長，以及最小影長所對應的時間。將對應的時間代入求得經(jīng)度。結合最小影長和經(jīng)度可求得緯度。4.2.3問題二的求解由所求的函數(shù)曲線可得出其最小影長所對應的時間為12.598，結合地球自轉360度，即每4分鐘地球自轉一度，用所求時間與正午12點之差即可求得該地的經(jīng)度與東八區(qū)的經(jīng)度差，再由東八區(qū)120度可求得兩個位置的經(jīng)度分別為111.03°E和128.97°E，再通過模型一所建立的模型以及公式可求得緯度為40.55°N和18.25°N，由GPS定位可得該地理位置在呼和浩特和海南三亞附近地區(qū)。4.3問題三4.3.1問題3的分析題目3與題目2較為相似，盡管題目3多出了另一個附件的數(shù)據(jù)，但是這兩個附件的數(shù)據(jù)都沒有給出確切測定的日期。然而緯度，經(jīng)度的值與最小影長有關，即因通過求得最小影長去求得相應的經(jīng)度和緯度。又因為影響赤緯角公式以及影長公式中有日期，即將求得的經(jīng)度，緯度代入公式中進而求得日期。4.3.2模型3建立在前兩問的基礎上，我們?nèi)匀徊捎枚囗検綌M合的方式對附件2，附件3的函數(shù)進行處理，同時通過求出擬合度的方式驗證其求得函數(shù)的可靠性。為求得該經(jīng)度，緯度，仍然是在第二問的基礎上求得影長的最小值以及所對應的時間。由赤緯角以及影長公式可知其與具體日期有關，故將經(jīng)度，緯度代入即可求得具體日期。4.3.3模型3的求解首先對附件2和附件3的數(shù)據(jù)利用MATLAB進行多項式擬合，Y=0.0981x^2-2.9883x+23.3691Y=0.0981x^2-2.9883x+23.3691圖圖4-1陰影部分為實際區(qū)，藍色部分為預測區(qū)YY=0.2964x^2-7.2839+48.2661圖圖4-2陰影部分為實際區(qū)，藍色部分為預測區(qū)然后對兩個曲線擬合模型進行regress擬合度分析，都得到的擬合度匹配程度較高。由附件2，附件3的函數(shù)曲線可求得最小的影長為15.23和12.29，并代入求得緯度為40.40°N，經(jīng)度為80.24°E，地點在新疆阿瓦提縣附近，日期為6月18日附近。4.4問題四4.4.1問題4的分析對問題四的分析中，題目中已知桿高為2，另外題目提供一個40分鐘視頻附件，我們能從視頻得到影子長度隨時間的變化規(guī)律，這樣已知數(shù)據(jù)包含時間，日期，桿子長度，以及桿影隨時間的變化曲線，根據(jù)這些已知，利用模型一公式，可以求出視頻拍攝地的位置。重要的是如何分析提取視頻里面的數(shù)據(jù)內(nèi)容，因為視頻中的是立體圖，所得到的矩陣數(shù)列要經(jīng)過修正才能代表桿子的長度，才能進行位置的估計。4.4.2模型4建立首先我們把視頻導入MATLAB中,用附件四的視頻導入MATLAB中，然后取每N幀大概每隔一分取一張圖片，存入數(shù)據(jù)庫矩陣A，再利用MATLAB編程將所得的所有彩色圖片簡化成二維矩陣的灰度圖片以png格式保存（如圖4-1左），下一步將灰度圖片化成二值圖片（如4-2右），得到簡化的矩陣組合，最后分析矩陣得到桿影與時間的大致變化曲線。.分析過程如圖4-1圖4-1圖4-2圖4-2二值矩陣的建立為防止灰度值流失，Matlab計算的灰度值會被限制在0-255之間。在模型計算中，為方便確定模型中桿長及桿影灰度矩陣，需將灰度值矩陣中桿長及桿影部分轉化成二值矩陣，具體操作：對于桿長的灰度矩陣：對于隨時間變化的桿影灰度矩陣：灰度矩陣的處理所以桿影長度與時間的函數(shù)曲線如圖4-3所示 利用MATLAB對灰度矩陣進行運算，求解出桿長、桿影二值化矩陣邊緣，從而確定合適的。邊界定位模型為了得到附件4中隨時間變化的桿影長度，采用邊界定位法。它是基于二值化矩陣的兩極化特性：二值化后的桿影灰度圖像的右上邊界的255灰度值與左下邊界的255灰度值即為桿影的邊界，即桿影的長度為兩灰度值之間的距離。 邊界定位模型的具體步驟如下： 步驟一：右上邊界模型的構建 將第張圖片的右上邊界進行右上邊界定位，即求第張圖片矩陣從第行由下向上，第列由左向右開始，檢測灰度值首先為255的元素，并提取其位置數(shù)據(jù)——。 步驟二：左下邊界模型的構建 將第張圖片的左下邊界進行左下邊界定位，即求第張圖片矩陣從第行由上向下，第列由右向左開始，檢測灰度值首先為255的元素，并提取其位置數(shù)據(jù)——。 步驟三：桿影長度的確定將提取出的第張圖片的右上邊界位置數(shù)據(jù)與左下邊界位置數(shù)據(jù)進行計算：(1)4.4.3模型4的求解以方向北為基準方向，結合模型一，已知條件為時間，日期，；可求得最短桿影對應時間為，進而求出視頻拍攝經(jīng)度設L（t）=再利用相應時間與表7.1中數(shù)據(jù)按二次曲線進行擬合，使用最小二乘法確定擬合方程為其中，L為視頻中桿影長度，t為當?shù)貢r間。因為桿影最短的時刻即為當?shù)卣鐣r刻，故t=14.939，所以初始經(jīng)度為其中，J為視頻所在地初始經(jīng)度，tmin為視頻所在地桿影最短時刻。之后代進模型一求得緯度為北緯38.8640°N或者41.325°N。從GPS搜索得拍攝地大約在寧夏回族自治區(qū)銀川市。4.5額外問題如果拍攝日期未知，你能否根據(jù)視頻確定出拍攝地點與日期？如果缺少日期條件，采取建立基于遺傳算法的太陽影子定位模型，可以對經(jīng)緯度矩陣增加日期緯度，引入3個0-1變量，對應著該矩陣的經(jīng)度、緯度列向量、日期列向量，考慮到可能有若干個地點和日期，以各0-1變量求和不超過3作為約束條件，以竿影長和實際竿影長的的相似度作為目標函數(shù)，擬建立基于太陽影子定位的目標優(yōu)化模型。五、結論通過以上四個問題的分析解答探討，本文研究了桿子的影子在太陽照射下的影長與時間，日期、桿長以及桿子所在地的位置五者之間的相互關系，進而研究分析視頻中物體的太陽影子變化，確定視頻拍攝的地點和日期或者反推，這幾個的關系如圖。問題一中，我們根據(jù)所知道的桿長，日期，時間，以及位置，求得影長。再問題二、三則分別在未知桿長，未知桿長和日期的情況下求解出桿子的地理方位，而問題四則利用視頻中的影長變化來推斷拍攝點位置。將這個模型推廣，就有利于視頻定位系統(tǒng)的建立。只要知道其他的三個或更多量就可以計算另外的量，這對本文得到的模型方法應用型較強。六、模型的改進與推廣分析幾個模型后可以發(fā)現(xiàn)，模型有一定缺陷性，下面逐一分析：其一，我們沒有考慮大氣對太陽光的折光率（大氣層從下到上可分為對流層、平流層、電離層和磁層四層。無線電波在大氣層中傳播時，由于在各層中的傳播速度變化產(chǎn)生影響）問題，這對相關量，所測得的目標角度、距離、高度都存在大氣折射誤差，如太陽高度角的計算是有相當大的影響的。改善的方法是利用折射系數(shù)對太陽光進行修正，然后再求相關角度，最后代進模型。第二是在問題二、三中，我們利用多項式擬合的方法有相當大的不準確性，在數(shù)據(jù)提供范圍外，預測曲線與實際曲線的擬合度不是太高，導致一定的誤差性，對此我們利用多次擬合度分析，并用實際曲線貼近擬合實際曲線，盡量減少誤差量。關于模型的推廣，此模型適用于太陽影子定位的相關計算中，將這個模型推廣，就有利于視頻定位系統(tǒng)的建立。只要知道其他的三個或更多量就可以計算另外的量，因此，此模型還可以應用于對已知視頻定位，求物體的高度，求相關日期時間等等問題。參考文獻：[1]司守奎，孫璽菁:數(shù)學建模算法與應用[S]:國防工業(yè)出版社,2013.[2]何銀濤，張梅黃華.垂直單軸跟蹤光伏支架方陣間距計算界面設計.太陽能報刊.2015.[3]鄭鵬飛,林大鈞,劉小羊,吳志庭.基于影子軌跡線反求采光效果的技術研究[].廣東:華南理工大學200237.[4]林根石,利用太陽視坐標的計算進行物高測量與定位[A].南京大學學報,1991.[5]吳濟廉,影端軌跡周年變化的實踐與分析.地理教學報，2013,50(5):9-20.[6]張麗梅,趙建立,喬立山.用二值矩陣表示研究格矩陣的{1}-廣義逆和{1,2}-廣義逆.中國知網(wǎng).附錄：桿子影子隨時間曲線圖：%%太陽赤角求解clc,clear;%清空環(huán)境Ts=(31-21+1)+30+31+30+31+31+30+21;%計算lmd經(jīng)度lmd=180+90*((Ts-186)/90);lmd=pi*lmd/180;%將角度轉化為弧度%計算直射點緯度角度K正弦值sink=0.39775*sin(lmd);k=asin(sink);%輸出太陽直射緯度k%緯度化為弧度fai=39+(54/60)+(26/3600);fai=pi*fai/180;%%時角問題求解x=9:0.1:15;fai1=116+(23/60)+(29/3600);st=x+(fai1-120)/15w=15*(st-12);%計算太陽高度w=w/180*pi;y=sin(fai)*sin(k)+cos(fai)*cos(k)*cos(w);z=asin(y);%輸出太陽高度角l=3./tan(z);x=9:0.1:15;figure(1);plot(x,y,'-b')%太陽角變化xlabel('時間點');ylabel('太陽角度變化');figure(2);plot(x,l,'r-')%繪制桿長變化圖像xlabel('時間點');ylabel('直桿的太陽影子長度');clc;clear;as=0.39795*cos(0.98563*(295-173));a=asin(as)b=0.7302;t=9:0.01:15;w=pi/12*(t-12);c=sin(b)*sin(a)+cos(b)*cos(a)*cos(w);l=asin(c);plot(t,l);%太陽高度角；s=3./tan(l);plot(t,s)xlabel('t');ylabel('L');title('直桿的太陽影子長度隨時間的變化曲線');利用附件1曲線擬合并用regress進行檢驗:clearclct=xlsread('附件1-3','D4:D24');x=xlsread('附件1-3','B4:B24');y=xlsread('附件1-3','C4:C24');l=sqrt(x.*x+y.^2);p=polyfit(t,l,2)x1=linspace(9,15);f1=polyval(p,x1);plot(x1,f1)holdonk=polyder(p)scatter(t,l,5)A=polyfit(t,l,2)z=polyval(A,t);[b,bint,r,rint,stats]=regress(l,[ones(21,1),t,t.^2]);statsplot(t,l,'k+',t,z,'r')xlabel('t');ylabel('L');title('直桿的太陽影子長度隨時間的變化曲線');利用附件2曲線擬合并用regress進行檢驗:clearclct=((12+42./60):0.05:(13+42./60))';x=xlsread('附件1-3',2,'B4:B24');y=xlsread('附件1-3',2,'C4:C24');l=sqrt(x.*x+y.^2);p=polyfit(t,l,2)x1=linspace(6,18);f1=polyval(p,x1);plot(x1,f1)holdonk=polyder(p)scatter(t,l,5)A=polyfit(t,l,2)z=polyval(A,t);[b,bint,r,rint,stats]=regress(l,[ones(21,1),t,t.^2]);statsplot(t,l,'k+',t,z,'r')xlabel('t');ylabel('L');title('直桿的太陽影子長度隨時間的變化曲線');利用附件2曲線擬合并用regress進行檢驗:clearclct=((12+42./60):0.05:(13+42./60))';x=xlsread('附件1-3',2,'B4:B24');y=xlsread('附件1-3',2,'C4:C24');l=sqrt(x.*x+y.^2);p=polyfit(t,l,2)x1=linspace(6,18);f1=polyval(p,x1);plot(x1,f1)holdonk=polyder(p)scatter(t,l,5)A=polyfit(t,l,2)z=polyval(A,t);[b,bint,r,rint,stats]=regress(l,[ones(21,1),t,t.^2]);statsplot(t,l,'k+',t,z,'r')xlabel('t');ylabel('L');title('直桿的太陽影子長度隨時間的變化曲線');緯度計算公式：利用附件3曲線擬合并用regre