太陽(yáng)影子定位技術(shù)2015高教社杯數(shù)學(xué)建模獲獎(jiǎng)?wù)撐膮R總

1、太陽(yáng)影子定位技術(shù)摘要本文以太陽(yáng)影子定位技術(shù)為背景，結(jié)合直桿影子軌跡的變化規(guī)律建立數(shù)學(xué)模型。并運(yùn)用視頻數(shù)據(jù)分析的方法，確定拍攝地點(diǎn)及日期等地理信息條件。第一問(wèn)給出了北京時(shí)間、拍攝日期，以及拍攝地點(diǎn)的經(jīng)緯度。我們可以結(jié)合太陽(yáng)赤緯、時(shí)角、直桿的經(jīng)緯度與太陽(yáng)高度角之間的關(guān)系建立模型，求出符合時(shí)間條件要求的太陽(yáng)高度角，再根據(jù)已知的桿的高度和三角公式求出影長(zhǎng)關(guān)于時(shí)間的變化曲線(xiàn)。第二、三問(wèn)在第一問(wèn)的基礎(chǔ)上增加難度，使部分變量未知。通過(guò)文獻(xiàn)查閱和方程推導(dǎo)，得出陰影運(yùn)動(dòng)軌跡形狀是雙曲線(xiàn)的一支，并且具體形狀和當(dāng)?shù)氐木暥纫约俺嗑曈嘘P(guān)，本文根據(jù)這點(diǎn)進(jìn)行模型假設(shè)與建立。附件中給出的坐標(biāo)并不一定是標(biāo)準(zhǔn)地理坐標(biāo)，通過(guò)對(duì)其進(jìn)

2、行坐標(biāo)變換，引入了實(shí)際坐標(biāo)系與標(biāo)準(zhǔn)地理坐標(biāo)系的偏角。在擬合多項(xiàng)高次變量組成的隱函數(shù)方程的過(guò)程中，為增加精確度，運(yùn)用最小二乘法進(jìn)行擬合求解未知參量時(shí)，可以利用直桿陰影頂點(diǎn)軌跡的形狀，建立參量和變量之間的關(guān)系，簡(jiǎn)化需擬合的隱函數(shù)方程。這樣就可以根據(jù)太陽(yáng)影子頂點(diǎn)橫縱坐標(biāo)以及對(duì)應(yīng)的時(shí)刻，把偏角、緯度、經(jīng)度、日期作為未知參數(shù)進(jìn)行擬合，得出要求的地理位置和相應(yīng)的日期。如通過(guò)對(duì)附件1數(shù)據(jù)的擬合求解可得到一組地理坐標(biāo)（東經(jīng)104.425度，北緯15.6578度），對(duì)附件2數(shù)據(jù)的擬合求解可得一個(gè)可能的日期6月21日，坐標(biāo)（東經(jīng)116度，北緯26度），由附件3得到的可能的日期地點(diǎn)為：6月21日，（東經(jīng)164.55

3、度，北緯71.26度）。為了便于定位，根據(jù)一般工程的實(shí)際需求，對(duì)美國(guó)天文學(xué)家紐康（New Comb）提出的太陽(yáng)公式作了綜合、簡(jiǎn)化，舍去了一些高階微小量。結(jié)合測(cè)量學(xué)的理論，用數(shù)學(xué)模型進(jìn)行非線(xiàn)性擬合求得直桿所處的經(jīng)緯度。第四問(wèn)給出一段視頻，實(shí)際是對(duì)前三問(wèn)模型的實(shí)際應(yīng)用。本問(wèn)對(duì)一些已有的論文以及專(zhuān)利進(jìn)行借鑒，創(chuàng)新與簡(jiǎn)化。首先對(duì)視頻中的圖像進(jìn)行取幀，在灰度處理中因?yàn)榧夹g(shù)限制，改為運(yùn)用Matlab二值化處理。并根據(jù)簡(jiǎn)單測(cè)量畫(huà)出運(yùn)行軌跡。運(yùn)用主元分析法求得陰影尖端坐標(biāo)與桿底坐標(biāo)的關(guān)系。確定影子的運(yùn)動(dòng)軌跡。之后借鑒已有成熟理論將2D圖像去畸變，恢復(fù)仿射的度量屬性，通過(guò)對(duì)3D圖形轉(zhuǎn)變2D過(guò)程的逆向推導(dǎo)，將坐標(biāo)

4、恢復(fù)為符合現(xiàn)實(shí)要求的坐標(biāo)。之后回歸前幾問(wèn)建立的的日晷數(shù)學(xué)模型進(jìn)行求解，得到一個(gè)可能的地理坐標(biāo)為（東經(jīng)104.9度，北緯25.33度）。并在最后進(jìn)行誤差修正。關(guān)鍵詞：日晷投影原理、桿影端點(diǎn)軌跡、非線(xiàn)性最小二乘法、主元分析法、 二值化處理、Floodfill圖論算法 一、問(wèn)題重述如何確定視頻的拍攝地點(diǎn)和拍攝日期是視頻數(shù)據(jù)分析的重要方面，太陽(yáng)影子定位技術(shù)就是通過(guò)分析視頻中物體的太陽(yáng)影子變化，確定視頻拍攝的地點(diǎn)和日期的一種方法。1.建立影子長(zhǎng)度變化的數(shù)學(xué)模型，分析影子長(zhǎng)度關(guān)于各個(gè)參數(shù)的變化規(guī)律，并應(yīng)用建立的模型畫(huà)出2015年10月22日北京時(shí)間9:00-15:00之間天安門(mén)廣場(chǎng)（北緯39度54分26秒

