趙珮瑤

摘要：該文針對太陽影子定位問題建立了依賴于拍攝地點(diǎn)經(jīng)度、緯度、日期的影子長度變化的數(shù)學(xué)模型，分析確定拍攝地點(diǎn)的日期。

關(guān)鍵詞：超定方程組；最小二乘；數(shù)據(jù)矯正

中圖分類號：TP311 文獻(xiàn)標(biāo)識碼：A 文章編號：1009-3044（2016）17-0260-03

1 問題重述

1.1 問題背景

在進(jìn)行視頻數(shù)據(jù)分析時(shí)，通過太陽影子定位技術(shù)，即分析視頻中的物體的太陽影子變化，是確定視頻的拍攝地點(diǎn)和拍攝日期的一種重要方法。

1.2 目標(biāo)任務(wù)

問題一：建立影子長度變化的數(shù)學(xué)模型，分析影子長度的變化規(guī)律，并畫出2015年10月22日北京時(shí)間9：00-15：00之間天安門廣場（北緯39度54分26秒，東經(jīng)116度23分29秒）3米高的直桿的太陽影子長度的變化曲線。

問題二：根據(jù)某固定直桿在水平地面上的太陽影子頂點(diǎn)坐標(biāo)數(shù)據(jù)，建立數(shù)學(xué)模型確定直桿所處的地點(diǎn)。并將所建模型應(yīng)用于附件1的影子頂點(diǎn)坐標(biāo)數(shù)據(jù)，給出若干個(gè)可能的地點(diǎn)。

問題三：附件4為一根直桿在太陽下的影子變化的視頻，并且已通過某種方式估計(jì)出直桿的高度為2米。請建立確定視頻拍攝地點(diǎn)的數(shù)學(xué)模型，并應(yīng)用你們的模型給出若干個(gè)可能的拍攝地點(diǎn)。

如果日期已知，該如何處理？

2 模型的建立與求解

2.1 問題一模型建立與求解

2.1.1 問題一模型的建立

我們采用太陽高度角和三角函數(shù)相結(jié)合的思路。直桿與太陽照射下成的影子可構(gòu)成直角三角形，其三角函數(shù)關(guān)系式可以表示為

[sinθ=HH2+L2] ⑴

其中，[H]表示直桿長度，[L]表示直桿影長，[sinθ]表示太陽高度角的正弦值。

由參考文獻(xiàn)[1]，太陽高度角隨著地方時(shí)和太陽赤緯的變化而變化。太陽高度角的計(jì)算公式為：

[sinθ=sinφsinδ+cosφcosδcost] ⑵

上式中，[θ]，。表示太陽高度角，[φ]表示觀測地理緯度（太陽赤緯和地理緯度都是北緯為正，南緯為負(fù)），[δ]表示太陽赤緯（與太陽直射點(diǎn)維度相等），[t]表示時(shí)角。

2.1.2 問題一模型的求解

（1）太陽赤緯[δ]與地理緯度[φ]正弦值的確定

根據(jù)太陽赤緯的公式：

[sinδ=0.39795cos0.98563N-173] ⑶

式中[N]表示觀測日期，自每年1月1日開始計(jì)算。即代入2015年10月22日，[N]=295。

計(jì)算可得：

[sinδ=0.39795cos0.98563295-173=-0.20046] ⑷

因?yàn)樘彀查T廣場地處（北緯39度54分26秒，東經(jīng)116度23分29秒），與代表北京時(shí)間的東八區(qū)經(jīng)度差為0，所以天不考慮安門廣場的經(jīng)度位置，得

[sinφ=sin[（39+54/60+26/3600）]=0.641546] ⑸

（2）太陽高度角的正弦值的求解

由公式⑷、⑸計(jì)算太陽赤緯[σ]與地理緯度[φ]的余弦值如下：

[cosδ=0.979702cosφ=0.767084] ⑹

由參考文獻(xiàn)[3]可知，時(shí)角公式為：

[t=（T-12）×15]° ⑺

