摘要通過太陽影子定位技術可以確定視頻的拍攝地點和時間，為拍攝出更好的視頻，掌握太陽影子的變化規律就變得尤為重要。本文主要綜合運用了地理學、幾何學、統計學、數學分析和高等代數等知識，并利用MATLAB,SPSS和mathematica等計算機軟件，通過建立數學模型來研究影子長度的變化特征,進一步確定視頻的拍攝地點和時間。針對問題一，首先我們通過分析影子長度的影響因素得到與影子長度的關系（見表達式六）整理計算之后，就得到了影子長度的數學模型。然后我們通過分析他們之間的關系，再利用MATLAB編程，得到了影子長度關于各個參數的變化規律（見圖3到圖7）。其次根據我們建立的模型，利用MATLAB編程畫出了給定時間天安門廣場3米高的直桿的太陽影子長度的變化曲線（見圖8），然后在考慮折射率的情況下又畫了一條變化曲線（見圖9），最后進行了誤差分析（見圖10）。針對問題二，我們采用了測試分析法（數據分析法和計算機仿真相結合），通過分析各個參量之間的關系，先以影長為目標做回歸，用模型一的模型，通過SPSS進行擬合得到多組數據，再用MATLAB進行檢驗得到符合的兩組經緯度。然后我們又以太陽方位角為目標做回歸,得到模型（見表達式12），其計算方法與影長做回歸目標時一樣。我們分步做了兩次擬合，先用MATLAB擬合出經度，再做回歸模型（見表達式14）最后得到經緯度和桿長。綜上可知，肯定有一地點是在海南，還有一個地點可能在云南。針對問題三，我們用問題二中的多項式回歸，得到回歸模型（見表達式17和20）利用附件二得到的經緯度為和桿長，得到天數。利用附件三得到的經緯度為和桿長，得到天數針對問題四，首先運用MATLAB軟件，根據畫面灰度，運用MATLAB軟件，把視頻轉化成二值圖，求得影子端點的像素坐標，然后根據相似原理，把像素坐標轉化成水平面上的坐標（消去了視角的影響），進而求得影子的長度。用以上方法求得的數據，運用多次擬合的方法，得到該地的經緯度為，日期未知時，得到的經緯度與其相似?！娟P鍵字】影子長度多項式擬合太陽方位角畫面灰度問題重述如何確定視頻的拍攝地點和拍攝日期是視頻數據分析的重要方面，太陽影子定位技術就是通過分析視頻中物體的太陽影子變化，確定視頻拍攝的地點和日期的一種方法。1.建立影子長度變化的數學模型，分析影子長度關于各個參數的變化規律，并應用建立的模型畫出已知時間天安門廣場3米高的直桿的太陽影子長度的變化曲線。2.根據某固定直桿在水平地面上的太陽影子頂點坐標數據，建立數學模型確定直桿所處的地點。將你們的模型應用于附件1的影子頂點坐標數據，給出若干個可能的地點。3.根據某固定直桿在水平地面上的太陽影子頂點坐標數據，建立數學模型確定直桿所處的地點和日期。將你們的模型分別應用于附件2和附件3的影子頂點坐標數據，給出若干個可能的地點與日期。4．先利用軟件提取視頻中的數據，再根據數據改善模型，求出若干個可能的拍攝點。當拍攝日期未知時，確定出拍攝地點與日期。問題分析由題可知，本題具體分析如下：問題一：本題要求建立影子長度變化的數學模型，這需要我們給出影子長度變化的影響因素。查閱文獻了解到各個參量的定義及其表達式，然后聯立即可得到影子長度變化的數學模型。分析影子長度關于各個參數的變化規律，首先我們要在保持其他參數不變的情況下，只改變一個參數，來研究影子長度的變化規律。對于具體問題的變化曲線，因為參數的值已經給出，帶入模型，利用軟件編程就可以畫出它的變化曲線。問題二：本題要求根據某固定直桿在水平地面上的太陽影子頂點坐標數據，建立數學模型確定直桿所處的地點，然后將模型應用于附件1的影子頂點坐標數據，給出若干個可能的地點。先用附件一中給的數據即頂點的與坐標，計算出影子的長度,然后用SPSS做回歸擬合，得出的數據再用MATLAB進行檢驗。問題三：問題三是問題二的拓展，建立數學模型確定直桿所處的地點和日期，比問題二多了一個未知量，我們可以采用問題二的模型和計算方法來解決本問題。問題四：問題四是前兩問的具體應用，只要求出視頻中影子的長度就可以運用前面的模型求解。對于求取視頻中的影長，可以用MATLAB軟件編寫程序，設定恰當的灰度閥值，把視頻轉化為二值圖像。從圖片右下角開始掃描桿子和方塊（桿子底座），求得影子右端、桿子底部的坐標。由于是在三維空間中拍攝的，圖片中物體的長寬比與實際的長寬比不同，可以根據桿子底盤的長寬比求得物體實際的長寬比。根據相似度原理，由桿子的實際長度、圖片中的像素維度等，求得像素與實際長度的比。最后，用影子右端、桿子底部的坐標、物體長寬比、像素與實際長度的比，求出影子的長度。得到上述數據之后，應用MATLAB進行多項式擬合，和應用SPSS軟件進行非線性回歸，兩次擬合得到經緯度和日期?；炯僭O一年是365天；地球表面是球表面；地球的公轉是正圓；大氣層有折射無厚度；視頻中桿的底盤是正方形，不考慮厚度；桿沒有厚度。四、主要變量的符號說明為了便于描述問題，本文將問題中涉及的主要變量用下列符號來表示(如下表1所示)，有些變量將在文中用到時陸續說明。符號代表的含義符號代表的含義桿長影子的長度太陽時太陽赤緯角當地的緯度時角太陽的高度角折射后的太陽高度角日期序號折射率太陽方位角 原始像素的高度球坐標系切平面的坐標系表1五、模型的建立與求解5.1問題一建模和求解5.1.1.影子長度變化的數學模型（1）太陽時：太陽時是指以太陽日為標準來計算的時間，時間的計量以地球自轉為依據，地球自轉一周，計24太陽時，當太陽達到正南處為12:00。太陽時可以分為真太陽時和平太陽時，平太陽時假設地球繞太陽是標準的圓形，鐘表所指的時稱為平太陽時（簡稱為平時），我國采用東經120度經圈上的平太陽時作為全國的標準時間，即“北京時間”。（2）時角：時角是以正午12點為0度開始算，每一小時為15度，上午為負下午為正，即10點和14點分別為-30度和30度（因計算需要，把度數換算成了弧度）。因此，時角的計算公式為(弧度）(1)(2)其中為太陽時，是當地時間與北京時間的時間差（單位：小時）。赤緯角:赤緯角也稱為太陽赤緯，即太陽直射緯度，其計算公式近似為(3)其中為日期序號，例如，1月1日為，10月22日為。太陽高度角：指某地太陽光線與該地作垂直于地心的地表切線的夾角。由圖1可知，三點的坐標為由此可得：即太陽高度角的計算公式(4)其中為太陽高度角，為時角，為當時的太陽赤緯，為當地的緯度(北京的緯度為北緯)。圖1由圖2可知，即（5）圖2（6）（6）將公式（1）（2）（3）帶入（4），得出的值之后再代入（5）就能得到影子長度變化的數學模型（7）5.1.2影子長度關于各個參數的變化規律根據建立的數學模型，利用MATLAB編程（見附錄），通過改變各個參數的值，觀察影子長度的變化規律。首先，我們保持其他的變量不變，改變的值，由圖3我們可以看到每天的最短影子長度在一年內呈周期性變化。圖3（2）接下來我們保持其他的變量不變，改變經度的值，我們選取了東經和東經。由圖4可以知道東經的度數越小，影子最短的時刻越往后。圖4（3）其次我們在保持其他的變量不變，改變緯度的值，我們選取了北緯和北緯。由圖5,圖6可以知道北緯的度數越大，影子越長。圖5圖6（4）同樣的，我們依然保持其他的變量不變，改變時間，由圖7可以知道影子最短的時刻在中午附近，越遠離正午，影子越長。圖7（5）最后，我們依然保持其他的變量不變，改變桿長，由可知桿長與影子的長度成正比，即桿長越長，影子越長。5.1.3應用模型求解具體問題根據我們建立的模型，利用MATLAB編程畫出2015年10月22日北京時間9:00-15:00之間天安門廣場3米高的直桿的太陽影子長度的變化曲線如圖8所示。圖8然后我們又考慮了折射率，根據（8），得到的變化曲線如圖9所示圖9由圖8和圖9對比可以知道，是否考慮折射率差別小于千分之三，所以畫了誤差圖10以便更好的觀察誤差。圖10折射率還與大氣層的厚度有關，但考慮到公式的復雜程度，這里不作考慮。5.2.問題二建模和求解5.2.1用一次擬合做回歸5.2.1.1以長度為目標做回歸根據問題一建立的數學模型（7）用附件一中給的數據即頂點的與坐標，計算出影子的長度,然后根據計算得到日期序號，然后利用SPSS以影長為因變量，不斷給定桿長、緯度和一組初值做多次回歸擬合(表達式見附錄)，得到多組數據（部分數據在附錄中給出），然后我們利用MATLAB進行檢驗。最后得到兩組符合的擬合數據,經緯度表示就是經查閱地圖可知這兩個地點分別在海南島上和云南省內。由于沒有考慮太陽方位角的問題，分不清方向，所以我們將會在接下來做改進。5.2.1.2以太陽方位角為目標做回歸太陽方位角是太陽在方位上的角度，它通常被定義為從北方沿著地平線順時針量度的角。由圖11可知太陽向量和法線的坐標為且（9）是正午的直射點，是桿的位置，是直射點，指向正東，指向正北，指向正上。在原坐標系上表示的是太陽的方向，切平面地表水平面（我們將其拿出來進行單獨分析由于坐標系，過渡矩陣是在平面上的坐標。又因為圖11由此可以得到然后結合各參量的表達式得到數學模型利用SPSS以為因變量不斷給定桿長、緯度和一組初值做多次回歸擬合,雖然我們考慮了太陽的方位角問題，但是畫出的圖形仍然存在誤差。5.2.2用多次擬合做回歸5.2.2.1確定該地的經度首先我們根據附錄一中所給的數據用MATLAB畫出散點圖，由高度角、方位角和天文知識可以知道擬合的多項式為偶函數（關于正午時刻對稱），為了得到更加精確的結果，我們采用了四次多項式模擬，然后用多項式擬合（程序見附錄）多項式為（13），模擬出全天（從早上九點到下午五點）桿的影子長度圖12然后利用mathematica求出曲線的根和頂點，得出最短影子長度出現的時間，進而確定該地的經度。5.2.2.2以影子長度為目標做回歸根據上面求出的經度，就可以確定時角的值，然后利用SPSS以影子長度為因變量不斷給定桿長和緯度初值做多次回歸擬合，擬合的回歸方程為（14）由以上公式得其中，可以把代入上式。由于我們版本的SPSS沒有函數，所以我們將上面的回歸方程化為（15）最后得到經緯度和桿長.綜上可知，肯定有一地點是在海南，桿也可能在云南。5.3.問題三建模和求解5.3.1.1確定經度問題三是問題二的拓展，與第二問相比，該題的日期是未知量，即日期序號是未知量。所以這個問題我們也可用多次擬合得到結果。首先我們根據附錄二中所給的數據用MATLAB畫出散點圖，我們依舊可以推得用偶次多項式進行擬合。然后用MATLAB求解得到經度。擬合的多項式為（16）模擬出全天（從早上九點到下午五點）桿的影子長度圖13然后利用mathematica求出曲線的根和頂點，得出最短影子長度出現的時間,進而確定該地的經度為。5.3.1.2以影子長度為目標做回歸根據上面求的經度，就可以確定時角的值，然后利用SPSS以影子長度為因變量不斷給定參變量桿長和緯度初值做多次回歸擬合，擬合的回歸方程為（17)轉化為SPSS用的方程為：（18）最后我們得到的經緯度為和桿長，得到天數即為5月20號，由對稱性可知7月23號也是該題的答案，然后根據經緯度地圖可以得到在新疆省內。5.3.2.1確定經度首先我們根據附錄三中所給的數據用MATLAB畫出散點圖，用偶次多項式進行擬合。然后用MATLAB求解得到經度，擬合多項式為(19)模擬出全天（從早上九點到下午五點）桿的影子長度圖14然后利用mathematica求出曲線的根和頂點，得出最短影子長度出現的時間進而確定該地的經度為5.3.2.2以影子長度為目標做回歸根據上面求的經度，就可以確定時角的值，然后利用SPSS以影子長度為因變量不斷給定參變量桿長和緯度初值做多次回歸擬合，擬合的回歸方程為(20)轉化為SPSS用的方程為：(21)最后我們得到的經緯度為和桿長，得到天數即為11月3號，同樣由對稱性可知2月8號也是該問題答案。然后根據經緯度地圖可以得到在湖北省內。5.4.問題四建模和求解5.4.1讀取視頻數據（1）讀取圖片運用MATLAB中的VideoReader函數可以知道視頻共有40分鐘41秒，每秒25幀，圖片像素是1920×1080，影子變化不是很明顯。為了分析結果明顯，每隔2分鐘取一次圖片，以桿子和影子為圖片中心，截取圖片的像素為961×749，并設定恰當灰度閥值，將彩色圖片轉化為二值圖像（程序見“將視頻變為圖片.m”），一共得到21張二值圖像，其中的一張圖片如下圖：圖15（2）求取影子端點像素坐標由圖15可知，影子端點大致分布在圖片的右下角，運用MATLAB軟件編程（程序見“視頻中影子端點的像素坐標.m”），從右下角開始掃描，尋找影子最右端，得到桿子的端點坐標如下圖：圖片序號1234567891011X’軸704704704706708708709710710712714Y’軸877869857849837827821809798790778圖片序號12131415161718192021X’軸716714716717718718719719719719Y’軸773763755741734725717707690691表2由于桿子底座近似方形的，以桿子底部為其中心點，運用MATLAB軟件編程（程序見“視頻處理（桿子底部像素坐標）.m”），在圖片中尋找方塊，求其中心點，得出桿子底部的的像素坐標為。（3）求取坐標系下，影子端點的坐標以二值圖片中桿子的底端為坐標原點，圖片像素矩陣的縱列為軸，建立坐標系（單位：像素）,坐標圖如下圖：影子影子Y’軸X’軸O圖16通過計算可以得出影子右端的坐標，如下表3：圖片序號1234567891011軸-11-11-11-9-7-7-6-5-5-3-1軸808800788780768758752740729721709圖片序號12131415161718192021軸1-112334444軸704694686672665656648638621622表3（4）求取XOY坐標系下，影子端點的坐標通過測量，得到原始圖片測量高度為，原始圖片中桿的測量高度為。根據相似原理，可得（22）其中原始圖片中桿的像素高為，原始圖片的像素高度。通過計算，可以得出像素。由題給的桿長，可得像素比，即圖片中每337像素為1m長。以桿子的底端為坐標原點，正對鏡頭方向為X軸正方向，建立水平面XOY坐標系（單位：米），坐標圖如下圖：Y軸Y軸X軸O影子圖17對于Y軸方向，物體在XOY坐標系中Y軸坐標、物體在坐標系中軸坐標和像素比的關系如下：由于圖片來自三維空間，攝影機拍攝角度的不同，各物體在圖片上的長度關系也會顯示不同。通過測量，可以得出桿子底座長寬比=長/寬=4.29。對于X軸方向，物體在坐標系中X軸坐標、物體在坐標系中軸坐標的關系如下：（23）通過計算可以得出水平面XOY坐標系上影子右端的坐標，見表4：圖片序號1234567軸-0.13989-0.13989-0.13989-0.11446-0.08902-0.08902-0.07630軸2.397632.373892.338282.314542.278932.249262.23145圖片序號891011121314軸-0.06359-0.06359-0.03815-0.012720.01272-0.012720.01272軸2.195852.163202.139472.103862.089022.059352.03561圖片序號15161718192021軸0.025430.038150.038150.050870.050870.050870.05087軸1.994071.973291.946591.922851.893181.842731.84570表4（5）求取影子長度通過影子右端的坐標可以計算出影子的長度，（24）影子長度（見附錄--視頻中影子長度）的散點圖如下圖：圖185.4.2模型的建立與求解5.4.2.1經度的確立由以上方法得到時間，影子的長度數據，為了減少誤差，我們根據散點圖，排除兩個異常點，并且為了數據的穩定性，我們對每兩個相鄰影子長度數據進行了取平均數處理。我們用上述處理過的數據用MATLAB畫出散點圖，用偶次多項式進行擬合。然后用MATLAB求解得到經度。擬合的多項式為（25）模擬出全天（從早上九點到下午五點）桿的影子長度圖19然后利用mathematica求出曲線的根和頂點，得出最短影子長度出現的時間,進而確定該地的經度為5.4.2.2以影子長度為目標做回歸根據上面求的經度，就可以確定時角的值，然后利用SPSS以影子長度為因變量不斷給定緯度初值做多次回歸擬合，擬合的回歸方程為（26）由其中得到且的值代入上式。轉化為SPSS用的回歸方程為：（27）最后我們得到的經緯度為，然后根據經緯度地圖可以得到在陜西省內。5.4.2.3當日期未知時的求解當日期未知時我們可以確定視頻的拍攝地點和日期。我們依舊用多次擬合的方法來求解。首先用上述的數據和MATLAB擬合的方法得到經度，因為MATLAB擬合的方法不涉及到日期，所以經度不變，還是，時間也不變。根據上面得到經度，然后利用SPSS以影子長度為因變量不斷給定緯度和日期序數初值做多次回歸擬合，擬合的回歸方程為（28）其中。轉化為SPSS用的回歸方程為：（29）最后我們得到的經緯度為，日期序號即6月22日。然后根據經緯度地圖可以得到依舊在陜西省內。和第一次沒有給日期的擬合結果相比較，地點基本相同，日期相差了22天。六、模型評價6.1模型的優點（1）在問題一的模型中，我們考慮了太陽的折射率，減小了影子長度的誤差；（2）給出了在水平面上太陽的實際方向的變化函數；（3）給出了計算并描繪影子軌跡的MATLAB程序；（4）在問題二的求解中用了多種回歸方法，互相驗證；(5)在問題二和問題三當中，運用一次回歸和二次回歸，得到相似的結果。在問題四中，用MATLAB程序處理視頻，使得誤差較小。在問題四中，有無日期，可以得到相似的結果，說明了模型的有效性。6.2模型的不足：提取視頻中的數據不夠精確，忽略了底盤的高度6.3模型的改進:對于問題二的模型，我們沒有考慮折射率的問題，接下來我們建立模型將折射率考慮在內。將影子長度分解到與軸上，得到（30）由表達式（n=1.000277）（7）可以得到（31）所以得到下面的數學模型（32）參考文獻[1]伍光和，自然地理學，北京：高等教育出版社，2000。[2]呂林根許子道，解析幾何，北京：高等教育出版社，2006.5。[3]茆詩松程依明濮曉龍，北京：高等教育出版社，2011.2。[4]華東師范大學數學系，數學分析上（下），北京：高等教育出版社，2010.7。[5]北京大學數學系幾何與代數教研室前代數小組，高等代數，北京：高等教育出版社，2003.9。附錄附錄1：問題一全年最短影子變化曲線phi=.6965;pi=3.1415926;t0=12+0.2406;date=[295:365+294+1];L0=3;t=12;omega=(t-t0).*15/360*2*pi;theta=0.4092797.*sin(2*pi*(284+date)/365);w=asin(cos(phi).*cos(theta).*cos(omega)+sin(phi).*sin(theta));l=L0./tan(w);%%A=asin(cos(theta).*sin(omega)./cos(w));y1=cos(theta)*cos(omega)*sin(phi)-cos(phi)*sin(theta);x1=cos(theta)*sin(omega);xy=sqrt(x1.^2+y1.^2);x=x1./xy.*l;y=-y1./xy.*l;plot(date,l)xlabel('最短影子長度');ylabel('最短影子長度');title('全年最短影子變化曲線');gridon;附錄2：問題一影子長度關于時間的變化曲線(沒考慮折射率)phi=0.696512;pi=3.1415926;t0=12.240574;date=295;t=[9:0.1:15];omega=(t-t0)*15/360*2*pi;theta=0.4092797*sin(2*pi*(284+date)/365);w=asin(cos(phi)*cos(theta)*cos(omega)+sin(phi)*sin(theta));l=3./tan(w);plot(t,l)xlabel('時間');ylabel('影子長度');title('太陽影子長度的變化曲線');gridon;附錄3：問題一考慮折射率后，影子長度關于時間的變化曲線