5、,東經(jīng)116度23分29秒）3米高的直桿的太陽(yáng)影子長(zhǎng)度的變化曲線(xiàn)。2.根據(jù)某固定直桿在水平地面上的太陽(yáng)影子頂點(diǎn)坐標(biāo)數(shù)據(jù)，建立數(shù)學(xué)模型確定直桿所處的地點(diǎn)。將模型應(yīng)用于附件1的影子頂點(diǎn)坐標(biāo)數(shù)據(jù)，給出若干個(gè)可能的地點(diǎn)。3. 根據(jù)某固定直桿在水平地面上的太陽(yáng)影子頂點(diǎn)坐標(biāo)數(shù)據(jù)，建立數(shù)學(xué)模型確定直桿所處的地點(diǎn)和日期。將模型分別應(yīng)用于附件2和附件3的影子頂點(diǎn)坐標(biāo)數(shù)據(jù)，給出若干個(gè)可能的地點(diǎn)與日期。4附件4為一根直桿在太陽(yáng)下的影子變化的視頻，并且已通過(guò)某種方式估計(jì)出直桿的高度為2米。請(qǐng)建立確定視頻拍攝地點(diǎn)的數(shù)學(xué)模型，并應(yīng)用模型給出若干個(gè)可能的拍攝地點(diǎn)。如果拍攝日期未知，能否根據(jù)視頻確定出拍攝地點(diǎn)與日期？2、 問(wèn)

6、題分析2.1問(wèn)題一的分析本問(wèn)以直桿影子長(zhǎng)度為研究對(duì)象，尋找影響影子長(zhǎng)度與各個(gè)參數(shù)的關(guān)系及其變化規(guī)律。為了使影長(zhǎng)的計(jì)算科學(xué)嚴(yán)謹(jǐn)，我們應(yīng)了解太陽(yáng)與地球之間相互的運(yùn)動(dòng)軌跡，由此計(jì)算出太陽(yáng)對(duì)地球上某一定點(diǎn)的相對(duì)位置。這主要由當(dāng)?shù)氐牡乩砭暥?、季?jié)(月、日)和時(shí)間三個(gè)因素決定，可以用地理緯度()、太陽(yáng)赤緯角()、太陽(yáng)高度角(h)、及時(shí)角t等參數(shù)進(jìn)行定量表達(dá)。2.2問(wèn)題二的分析由于附件中給出的直桿陰影頂點(diǎn)的坐標(biāo)系的x軸和y軸并不一定垂直或重合，有可能坐標(biāo)軸向與南北方向存在一定的偏角（）。本問(wèn)根據(jù)太陽(yáng)影子頂點(diǎn)橫縱坐標(biāo)的21組數(shù)據(jù)，對(duì)x和y的坐標(biāo)進(jìn)行旋轉(zhuǎn)坐標(biāo)變換，得到陰影頂點(diǎn)在新坐標(biāo)系下的坐標(biāo)，該新坐標(biāo)系以正東

7、方向?yàn)閤軸正向，正北方向?yàn)閥軸正向，直桿底端為坐標(biāo)原點(diǎn)，由原坐標(biāo)系旋轉(zhuǎn)得到。這樣就可以由x和y求出相應(yīng)時(shí)刻的太陽(yáng)方位角（A），再結(jié)合太陽(yáng)高度角(h)的計(jì)算公式，經(jīng)過(guò)一系列化簡(jiǎn)，可以得到經(jīng)緯度之間的關(guān)系。將經(jīng)緯度作為參數(shù)，利用Matlab進(jìn)行非線(xiàn)性擬合，選取合適的初值，即可得直桿所處地點(diǎn)的經(jīng)緯度。2.3問(wèn)題三的分析該問(wèn)題的求解可利用問(wèn)題二建立起來(lái)的模型，將由日期確定的太陽(yáng)赤緯作為未知參數(shù)，在Matlab中對(duì)時(shí)間和直桿影子長(zhǎng)度進(jìn)行非線(xiàn)性擬合，選取合適的初值得到經(jīng)緯坐標(biāo)和日期的值。2.4問(wèn)題四的分析本問(wèn)考察基于視頻數(shù)據(jù)分析方法進(jìn)行太陽(yáng)影子定位。從圖像或視頻中估算經(jīng)緯度是目前計(jì)算機(jī)視覺(jué)領(lǐng)域的研究熱點(diǎn)問(wèn)

