基于變量控制的太陽影子定位

田 冰，朱英豪，彭楊楊

(1. 華北理工大學數學建模創新實驗室，河北 唐山 063210；2. 華北理工大學冶金與能源學院， 河北 唐山 063210；3. 華北理工大學以升教育創新基地，河北 唐山 063210)

0 引言

影子是一種光學現象。光線在同種介質中沿直線傳播，不能穿過不透明物體而形成的較暗區域，就是影子。太陽影子定位技術逐漸應用于各個領域，其原理為通過觀察物體的影子變化來確定日期時間和物體所在的地理位置。通過建立太陽高度角與影長的物理模型，以及最小二乘法擬合模型，研究出直桿隨各參數的變化規律，利用影子頂點坐標，建立數學模型求解直桿所處位置。利用MATLAB 進行最小二乘法擬合，得到影長關于時間變化的函數方程，利用該方程求解圖像最低點，從而確定觀測點的精確經緯度。

1 太陽影子長度相關函數方程的建立

物體在光線的照射下，會在地面上留下它的投影，這就是影子。物體的影子長短不僅與物體自身體量有關系，還與許多外界因素有關，即位置、日期、時間、太陽高度角等[1]。

對于地球表面上的一個點，太陽高度角指的是太陽光的入射方向與地平面切線的夾角，記作h。隨著地球的自轉運動，太陽光線與地平面切線的夾角改變，太陽高度角也隨之變化。太陽赤緯是地球赤道平面與太陽和地球中心的連線之間的夾角，等于太陽直射點的緯度值，記作σ 。太陽時角時太陽所處的位置與正午太陽位置之間的角度差，記作ε 。地理緯度記作δ 。

在進行模型建立前，首先建立合理的假設條件，簡化模型：

（1）假設影子的長度與地形無關；（2）假設大氣層分布均勻，對陽光折射率不變。

太陽赤緯δ 的計算公式為[2]：

圖1 太陽相關角的圖像 Fig.1 Image of solar correlation Angle

其中N 為積日，即日期在年內的序號，int 表示取整數。

太陽時角ε 的計算公式為[3]：

太陽高度角的計算公式[3]為：

設物體自身高度為l，影子的長度m，由三角函數可知：

由計算公式可以看出物體影子的長度與太陽高度角有著直接的關系，而隨日期變化的太陽赤緯、隨時間和經度變化的時角都影響著太陽高度角，所以物體影子的長短與自身長度、日期和時間、當地經緯度有關。

2 影長隨各參數變化的模型建立

2.1 影子長度隨地理緯度的變化規律

在影響太陽高度角的這些因素中，首先控制日期、時間、經度不變，僅物體所處的地理緯度改變。為便于研究，設定日期為問題1 給出的2015 年10 月22 日，物體長度為3 米，時間選取正午12 點整，代入公式求解得到影子長度隨緯度變化的關系式：

利用MATLAB 繪制出影長隨地理緯度變化的曲線為：

圖2 影長關于緯度的變化曲線圖 Fig.2 Change curve of shadow length with respect to latitude

總結規律為：當物體的日期時間和經度確定，隨著地理緯度的升高，影長逐漸增加。

2.2 影子長度隨時間的變化規律

研究了影子長短隨緯度變化的規律后，控制日期和經緯度不變，研究物體影長關于時間的變化規律。仍以問題1 中2015 年10 月22 日天安門廣場3 米高的直桿為參照，研究北京時間9：00-15：00 的影子長短變化，計算得到影子長度隨時間的變化規律為：

當ε∈[9,15]時，計算出部分影長如下表：

表1 部分時間對應的影長 Tab.1 Corresponding shadow lengths of part time

利用MATLAB 繪制出影長及太陽高度角隨時間變化的曲線為：

圖3 影長及太陽高度角隨時間的變化曲線圖 Fig.3 Curve of shadow length and solar altitude Angle over time

總結規律為：當物體的高度和地點確定，在一天之內，從早上開始影長逐漸減小，減小的速度越來越慢，到正午12 點時影長最短，從正午到傍晚時，影長不斷增加，增加的速度越來越快。

2.3 基于曲線擬合的經度求解

根據某一時刻物體影子的頂點坐標x 和y，可以得出這一時刻影子長度為：