phi=0.696512;pi=3.1415926;t0=12.240574;date=295;n=1.000277

t=[9:0.1:15];

omega=(t-t0)*15/360*2*pi;

theta=0.4092797*sin(2*pi*(284+date)/365);

w=asin(cos(phi)*cos(theta)*cos(omega)+sin(phi)*sin(theta));

l=3./tan(w);

plot(t,l)

xlabel('時間');

ylabel('影子長度');

title('影子長度變化曲線(考慮了折射率)');

grid

on;附錄4：問題一影子長度變化誤差圖

phi=0.696512;pi=3.1415926;t0=12.240574;date=295;n=1.000277

t=[9:0.1:15];

omega=(t-t0)*15/360*2*pi;

theta=0.4092797*sin(2*pi*(284+date)/365);

w=asin(cos(phi)*cos(theta)*cos(omega)+sin(phi)*sin(theta));

w1=pi/2-asin(sin(pi/2-w)/n)

l=3./tan(w);

l1=3./tan(w1);

plot(t,l)

xlabel('時間');

ylabel('未考慮折射率和考慮了折射率的影子差值');

title('誤差圖');

grid

on;附錄5：問題二以桿長的影子做回歸SPSS用的表達式l*SQRT(1/((COS(N)*0.98321997*COS(0.2618*(t-12-t1))+SIN(N)*0.18242391)**2)-1)附錄6：問題二以桿長的影子做回歸SPSS得到的部分數據賦予初值擬合的結果Nt1lNt1l0.3-30.56.356-10.8402.2500.3-20.51.415-7.2550.0450.3-10.5-0.0731.1602.2570.300.5-0.0731.1602.2570.310.50.3350.7812.0500.320.50.3350.7812.0500.4-20.50.3350.7812.0500.5-20.51.415-7.2550.0450.6-20.50.3350.7812.2501-230.3350.7812.0500.4-130.3350.7812.0500.5-13-0.0731.1602.2571030.4191.4392.7290.713-0.0731.1602.2571130.4191.4392.7290.6-33-2.8070.7812.0500.9-33-2.8070.7812.0501.3-110.3350.7812.0501.2110.3350.7812.0501.21.520.3350.7812.0501.2121.7274.7450.04500.10.50.4191.4392.7290.60.220.3350.7812.0500.7-11-0.0731.1602.2570.8230.4191.4392.7290.813-0.0731.1602.2571130.4191.4392.729012-0.0731.1602.2571.40.121.415-7.2550.450附錄7：問題二以太陽方位角做回歸SPSS用的表達式0.98321997*SIN(0.2618*(t-12-t1))/SQRT(1-(COS(N)*0.98321997*COS(0.2618*(t-12-t1))+SIN(N)*0.18242391)**2)附錄8：問題二多次擬合做回歸求經度，用MATLAB模擬出全天（從早上九點到下午五點）桿的影子長度%%用多項式擬合，y0=[1.149625826 1.182198976 1.215296955 1.249051052 1.28319534 1.317993149 1.353364049 1.389387091 1.426152856 1.463399853 1.501481622 1.540231817 1.579853316 1.620144515 1.661270613 1.703290633 1.74620591 1.790050915 1.835014272 1.880875001 1.927918447];x0=[14.7 14.75 14.8 14.85 14.9 14.95 15 15.05 15.1 15.15 15.2 15.25 15.3 15.35 15.4 15.45 15.5 15.55 15.6 15.65 15.7];plot(y0,x0,'*');x=9:0.01:17;p1=polyfit(x0,y0,4);%p1為影子長度y關于時間x的二項式函數y1=polyval(p1,x);plot(x,y1)附錄9：問題二的多次擬合做回歸用SPSS以影子長度為目標做回歸的表達式l*SQRT(1/((COS(N)*0.98321997*COS(0.2618*(t–12.71))+SIN(N)*0.18242391)**2)-1)附錄10：問題三附錄二中所給的數據用MATLAB求經度y0=[1.247256205 1.22279459 1.198921486 1.175428964 1.152439573 1.12991747 1.10783548 1.086254206 1.065081072 1.044446265 1.024264126 1.004640314 0.985490908 0.966790494 0.948584735 0.930927881 0.91375175 0.897109051 0.880973762 0.865492259 0.850504468];x0=[12.68333333 12.73333333 12.78333333 12.83333333 12.88333333 12.93333333 12.98333333 13.03333333 13.08333333 13.13333333 13.18333333 13.23333333 13.28333333 13.33333333 13.38333333 13.43333333 13.48333333 13.53333333 13.58333333 13.63333333 13.68333333];plot(y0,x0,'*');x=9:0.001:17;p1=polyfit(x0,y0,4);%p1為影子長度y關于時間x的二項式函數y1=polyval(p1,x);plot(x,y1)附錄11：問題三附錄二用多次擬合做回歸，用SPSS以影子長度為目標做回歸的表達式l*SQRT(1/((COS(N)*COS(0.4028*SIN(2*3.1415*(284+c)/365))*COS(0.2618*(t-14.7))+SIN(N)*SIN(0.4028*SIN(2*3.1415*(284+c)/365)))**2)-1)附錄12：問題三附錄三中所給的數據用MATLAB求經度%%用多項式擬合y0=[3.533142184 3.546768029 3.561797643 3.578100715 3.595750783 3.61493428 3.635425983 3.657218272 3.680541115 3.705167836 3.731278025 3.758917911 3.788087888 3.818701015 3.850809619 3.88458522 3.919911828 3.956875992 3.99553479 4.035750835 4.077863059];x0=[13.15 13.2 13.25 13.3 13.35 13.4 13.45 13.5 13.55 13.6 13.65 13.7 13.75 13.8 13.85 13.9 13.95 14 14.05 14.1 14.15];plot(y0,x0,'*');x=9:0.001:17;p1=polyfit(x0,y0,4);%p1為影子長度y關于時間x的二項式函數y1=polyval(p1,x);plot(x,y1)附錄13：問題三附錄三用多次擬合做回歸，用SPSS以影子長度為目標做回歸的表達式l*SQRT(1/((COS(N)*COS(0.4028*SIN(2*3.1415*(284+c)/365))*COS(0.2618*(t-12.65))+SIN(N)*SIN(0.4028*SIN(2*3.1415*(284+c)/365)))**2)-1)附錄14：問題四將視頻變為圖片的MATLAB程序