8、題，估算經(jīng)緯度不僅自身具有重要理論意義，而且它對(duì)計(jì)算機(jī)視覺(jué)問(wèn)題也有積極的啟示意義。第四問(wèn)提供的視頻，其中體現(xiàn)了標(biāo)志物的影子在一段時(shí)間內(nèi)的移動(dòng)軌跡。在這種條件下求經(jīng)緯度，實(shí)際上就是基于視頻中太陽(yáng)影子軌跡來(lái)估計(jì)經(jīng)緯度的實(shí)際應(yīng)用。首先我們把視頻取幀處理得到影子的軌跡點(diǎn)，在擬合出地平線(xiàn)后，運(yùn)用計(jì)算機(jī)作圖的相關(guān)知識(shí)對(duì)坐標(biāo)進(jìn)行糾正后，把圖片上的2D坐標(biāo)恢復(fù)為實(shí)際中真實(shí)的3D坐標(biāo)，最后把問(wèn)題回歸日晷模型算出經(jīng)緯度，并對(duì)因?yàn)榈胤綍r(shí)標(biāo)準(zhǔn)時(shí)引起的誤差進(jìn)行修正。三、模型假設(shè)1.將太陽(yáng)光近似地看成平行光投射到地球；2.忽略太陽(yáng)光受地球大氣層折射和漫反射的影響；3.地球和太陽(yáng)在運(yùn)行中不規(guī)則變化和周期性變化產(chǎn)生的誤差忽略

9、不計(jì)；4.假設(shè)直桿所在的地面為水平的；5.忽略地球形狀對(duì)試驗(yàn)結(jié)果的影響；6.假設(shè)每天的時(shí)間為24小時(shí)整；7.假設(shè)所求日期均均為2015年日期。四、符號(hào)說(shuō)明符號(hào)符號(hào)說(shuō)明物體的地理緯度物體的地理經(jīng)度太陽(yáng)對(duì)應(yīng)日期的赤緯角h太陽(yáng)高度角A太陽(yáng)方位角T時(shí)角H直桿長(zhǎng)度五、模型建立及求解5.1問(wèn)題一模型的建立及求解5.1.1模型的建立(1)影響影長(zhǎng)的因素：其中h為太陽(yáng)光線(xiàn)與水平地面的夾角，即太陽(yáng)高度角；通過(guò)查閱目前國(guó)內(nèi)的大部分天文學(xué)文獻(xiàn)，可得太陽(yáng)高度角的計(jì)算公式1，即：(2) 其中T為時(shí)角，可用計(jì)算，t為24小時(shí)制的當(dāng)?shù)貢r(shí)間;為赤緯角，其較精確公式為2：(3)式中稱(chēng)日角，即(4)這里n又分兩部分組成，即n=N

10、-N0；式中N為積日，即日期在年內(nèi)的順序號(hào)；N0的計(jì)算公式如下：(5)其中INT(X)為取不大于X的最大整數(shù)。5.1.2模型的求解 根據(jù)問(wèn)題一的題目已知，10月22日是2015年的第295天，即N=295。對(duì)于日照計(jì)算來(lái)說(shuō)，地球和太陽(yáng)在運(yùn)行中不規(guī)則變化和周期性變化產(chǎn)生誤差的數(shù)值相當(dāng)小，可以忽略不計(jì)。所以此赤緯角計(jì)算公式的結(jié)果較為精確，可以滿(mǎn)足計(jì)算影長(zhǎng)變化曲線(xiàn)精度要求。由此可建立影長(zhǎng)和北京時(shí)間的數(shù)學(xué)函數(shù)，并作出圖形。Matlab程序見(jiàn)附錄一，影子長(zhǎng)度變化曲線(xiàn)見(jiàn)圖1。圖1 直桿的太陽(yáng)影子長(zhǎng)度的變化曲線(xiàn)5.2問(wèn)題二模型的建立及求解5.2.1模型的建立（1）坐標(biāo)變換第二問(wèn)的附件雖然在直桿所在的地平面處

11、建立了正交分解的x-y軸，但并沒(méi)有指明坐標(biāo)軸和地理正北方向的夾角。據(jù)此分析，可設(shè)為地平面坐標(biāo)系的x軸正向與正東方向的夾角。對(duì)x和y進(jìn)行下式給出的變換，可得直桿頂點(diǎn)在新坐標(biāo)系下的投影坐標(biāo)。變換如下：(6)（2）桿影端點(diǎn)移動(dòng)軌跡 關(guān)于桿影端點(diǎn)移動(dòng)軌跡的周年變化規(guī)律如圖2所示，該圖描繪了一年內(nèi)不同日期的軌跡。夏至日和冬至日是太陽(yáng)直射點(diǎn)移動(dòng)方向發(fā)生轉(zhuǎn)換的日期，春分日和秋分日是直射點(diǎn)的南北半球位置轉(zhuǎn)換的日期。因此，選擇的這幾個(gè)日期具有較強(qiáng)的代表性。在圖中，AA、CC、EE依次表示冬至日、兩分日、夏至日的軌跡，BB、DD表示介于二分二至日的軌跡（如立冬和立夏），該圖清晰的反映除兩分日之外，其余日期的軌跡圖