我們在計(jì)算時(shí)角時(shí)發(fā)現(xiàn)，北京時(shí)間并不是北京（東經(jīng)116.4°）地方的時(shí)間，而是東經(jīng)120°地方的地方時(shí)間。根據(jù)經(jīng)度變化15°，時(shí)間變化1小時(shí)這一地理規(guī)則，我們推算出北京當(dāng)?shù)氐恼嬲龝r(shí)間比代表東八區(qū)的北京時(shí)間晚約14分鐘。計(jì)算[sinθ]如下表1所示：

（3）計(jì)算直桿的太陽影子的長度變化

結(jié)合公式⑴的三角函數(shù)公式以及表1數(shù)據(jù)，計(jì)算不同時(shí)刻影長如下表2所示：

通過MATLAB擬合曲線，我們可知直桿的太陽影子長度[L]最短的時(shí)間并不是在12時(shí)，而是在11時(shí)至12時(shí)之間。

2.2 問題二模型建立與求解

2.2.1 問題二模型的建立

（1）直桿直角坐標(biāo)系的建立

假設(shè)可能地區(qū)在太陽直射點(diǎn)的以北地區(qū)數(shù)據(jù)為正。2015年4月18日的太陽直射點(diǎn)在赤道和北回歸線之間，且相對靠近赤道的位置。我們假設(shè)，在附表中所給14：42-15：42的時(shí)間段內(nèi)，太陽在直桿的南面。并且在14：42-15；42這一時(shí)間段，已經(jīng)過正午太陽高度角所處的12時(shí)，太陽從直桿的西面進(jìn)行照射。綜上可得到假設(shè)，直桿的影子在它的東北方向。以直桿底端為坐標(biāo)原點(diǎn)，分別以東向、北向?yàn)閇x]軸、[y]軸的方向，建立直角坐標(biāo)系，如圖3所示：

圖1 直桿直角坐標(biāo)系的建立

（2）經(jīng)度定位模型的建立

由參考文獻(xiàn)[5]，太陽方位角與觀測地點(diǎn)經(jīng)度、緯度以及觀測時(shí)刻的關(guān)系如下：

[A=arcsin（cosδsintcosθ）] ⑻

根據(jù)太陽高度角的余弦值為

[cosθ=LH2+L2] ⑼

結(jié)合公式⑻和公式⑼，化簡為

[LsinA=H2+L2cosδsint] ⑽

在上式⑽中，我們可根據(jù)附件一的數(shù)據(jù)求解太陽方位角的正弦值，由表3和公式⑶可知不同時(shí)間的影長和太陽赤緯的余弦值。最終我們得到了未知量為直桿高度[H]、時(shí)角[t]的經(jīng)度的定位模型。

（3）緯度定位模型的建立

結(jié)合問題一中有關(guān)太陽高度角的公式⑴⑵⑺，得到下列方程：

[H=H2+L2C1sinφ+C2cosφcos[（T-12）×15]] ⑾

上式中，[C1]表示太陽赤緯的正弦值，[C2]表示太陽赤緯的余弦值，因?yàn)闀r(shí)間是確定的，所以太陽赤緯的正弦值、余弦值都是常數(shù)。從而我們得到了未知量為直桿高度[H]、地理緯度[φ]的緯度的定位模型。

2.2.2 問題二模型的求解

（1）太陽赤緯的正弦值、余弦值的確定

太陽方位角是指太陽光線在地平線上的投影與當(dāng)?shù)亟?jīng)線的夾角。根據(jù)所給數(shù)據(jù)，由勾股定理[L=x2+y2]及太陽方位角計(jì)算公式：

[A=arctanxy] ⑿

我們計(jì)算出不同時(shí)間的影子長度，見下表3。

（2）利用解超定方程組的最小二乘解思想求解經(jīng)緯度的定位模型

對于經(jīng)度的定位模型，我們所給數(shù)據(jù)的21個(gè)不同組的數(shù)據(jù)，得到21組不同長度的直桿影長，建立了含有21個(gè)不同的方程，只有直桿高度[H]、時(shí)角[t]兩個(gè)未知量。利用[MATLAB]中的[fsolve]函數(shù)，假定直桿高度[H]的初值為0米，時(shí)間初值為上午9時(shí)。計(jì)算得直桿的高度[H1]為1.0134米，推測測量地點(diǎn)的時(shí)間為14時(shí)55分5秒。已知代表東八區(qū)的北京時(shí)間（東經(jīng)120）為14時(shí)42分，對比14時(shí)55分5秒和14時(shí)42分，根據(jù)地區(qū)經(jīng)度與時(shí)間差的地理知識，計(jì)算出該地的經(jīng)度為東經(jīng)1231612。