[filename,pathname,fileindex]=uigetfile('*.avi','選擇視頻文件','video.avi','Multiselect','on');

if

ischar(filename)

%只有選擇了文件才進行以下計算

video=VideoReader([pathname

filename]);

%%General

Settings:

%Duration

時間

%Name

視頻名稱

%Path

路徑

%Tag

=

VideoReader

%Type

=

mmreader

%UserData

=

[]

%%Video

Settings:

%BitsPerPixel=24

%FrameRate

視頻采集速率

%Height

高度

%NumberOfFrames

總幀數

%VideoFormat

圖像模式

%Width

=

寬度

LEN=video.NumberOfFrames;

%獲得視頻長度

dir=strcat(pathname,strrep(filename,'.avi',''),'\pic');

mkdir(dir);

fn=strrep(filename,'.avi','');

%thresh

=

graythresh(200);

for

k=1:3000:LEN-1

%若read到len，常會報錯如下???

MATLAB:read:readTimedOut，read到len-1就好了

frame=im2bw(read(video,k),0.8);

imwrite(frame(166:914,800:1760),strcat(dir,'\',fn,'-avi-',int2str(k),'.bmp'),'bmp');%把每幀圖像存入硬盤

end

elseif

iscell(filename)

navi=length(filename);

for

n=1:navi

video=VideoReader([pathname

filename{n}]);

%%General

Settings:

%Duration

時間

%Name

視頻名稱

%Path

路徑

%Tag

=

%Type

=

mmreader

%UserData

=

[]

%%Video

Settings:

%BitsPerPixel=24

%FrameRate

視頻采集速率

%Height

高度

%NumberOfFrames

總幀數

%VideoFormat

圖像模式

%Width

=

寬度

LEN=video.NumberOfFrames;

%獲得視頻長度

dir=strcat(pathname,strrep(filename{n},'.avi',''),'\pic');

mkdir(dir);

fn=strrep(filename{n},'.avi','');

for

k=1:3000:LEN-1

%若read到len，常會報錯如下???

MATLAB:read:readTimedOut，獨到len-1就好了

frame=im2bw(read(video,k),0.8);

imwrite(frame(166:914,800:1760),strcat(dir,'\',fn,'-avi-',int2str(k),'.bmp'),'bmp');%把每幀圖像存入硬盤

end

disp(strcat(num2str(n),'/',num2str(navi),'

:

"',filename{n},'"

Finished!',datestr(now,13)));

end

else

return

end

msgbox('所有幀提取完畢,已寫入磁盤！','提示');附錄15：問題四求視頻中影子端點的像素坐標的MATLAB程序%%在運行完“視頻變為圖片.m”后，將圖片的名稱依次變為1,2,3...21.bmp，方便讀取%%把圖片的存儲路徑改為d:\Appendix4\pic%%讀取視頻轉化后圖片中影子端點的像素坐標%%需要多次手動改變圖片的名稱序1,2,3...21.bmpimax=749;jmax=961;i0=0;j0=0a=imread('d:\Appendix4\pic\21.bmp');%imshow(a)%通過顯示圖片，可知影子端點處于右下角p=0;q=0;forj=1:1500;%從圖片右下角開始讀fori=30:70;ifa(imax-i,jmax-j)+a(imax-i,jmax-j-1)++a(imax-i,jmax-j-2)==0%在二值圖中，1表示“黑色”p=i;q=j;break;endendifp>0break;endendi0=imax-pj0=jmax-qfori=1:20forj=1:20a(i0+i,j0+j)=0;endendimshow(a)附錄16：問題四用MATLAB做視頻處理（桿子底部像素坐標）

%%在運行完“視頻變為圖片.m”后，將圖片的名稱依次變為1,2,3...21.bmp，方便讀取

%%把圖片的存儲路徑改為d:\Appendix4\pic

%%需要多次手動改變圖片的名稱序1,2,3...21.bmp

imax=749;jmax=961;isum=0;jsum=0;ii=0;jj=0;

bimax=100;bjmax=140;cc(1:100)=0;dd(1:140)=0;

a=imread('d:\Appendix4\pic\21.bmp');