12、都是雙曲線(xiàn)中的一條。圖2 桿影端點(diǎn)移動(dòng)軌跡的周年變化規(guī)律而對(duì)于具體的某一天中直桿陰影頂點(diǎn)的運(yùn)動(dòng)軌跡的描繪，則可以借助Analemmatic日晷的模型來(lái)描繪。在這種日晷的模型中，日晷采用了垂直的指時(shí)針，它的時(shí)間線(xiàn)是赤道日晷的時(shí)間線(xiàn)在地平面的投影水平方向指示東西軸，垂直方向指向南北軸，投影日晷的位置沿著短軸移動(dòng)，如圖3所示。因此具有時(shí)間一致對(duì)應(yīng)關(guān)系的陰影軌跡如果投影到水平平面上將得到一個(gè)二次曲線(xiàn)，二次曲線(xiàn)的對(duì)稱(chēng)軸位于南北方向指示軸上。軌跡的形狀與當(dāng)?shù)氐木暥纫约疤?yáng)直射點(diǎn)的緯度即太陽(yáng)赤緯有關(guān)，而太陽(yáng)赤緯又與日期相關(guān)，故軌跡的形狀與當(dāng)?shù)氐木暥群腿掌谙嚓P(guān)聯(lián)。圖3 日晷模型及影子變化規(guī)律變化規(guī)律（3）太陽(yáng)

13、方位角和太陽(yáng)高度角有了（1）得到的新坐標(biāo)，就可以用來(lái)表示太陽(yáng)方位角的正切值。由天文學(xué)的相關(guān)知識(shí)3，可以得到以下等式：(7)A為太陽(yáng)方位角，并且(8)又因太陽(yáng)高度角公式為(9)以及(10)聯(lián)立化簡(jiǎn)可得到如下等式：(11)其中(12)(13)5.2.2模型的求解根據(jù)表格附件一做出桿影的軌跡圖如下圖4 附件一桿影頂點(diǎn)軌跡圖對(duì)上述結(jié)果進(jìn)一步化簡(jiǎn)可得：(14)其中(15)上式建立了陰影頂點(diǎn)坐標(biāo)與時(shí)間（t）、經(jīng)度（）、緯度（）以及所用坐標(biāo)系與地理坐標(biāo)系之間偏角（）之間的關(guān)系模型。由于附件已給出直桿陰影頂點(diǎn)在每一時(shí)刻的實(shí)際坐標(biāo)，可以將其它待求的量作為未知參量，基于非線(xiàn)性最小二乘法在Matlab中對(duì)Y和t進(jìn)行

14、擬合，結(jié)合地理資料選取合適的初值求得直桿可能的地點(diǎn)。我們得到了（104.4249E，15.65784N）、（33.0454E，16.7577N）等地點(diǎn)。5.3問(wèn)題三模型的建立與求解5.3.1模型的建立基于問(wèn)題二的模型，我們很自然地得出問(wèn)題三的模型：(16)其中為太陽(yáng)赤緯，是與日期有關(guān)的量，問(wèn)題三相對(duì)于問(wèn)題二的變化是日期作為參量是未知數(shù)。結(jié)合問(wèn)題一中赤緯的計(jì)算公式，可用日期N表示出。這樣就得到了含緯度、經(jīng)度、日期N、偏角四個(gè)未知參量的隱函數(shù)方程，也就是問(wèn)題三的模型。5.3.2模型的求解類(lèi)似于問(wèn)題二的求解，先畫(huà)出附件二、三中桿影端點(diǎn)軌跡的變化規(guī)律如下圖： 圖5 附件二桿影端點(diǎn)軌跡圖圖6 附件三桿影

15、端點(diǎn)軌跡圖關(guān)于時(shí)間的已知參量在第三問(wèn)中變?yōu)槲粗獏⒘?，因此可以在Matlab中調(diào)用非線(xiàn)性擬合函數(shù)擬合隱函數(shù)方程來(lái)求解緯度、經(jīng)度、日期N、偏角四個(gè)未知參量。對(duì)于附件二日期和地理坐標(biāo)的求解，通過(guò)選取合適的不同初值，我們得到了以下不同的地點(diǎn)和日期，如下所示：（1）6月1日，（21N,101E）（2）6月21日，（30N,122E）（3）6月21日，（26N,116E）其中括號(hào)內(nèi)為該地點(diǎn)的經(jīng)緯度。對(duì)于附件三，我們也求解出了一系列點(diǎn)，如下所示：（1）6月21日，（71.26N，164.55E）（2）6月21日，（52.85N，140.11E）5.4問(wèn)題四模型的建立與求解5.4.1模型的建立問(wèn)題四實(shí)質(zhì)與問(wèn)題

16、三的原理相同，都是給定時(shí)刻和該時(shí)刻下現(xiàn)實(shí)坐標(biāo)系中的橫縱坐標(biāo)，求解日期與地理坐標(biāo)；不同點(diǎn)在于本題的坐標(biāo)并不是直接給出的，需要借助視頻處理技術(shù)從中提取出桿影端點(diǎn)的坐標(biāo)。因此，基于對(duì)問(wèn)題二、三的分析可以建立問(wèn)題四的模型，即：(17)模型與問(wèn)題三相同。5.4.2模型的求解（1）從視頻中提取桿影端點(diǎn)坐標(biāo)觀看時(shí)長(zhǎng)為40分鐘的視頻錄像，根據(jù)直桿的高度為2m，并結(jié)合圖中桿長(zhǎng)和影長(zhǎng)的比例關(guān)系，即可求得實(shí)際的桿影的長(zhǎng)度隨時(shí)間的變化情況。但眾所周知，廣角鏡頭所產(chǎn)生的圖片或視頻都存在較為明顯的透視畸變，即被攝體離鏡頭越遠(yuǎn)，在屏面的成像越小。為了避免透視畸變對(duì)實(shí)際影長(zhǎng)測(cè)量帶來(lái)的影響，引入主元分析法解決桿影尖點(diǎn)的確定問(wèn)題