同樣，對于緯度的定位模型，我們可以根據(jù)附件一的21個(gè)組不同的數(shù)據(jù)，建立了含有21個(gè)不同的方程，只有直桿高度[H]、地理緯度[φ]兩個(gè)未知量。計(jì)算可得，直桿的高度[H3]為1.0440米，地理緯度[φ]的弧度角為-0.1831[rad]。根據(jù)緯度算法，得到以下的計(jì)算公式⒀：

當(dāng)?shù)乩砭暥萚φ]的弧度角為0.3025[rad]時(shí)，將其數(shù)值代入上述公式⑿，算的第一組結(jié)果的緯度值為北緯17°18；當(dāng)?shù)乩砭暥萚φ]的弧度角為-0.1831[rad]，表示為弧度制2.9585[rad]，將其數(shù)值代入公式⒀，算的第二組結(jié)果的緯度值為南緯10°30。

結(jié)合所求經(jīng)度，我們得到符合條件的兩個(gè)地區(qū)：菲律賓的呂宋島附近（北緯17°18，東經(jīng)1231612）和印度尼西亞努沙登加拉群島附近（南緯10°30，東經(jīng)1231612）。

（3）將直桿高度H作為檢驗(yàn)因子進(jìn)行模型檢驗(yàn)

前面我們分別求出直桿的高度[H1]、[H2]、[H3]依次為1.0134米、1.2442米以及1.0440米。[H2]和[H3]的平均值為[E（H）=12（H2+H3）=1.1441]，將[E（H）]記為[H4]。則[H1]與[H4]平均值為[E（H）=12（H1+H4）=1.07875]，則[H1]與[H4]方差為[DH=12i=12（EH2-Hi）2=4.27×10-3]，可知，方差值很小，認(rèn)為兩種方法下計(jì)算的直桿高度近似相等，模型誤差很小。

2.3 問題四模型建立與求解

2.3.1 問題四模型的建立

（1）對視頻進(jìn)行處理

首先利用[KMPlayer]軟件對視頻進(jìn)行取幀，得到一組像素值為[1920×1080]的圖片。因?yàn)橐曨l時(shí)長40分鐘，為保證截取的圖片能準(zhǔn)確反映整個(gè)視頻中直桿影子的變化，我們對視頻一分鐘截取一幀，從8時(shí)56分3秒至9時(shí)34分4秒共截取39張圖片。然后，對這一組圖片進(jìn)行灰度處理，這樣可以使顏色細(xì)節(jié)更加細(xì)致，更容易辨認(rèn)物體邊緣。

（2）對直桿頂端、直桿影子頂端、直桿底端坐標(biāo)的求解

利用[MATLAB]中的[Data Cursor Tools]分別選中你所觀察到的直桿頂端、影子頂端、直桿底端記錄它的像素點(diǎn)坐標(biāo)，我們得到39組灰色矩陣，每一組灰色矩陣中存儲了圖片中的所有像素點(diǎn)的灰度值。以直桿頂端為例，在所記錄的直桿頂端的像素點(diǎn)坐標(biāo)附近尋找灰度值最小的像素點(diǎn)，記錄它的坐標(biāo)，即為直桿頂端的真實(shí)坐標(biāo)點(diǎn)。求解直桿影子頂端、直桿底端的坐標(biāo)點(diǎn)時(shí)也采用相同的方法。但對于直桿底端，我們認(rèn)為底端相鄰兩個(gè)直角的y坐標(biāo)的平均值為直桿底端的坐標(biāo)。已知桿長可估計(jì)為2m，則可根據(jù)直桿頂端、直桿影子頂端、直桿底端的像素坐標(biāo)得到前五幅圖影子的坐標(biāo)如下表所示：

表4 前五幅幀圖中影子的坐標(biāo)