利用附件4 中頂點坐標數據可以計算出14：42 到15：42 這一個小時內不同時刻的影子長度。在問題1 中，最終求得的影長隨時間的變化曲線為拋物線，于是利用這一個小時的影長與時間的數據，在MATLAB 中對影長和時間進行拋物線擬合[4]，得到曲線方程為：

圖像最底點處即為導數值為0 的點，也是當地正午12 點對應的影長，對方程（8）求導可知在北京時間12.6 小時的時候影長取最小值。

圖4 影長隨時間變化的擬合曲線圖 Fig.4 Fitting curve of shadow length changing with time

由地理知識可知，影子最短的時刻對應于一天的正午12 點，此點對應的北京時間為12.6，即該地與北京的時差為0.6 小時。所以直桿所處地點的經度為：

3 基于最小二乘法的緯度迭代求解

確定了物體所在地的經度之后，為了確定其緯度，引入太陽方位角，即直桿在太陽光線下的投影與正南方向的夾角γ，其計算公式[5]為：

根據影子定點坐標數據確定實測量方位角為

在物體所在地的經緯度未知時，時角的值不能確定，若利用問題1 建立的若干太陽角數學模型求解緯度，涉及的未知量太多，無法直接求解方程得出結果。考慮到變量個數不為一且范圍廣，采用尋找最優解的思想，模擬出各種可能的緯度，得到它的方位角，然后與實際值作誤差分析，若誤差足夠小，尋找使模擬出的數據與真實數據為誤差最小值時，即為緯度最優解。

在這里采用最小二乘法對緯度進行優化求解[6]，太陽方位角和實際方位角的相對誤差平方和為：

目標函數為：

約束條件為：

在MATLAB 中進行緯度的迭代[7]，設定初始緯度范圍是[-90°, +90°]，設定步長為20°，計算出滿足條件的大致范圍，在這些符合條件的區域中重新設定步長為0.5°，繼續進行迭代選取出精確到更高的區域后，在這些區域中重復進行迭代，就這樣一直擇優選取迭代區域，使精度越來越高，直到符合各參數精度要求后結束，確定了誤差平方和最小時的緯度為19.50°N 和19.68°N。

所以附件1 影子的數據對應的地點可能為（19.50°N，110°E）和（19.68°N，110°E），在地圖上進行定位后如下圖：

圖5 （19.68°N，110°E）定位圖 Fig.5 (19.68°N, 110°E) location diagram

圖6 （19.50°N，110°E）定位圖 Fig.6 (19.50°N, 110°E) location diagram

4 結論

（1）太陽影子長度變化規律：正午之前，隨著時間的推移，影長逐漸減小；正午以后，隨時間的推移，影長逐漸增加。

（2）在已知影長隨時間變化的圖像近似為拋物線前提下，于是根據已有數據確定影長，并對影長和時間進行拋物線擬合，求得影長最小值處所對應的北京時間為12：36，根據兩地的時差為36 分鐘，求得影子所在地的經度為110°E。以太陽方位角和實際測量方位角的誤差為切入點，用最小二乘法逐步選取出更精確的緯度范圍，進行迭代確定出符合條件的緯度最優解為19.50°N 和19.68°N。

（3）在建立影子長度變化的數學模型中，分析影子長度關于各參數的變化規律，并在建立3m 高的直桿的太陽影子長度變化曲線中，采用了控制變量法，將多因素問題變為單因素求解問題，使要研究的參數對影長的影響更加顯著。

（4）根據某固定直桿在水平面上的太陽影子頂點坐標數據，建立數學模型確定其所處的具體地點和日期的推理過程中，利用曲線擬合思想由已知推未知，并利用迭代求解目標函數最優解，避免了參數不確定性的討論。問題3 中，未知參數過多，在第2 問的基礎上進行多參數、多層迭代，設置的變量范圍廣，求解的結果也越精確。問題4 中，利用圖像分析技術，將圖像中有用的數據提取出來，結合現有的問題3 中的模型直接進行求解，使求解過程更加簡便。

（5）曲線擬合的思想還可以用在科學和工程上，利用采集的離散數據樣本，能都得到一條連續函數或者更加密集的離散樣本。圖像分析中的圖像特征提取能夠運用在分類器上分類紋理圖像，在信號處理上也有重要作用。