17、。本題模型引用一種半自動(dòng)輸入陰影軌跡檢測(cè)的方法，首先以1分鐘為單元對(duì)所給視頻做取幀處理，形成一組圖片 并構(gòu)建背景圖片B：(18)在B中，每一個(gè)像素點(diǎn)（x,y）是V中最亮的像素點(diǎn)，且是基于灰度級(jí)的。設(shè)置背景絕對(duì)差值，陰影點(diǎn)就可以顯示出來(lái)。然后在陰影區(qū)域運(yùn)用floodfill圖論算法進(jìn)行陰影點(diǎn)計(jì)算，最后的突出陰影點(diǎn)用主成分技術(shù)（PCA），我們同時(shí)運(yùn)用MATLAB對(duì)圖像進(jìn)行二值化處理獲得效果如圖8所示。圖8 二值化圖求相對(duì)坐標(biāo)以直桿底端為原點(diǎn)建立正交的x-y坐標(biāo)系，單位為像素點(diǎn)。結(jié)合主元分析法（PCA）合理得到準(zhǔn)確的陰影運(yùn)動(dòng)軌跡。（2）相機(jī)模型，消隱點(diǎn)的確定以及對(duì)圖像的糾正處理世界坐標(biāo)中陰影點(diǎn)的位置

18、取決于投影物所在緯度位置以及太陽(yáng)光朝向與投影平面之間的幾何關(guān)系。第四問(wèn)進(jìn)行的緯度估計(jì)技術(shù)是基于計(jì)算機(jī)視圖的研究算法，采用幾何分析算法來(lái)進(jìn)行視頻地理位置的估算。由于存在透視畸變，照片中的景物會(huì)和實(shí)際景物有很大的差別。我們生活中有許多透視變化的例子：照片中的景物輪廓形狀改變，但透視變化下只有直線(xiàn)被保留。所以，本問(wèn)的數(shù)學(xué)模型我們?cè)跉W式幾何中加入一些無(wú)窮遠(yuǎn)的理想點(diǎn)形成透視空間，定義一個(gè)平面透視幾何變換作為任何平面中點(diǎn)以保持直線(xiàn)性的映射。在進(jìn)行地理經(jīng)緯度定位中，視頻和相機(jī)校準(zhǔn)是十分重要的一步，要進(jìn)行平面但因性變換。在透視相機(jī)的作用下，有些幾何屬性是保持的，例如共線(xiàn)性，即一條直線(xiàn)透視變換后仍為一條直線(xiàn)，然

19、而一般的透視情況下平行直線(xiàn)將6不再保持平行。透視幾何模型決定了透視效果同時(shí)也提供了相應(yīng)計(jì)算的數(shù)學(xué)表達(dá)式。在計(jì)算機(jī)視覺(jué)中，視覺(jué)幾何是用來(lái)研究在各種變化下仍然保持不變的屬性。從這個(gè)角度來(lái)分析，2D的透視幾何就是來(lái)研究在2D平面下經(jīng)過(guò)各種變換之后仍然保持不變的屬性特征。這種2D平面變換對(duì)于點(diǎn)或者直線(xiàn)來(lái)說(shuō)是不可逆的。另外，透視變換的逆變換也是一種透視變換，所以這種變換同時(shí)存在兩種變換：共平面變換和透視變換，也叫做平面單應(yīng)性。為了理解這種變換，可以假設(shè)用一個(gè)三維向量來(lái)表示一個(gè)二維空間的點(diǎn)x，那么Hx就是這個(gè)點(diǎn)經(jīng)過(guò)平面線(xiàn)性變換后得到的齊次坐標(biāo)。這種變換法則認(rèn)為任何透視效果都會(huì)引起齊次坐標(biāo)的線(xiàn)性變換，并且逆

20、變換也是這種線(xiàn)性透視。一個(gè)平面透視變換是一個(gè)采用三維向量的線(xiàn)性變換，可以用一個(gè)33的矩陣表示：(19)一個(gè)透視變換能把每一個(gè)照片投影成為一個(gè)透視等價(jià)的照片，并且保留所有的不變屬性。由此可見(jiàn)，經(jīng)過(guò)相機(jī)中心的直線(xiàn)定義了一個(gè)平面與另一個(gè)平面的變換關(guān)系。實(shí)際上，如果兩個(gè)坐標(biāo)系統(tǒng)都是歐式空間系統(tǒng)，那么這種由中心透視變換決定的映射就會(huì)比任意的透射變換要產(chǎn)生更多的限制，而這就叫做透視變換而不是完全的投影變換。2D平面的等度量變換能夠保留歐氏距離不變，一個(gè)等度量變換可以表示成：(20)相比之下，給定一個(gè)投影變換面積的等級(jí)大小變換隨著位置的不同而發(fā)生變化。例如在透視變化下,相同平面中一個(gè)較遠(yuǎn)的正方形會(huì)比較近的正