[圖片＼&1＼&2＼&3＼&4＼&5＼&[x]坐標(biāo)（米）＼&2.307463＼&2.298507＼&2.265672＼&2.241791＼&2.229851＼&[y]坐標(biāo)（米）＼&0.016418＼&0.01791＼&0.016418＼&0.016418＼&0.013433＼&]

以直桿底端為原點(diǎn)，以圖片右側(cè)方位為[x]軸，垂直圖片向里的方向?yàn)閇y]軸，建立直角坐標(biāo)系。

圖2 直桿底端坐標(biāo)系的建立

（3）影子[y]軸數(shù)據(jù)加權(quán)確定

因?yàn)閿z影機(jī)在錄制視頻的過程中，與地面存在一個(gè)固定的角度，因?yàn)殓R像成像的原理，[y]軸方面人眼看到的影子長度小于實(shí)際中影子的[y]軸長度。因此引入[y]軸加權(quán)系數(shù)[C]，當(dāng)攝像機(jī)角度不變時(shí)，人眼看到的影子的[y]軸長度與實(shí)際影子的[y]軸長度形成一個(gè)直角三角形的斜邊和直角邊，且比例保持不變，則在對本問題運(yùn)用問題二的數(shù)學(xué)模型進(jìn)行求解時(shí)，[y][坐標(biāo)軸=通過灰色圖像得出的坐標(biāo)值×設(shè)定的權(quán)值C]。由于[x]軸坐標(biāo)與攝像機(jī)觀測方向垂直，[x]軸坐標(biāo)可看作實(shí)際中的[x]軸坐標(biāo)。

2.3.2 問題四模型的求解

當(dāng)日期已知時(shí)：可求出當(dāng)?shù)貢r(shí)間為9時(shí)9分48秒，與之對應(yīng)的北京時(shí)間為8時(shí)56分3秒。[經(jīng)度=116+9+9/60+48/3600-8-56/60-3/360015=119.4o；][北緯=0.5537×180°π=31.7°]，所以該地點(diǎn)北緯314329，東經(jīng)11924，約為江蘇常州。

當(dāng)日期未知時(shí)：可求出當(dāng)?shù)貢r(shí)間為9時(shí)06分45秒。此時(shí)，[N=173]，為2015年6月22日。地處北緯515313，東經(jīng)1184030。

3 模型的評價(jià)與改進(jìn)

問題一中，突出優(yōu)點(diǎn)是考慮到北京當(dāng)?shù)貢r(shí)間和代表東八區(qū)的北京時(shí)間之間存在約14分的時(shí)差，并且應(yīng)用在時(shí)角的計(jì)算中。

問題二中，能夠運(yùn)用解超定方程組的最小二乘法的思想，利用[MATLAB]的[fsolve]函數(shù)求得最優(yōu)解。缺點(diǎn)是在尋求最優(yōu)解的過程中，運(yùn)行程序?qū)Τ踔凳置舾?，?dāng)初值稍作變化時(shí)，所求的數(shù)值變化很大，這點(diǎn)不利于尋求最優(yōu)解。在對此模型進(jìn)行改進(jìn)的過程中，考慮遺傳算法和模擬退火模型，減小初值對結(jié)果的影響。

問題四中，突出的優(yōu)點(diǎn)是考慮到人眼成像的特點(diǎn)，在直桿影子的[y]軸方向引入加權(quán)系數(shù)[C]，使數(shù)據(jù)更加準(zhǔn)確，缺點(diǎn)是運(yùn)行程序的過程中，不同的初值對結(jié)果幾乎沒有影響，觀察不到某些影響因子的作用力。

參考文獻(xiàn)：

[1] 談小生， 葛成輝. 太陽角的計(jì)算方法及其在遙感中的應(yīng)用[J]. 國土資源遙感，1995（6）.

[2] 王昌明. 可照時(shí)數(shù)和太陽高度角計(jì)算公式的簡化證明[J].山東氣象，1989（2）.

[3] 閆利， 胡修兵， 陳長軍， 等. 中波紅外幾何標(biāo)定場靶標(biāo)定制問題研究[J].測繪地理信息，2013（8）.