21、方形具有更小的透視結(jié)果。并且經(jīng)過(guò)投影變換的直線(xiàn)的朝向不僅依賴(lài)于原始直線(xiàn)的朝向而且依賴(lài)于它的位置。在本文實(shí)現(xiàn)平面單應(yīng)性計(jì)算的過(guò)程中,需要一個(gè)重要的步驟：對(duì)投影變換進(jìn)行分解。一個(gè)投影變化可以被分解成一系列的變化,其中每一個(gè)矩陣都比前一個(gè)具有更高級(jí)的變換。(20)其中A是一個(gè)非奇異矩陣，K是一個(gè)上三角矩陣并且它的行列式為1。這種分解只有在v不等于的時(shí)候才成立，并且當(dāng)s被確定才會(huì)有唯一的分解結(jié)果。分解出來(lái)的三個(gè)矩陣代表了相應(yīng)類(lèi)型變換的實(shí)質(zhì)變化?？紤]到從透視照片中進(jìn)行校準(zhǔn)，HP將直線(xiàn)恢復(fù)到無(wú)窮遠(yuǎn)；HA影響了仿射屬性，但是沒(méi)有將直線(xiàn)恢復(fù)到無(wú)窮遠(yuǎn)；最后HS代表一般相似變換，并且不會(huì)影響仿射或者透視屬性。從照

22、片中進(jìn)行透視糾正的目的是消除透視照片中的投影畸變，以便得到相似屬性，即角度和長(zhǎng)度比值等可以直接從照片平面中測(cè)量。投影畸變可以通過(guò)照片平面中的四組對(duì)應(yīng)點(diǎn)來(lái)消除，并且能夠明確的計(jì)算出參考點(diǎn)與相應(yīng)的像點(diǎn)之間的映射。而這實(shí)際上過(guò)度具體化了幾何關(guān)系，這是因?yàn)橐粋€(gè)投影變換與相似變換比較起來(lái)只有四個(gè)自由度，所以?xún)H僅需要指定四個(gè)自由度而非八個(gè)就可以確定度量屬性了。在投影幾何中，這四個(gè)自由度給出了幾何物體的物理形狀：消隱線(xiàn)提供了兩個(gè)自由度，兩個(gè)虛圓點(diǎn)也提供了兩個(gè)自由度。因此，當(dāng)消隱點(diǎn)直線(xiàn)的像被確定了，那么透視畸變就可以被消除了，同時(shí)如果虛圓點(diǎn)被確定了，仿射畸變也會(huì)被消除。因此，唯一的畸變就是相似變換了。在投影變

23、換下，一個(gè)理想點(diǎn)被映射到一個(gè)有限的點(diǎn)，因此消隱線(xiàn)被映射到一條有限的直線(xiàn)。但是，如果這種變換是仿射的，那么消隱線(xiàn)不會(huì)被映射為有限的直線(xiàn)而是仍然保持在無(wú)窮遠(yuǎn)的位置。這種變換的逆變換也是成立的，即如果仿射變換是保持無(wú)窮遠(yuǎn)直線(xiàn)的最一般的線(xiàn)性變換。但是消隱線(xiàn)上的點(diǎn)對(duì)于仿射變換不是點(diǎn)對(duì)點(diǎn)的保持。也就是說(shuō)消隱點(diǎn)經(jīng)過(guò)仿射變換后得到的點(diǎn)雖然位于消隱線(xiàn)上，但是卻不是原來(lái)的點(diǎn)。一旦照片平面中無(wú)窮遠(yuǎn)處的直線(xiàn)被確定了，那么就有可能實(shí)現(xiàn)原始平面上的仿射測(cè)量。例如，原始平面上的平行直線(xiàn)可以被識(shí)別為平行，如果這些直線(xiàn)相交于無(wú)窮遠(yuǎn)直線(xiàn)的話(huà)。這是因?yàn)闅W式平面中的平行線(xiàn)相交于無(wú)窮遠(yuǎn)的直線(xiàn)，在經(jīng)過(guò)仿射變換后這些直線(xiàn)仍然相交于無(wú)窮遠(yuǎn)直

24、線(xiàn)的像，因?yàn)橥敢曌儞Q下的相交是保持的。類(lèi)似的，一旦無(wú)窮遠(yuǎn)的直線(xiàn)被確定了，一條直線(xiàn)上的長(zhǎng)度比值就可以通過(guò)三個(gè)點(diǎn)確定的交比來(lái)確定。但是，稍微彎曲的直線(xiàn)卻是更適合這種計(jì)算上的算法，這是因?yàn)閮H僅需要簡(jiǎn)單地將消隱線(xiàn)變換到它的原始位置。實(shí)現(xiàn)了這種變換的投影矩陣可以應(yīng)用到任何點(diǎn)中來(lái)實(shí)現(xiàn)仿射糾正，也就是經(jīng)過(guò)了這種變換，仿射測(cè)量可以直接從校正后的照片中獲得。根據(jù)以上分析編寫(xiě)程序，對(duì)圖像處理得到視頻對(duì)應(yīng)的三維空間內(nèi)的數(shù)據(jù)，并用Matlab擬合，結(jié)果如圖10。圖10 視頻中桿影端點(diǎn)軌跡變化的俯視圖這樣，我們就從視頻中提取出了一組點(diǎn)如下，由于篇幅原因，在此只羅列其中的某幾個(gè)點(diǎn)，詳細(xì)數(shù)據(jù)在附錄六中。時(shí)間/h8.928.

25、938.958.978.989.019.029.03x值804795786780780772766760y值-15-15-15-15-13-13-12-12（3）編程求解求解過(guò)程與問(wèn)題三的求解過(guò)程類(lèi)似，賦初值后得到一個(gè)可能的拍攝地點(diǎn)為（25.33N，104.90E）詳細(xì)的程序代碼見(jiàn)附錄六。5.5問(wèn)題五 此問(wèn)題與問(wèn)題三相同，可將日期作為未知參量，通過(guò)擬合得到日期，詳細(xì)過(guò)程類(lèi)似于問(wèn)題三。6、 誤差分析及敏感度分析將通過(guò)數(shù)學(xué)模型求解的數(shù)值與使用經(jīng)緯儀和有水準(zhǔn)器的支架所實(shí)際觀測(cè)到的直桿影長(zhǎng)相比較。本模型方法的誤差約為，對(duì)50m內(nèi)的物高的測(cè)量精度可達(dá)厘米級(jí)要求。由于擬合的函數(shù)過(guò)于復(fù)雜，參量與變量之間并不

26、能用顯性的關(guān)系表示出來(lái)，這就給對(duì)模型進(jìn)行的靈敏性分析造成困難，因此，我們通過(guò)定性的分析來(lái)說(shuō)明某一參量的變化對(duì)其它量造成的影響。在范圍和步長(zhǎng)給定的情況下，下圖為遍歷所有可能的初始點(diǎn)得到的擬合出的參數(shù)的分布情況。從圖中可以看出，擬合得到的點(diǎn)大概率地分布于某一特定區(qū)域，這說(shuō)明建立的模型較為合理。圖11 不同初值下的所求參數(shù)的分布 七、模型的評(píng)價(jià)及改進(jìn)6.1視頻數(shù)據(jù)分析的優(yōu)點(diǎn)本文提供的太陽(yáng)影子定位的數(shù)學(xué)模型，對(duì)于光照研究以及航空攝影、精密水準(zhǔn)測(cè)量的最佳時(shí)間段的選擇有一定的意義。此數(shù)學(xué)模型還可以推廣到生活中高層建筑群的合理布局和農(nóng)林間種的最佳距離及林帶走向的優(yōu)化設(shè)計(jì)中。第四問(wèn)建立的視頻經(jīng)緯度分析模型也有

27、廣泛的應(yīng)用。即使粗略的經(jīng)緯度估計(jì)也能夠提供有用的線(xiàn)索，預(yù)測(cè)當(dāng)?shù)貧鉁?、平均降雨量等大量背景信息?.2模型求解的不足 由于求解參數(shù)采用的是曲線(xiàn)擬合的方法，其對(duì)初值的依賴(lài)性較大，初值設(shè)置的不合理就有可能造成結(jié)果誤差較大，導(dǎo)致精確性較差。八、參考文獻(xiàn)1 陳曉勇,鄭科科.對(duì)建筑日照計(jì)算中太陽(yáng)赤緯角公式的探討,浙江建筑,第28卷,2011.2建筑設(shè)計(jì)資料集編委會(huì).建筑設(shè)計(jì)資料集M.2版.北京:中國(guó)建筑工業(yè)出版社,1994:179-185.3 林根石,利用太陽(yáng)視坐標(biāo)的計(jì)算進(jìn)行物高測(cè)量與定位,南京林業(yè)大學(xué)學(xué)報(bào),第15卷,1991.4吳濟(jì)廉,影端軌跡周年變化的實(shí)踐與分析以北溫帶地區(qū)為例.地理教學(xué),2013年第

28、10期.5吳濟(jì)廉.關(guān)于桿影端點(diǎn)移動(dòng)軌跡的討論與求證.中學(xué)地理,2011年第3期.6天津大學(xué).基于視頻中太陽(yáng)影子軌跡的經(jīng)緯度估計(jì)方法:中國(guó),200910067817.2012-04-11.九、附錄附錄一：Year=2015;N=295; %N為積日t0=9:0.1:15;T=(t0-12-0.24)*15*pi/180; %T為時(shí)角fai=39.9072*pi/180; %fai為物體的地理緯度H=3;N0=79.6764+0.2422*(Year-1985)-fix(Year-1985)/4);t=N-N0;c=2*pi*t/365.2422; %c為日角delta=0.3723+23.256

29、7*sin(c)+0.1149*sin(2*c)-0.1712*sin(3*c)-0.758*cos(c)+0.3656*cos(2*c)+0.0201*cos(3*c);delta0=delta*pi/180;sinh=sin(fai)* sin(delta0)+cos(fai)*cos(delta0)*cos(T);L=H./tan(asin (sinh);plot(t0,L);title(影長(zhǎng)變化曲線(xiàn))xlabel(北京時(shí)間t)ylabel(影長(zhǎng))附錄二：%本函數(shù)求解赤緯角function delta = chiwei(N,Year)N0=79.6764+0.2422*(Year-198

30、5)-fix(Year-1985)/4);%求解積日a=N-N0;c=2*pi*a/365.2422;delta=0.3723+23.2567*sin(c)+0.1149*sin(2*c)-0.1712*sin(3*c)-0.758*cos(c)+0.3656*cos(2*c)+0.0201*cos(3*c); %求解赤緯角end附錄三：附件一求解程序：clear;%clc a=10.6306*pi/180;F=(p,x)(1-x(:,2).*tan(p(1)./(x(:,2)+tan(p(1)-sin(p(2)*pi/180).*cot(p(3)/15-x(:,1)-20)*15*pi/18

31、0)+tan(chiwei(108,2015)*pi/180)*cos(p(2)*pi/180)./sin(p(3)/15+x(:,1)-20)*15*pi/180);x=14.7 0.59093101814.75 0.57443063214.8 0.55867608914.85 0.5435737514.9 0.52917945214.95 0.51536085615 0.50208868915.05 0.48933541415.1 0.47703894715.15 0.46524911615.2 0.45383972315.25 0.4428299315.3 0.4321740915.35

32、 0.42188874215.4 0.41193212215.45 0.40226832515.5 0.39289045715.55 0.38379141515.6 0.3749206915.65 0.36631837215.7 0.357917966;p0=8 36 115%賦初值warning off p=nlinfit(x,zeros(size(x,1),1),F,p0); disp(num2str(p); plot(x(:,1),x(:,2),ro);hold on; ezplot(x,y)F(p,x,y),0,1,-1e-3,1e-3); title(擬合曲線(xiàn));legend(樣本點(diǎn)

33、)附錄四：附件二求解程序：clear;clc F=(p,x)(1-x(:,2).*tan(p(1)./(x(:,2)+tan(p(1)-sin(p(2)*pi/180).*cot(p(3)/15+x(:,1)-20)*15*pi/180)+tan(chiwei(p(4),2015)*pi/180)*cos(p(2)*pi/180)./sin(p(3)/15+x(:,1)-20)*15*pi/180);x=12.68 -773 -6.39206349212.78 -5.76806640612.83 -5.24103495212.88 -4.7882003412.93 -

34、4.39840319413.98 -4.05427817613.03 -3.75125089313.08 -3.48197278913.13 -3.2402597413.18 -3.02175264113.23 -2.82349433513.28 -2.64248853213.33 -2.47720364713.38 -2.32497331913.43 -2.1836945313.48 -2.05323669113.53 -1.9318622713.58 -1.81794630213.63 -1.71163536413.68 -1.611730598;p0=5 40 116 250%賦初值wa

35、rning off p=nlinfit(x,zeros(size(x,1),1),F,p0)plot(p(3),p(2),.);hold ondisp(num2str(p); plot(x(:,1),x(:,2),ro);hold on; ezplot(x,y)F(p,x,y),0,1,-1e-3,1e-3); title(擬合曲線(xiàn));legend(樣本點(diǎn))附錄五：附件三求解程序：clc F=(p,x)(1-x(:,2).*tan(p(1)./(x(:,2)+tan(p(1)-sin(p(2)*pi/180).*cot(p(3)/15+x(:,1)-20)*15*pi/180)+tan(chi

36、wei(p(4),2015)*pi/180)*cos(p(2)*pi/180)./sin(p(3)/15+x(:,1)-20)*15*pi/180);x=13.15 0.34883093513.2 0.3667377413.25 0.38478430913.3 0.40294684813.35 0.42128134713.4 0.43975703413.45 0.45836359213.5 0.47717150913.55 0.49610263613.6 0.51522603813.65 0.53450633613.7 0.55396733713.75 0.57360236213.8 0.59

37、342265513.85 0.61345357113.9 0.63365129713.95 0.65407712214 0.67470503914.05 0.69555806214.1 0.71661992414.15 0.737931245;p0=5 26 115 222;%賦初值warning off p=nlinfit(x,zeros(size(x,1),1),F,p0); disp(num2str(p); plot(x(:,1),x(:,2),ro);hold on; ezplot(x,y)F(p,x,y),0,1,-1e-3,1e-3); title(擬合曲線(xiàn));legend(樣本點(diǎn)

38、)附錄六：?jiǎn)栴}四求解程序：clear;clc F=(p,x)(1-x(:,2).*tan(p(1)./(x(:,2)+tan(p(1)-sin(p(2)*pi/180).*cot(p(3)/15+x(:,1)-20)*15*pi/180)+tan(chiwei(p(4),2015)*pi/180)*cos(p(2)*pi/180)./sin(p(3)/15+x(:,1)-20)*15*pi/180);x=8.9 8.916666667 8.933333333 8.95 8.966666667 8.983333333 9 9.016666667 9.033333333 9.05 9.066666667 9.083333333 9.1 9.116666667 9.133333333 9.15 9.166666667 9.2 9.216666667 9.233333333 9.283333 9.333333 9.42 9.5 9.57-1.552118549 -1.552141775 -1.551930641 -1.551714674 -1.554131203 -1.554131203 -1.555253537 -1.555131812 -1.555008165 -1.55760448 -1.557499238 -1.558764822 -1.55870