何偉怡 李向春 明德烈
(1.華中科技大學自動化學院 武漢 430074)(2.火箭軍研究院 北京 100094)
海面紅外輻射特性與本身的物理特性和環境因素密切相關,其中海面在紅外輻射特性分析下,涉及到的物理特性主要包括海面溫度、海面發射率和反射率[1]。海面溫度和發射率影響著自身紅外輻射能量,海面反射率則影響其對外界太陽輻射和天空背景輻射的反射比例[2],兩者最終影響海面在紅外波段下輻射到探測器的能量值,所以在對海面紅外輻射特性的研究中,這三個參數就顯得至關重要。本文主要介紹海面溫度模型的建立和求解。
海面海水溫度與測量點所處經緯度、季節天氣均有緊密聯系,因而具有多變性。目前,根據不同的要求,測量海水溫度主要有下面幾種方法[3]:一是使用投入式測溫儀進行測量,這種方法主要適用于對特定海域固定深度的多年溫度變化檢測;二是采用遙感技術對海水溫度進行測量,此方法主要用于對大面積海域或者全球海域溫度的檢測。目前,包括中國在內的許多國家已經發射了多顆極軌氣象衛星,為全球海域全方位、多角度觀測提供了便利。
然而,由于海洋極端天氣條件對船只及人員安全性的威脅,以及特定海域局勢軍事敏感性的限制,使用投入式測溫儀對特定海域進行海水溫度的測量就受到較大的限制。基于遙感技術的海溫測量對于反應海水表層24小時內的溫度變化精確度以及時效性有待考量。所以,前述測量方法適用于氣象學研究,但無法為紅外實時仿真所使用。
針對海面紅外仿真實時性高、對海水表層溫度變化精度要求高的特點,需要建立專門的模型計算海水表皮在一天之中隨時間推移的變化,其中一種常用的方法是查閱當地海域的實測資料(海面溫度、大氣溫度、相對濕度和風速等參數),然后對實測數據進行分析,研究水溫、氣溫以及其它氣象參數之間的數學關系,推導包含以上參數的回歸方程和回歸曲線,得到一個計算海面溫度的經驗公式[4]。這種方法的優點在于計算方便,能滿足一些實時性要求較高的系統需求,但是缺點也很明顯,因為經驗公式是通過當地實測數據的分析推導而來,存在一定的誤差,并且公式的局限性較大,如果要計算另一片海域的海面溫度,就需要重新分析實測數據和推導回歸方程。
計算海面溫度的另一種方法就是對海面溫度進行數學建模,即利用熱平衡理論對海面建立熱平衡方程,把海面溫度當作熱平衡方程的未知量,通過求解平衡方程中其它分項的值,最后再計算得到海面溫度。
本文將采用熱平衡方程來計算海面溫度,并針對海面熱平衡方程中的海面熱傳導分量的計算方法進行研究。在求解方程過程中,首先通過隱式差分方法對方程進行空間和時間維度上離散化,接著將方程轉換成矩陣運算形式,最后利用追趕法求解海面溫度。
根據能量守恒定律,對于整個海水表面薄層來說,任意時刻,表面從外界吸收的能量與表面向外界散失的能量滿足以下熱平衡方程[5]:

其中各項含義如下:
R為凈入射輻射通量,即表面從外界(太陽、天空等)獲取的凈輻射通量,方向為從上向下指向表面為正。不考慮外部熱源的情況下,其值等于太陽入射通量和天空入射通量與表面自發輻射散失的輻射通量之差:

G為海面向下熱傳導熱量,即地面表面朝下向更深層熱傳導散失的熱量,方向為從表面向下為正;
H為海面與大氣交換熱量,即表面與大氣發生熱交換散失的能量,方向為從表面向上為正;
LE為海面蒸發熱量,即表面由于水分蒸發作用而散失的熱量,方向為從表面向上為正。
除了G分量,方程剩余分量與傳統方法中的熱平衡方程一致。傳統方法計算得到的G分量表達式較為簡單,不能正確反應出海水表面向下傳導的熱量,得到的結果誤差較大。
本文采用一維傅里葉熱傳導方程來計算海面向下傳導的熱量,求解方程過程中,對方程進行空間和時間維度上進行離散化,通過設定初始值迭代計算得到逐小時海面溫度。得到海面溫度結果較傳統方法更準確,因此計算復雜度也較大,可以通過離線計算和查表方式來解決實時性問題。
接下來具體介紹熱平衡方程中每個分項的具體計算方式。
1)凈入射輻射熱量
凈入射輻射通量包括太陽入射熱量、天空入射熱量和海面自發輻射熱量,其中太陽輻射[6]為短波輻射,到達地球大氣層外的太陽輻射可以通過日地距離修正太陽常數精確計算得到;天空輻射[7]分為大氣輻射和云層輻射,均為長波輻射;物體通過輻射向外散發能量,由于地面物體的溫度不高,自發輻射也主要是長波輻射[7]。以上三種輻射能量的計算比較常用,這里就不詳細展開。
2)海面與大氣交換熱量
海面與大氣之間交換的熱量也稱感熱通量,熱量交換的速度和方向與表面與近地層空氣的溫差以及風速等因素相關。通常采用空氣動力學公式進行計算[8]:

其中,ρa為濕空氣的密度,Cp為空氣的定壓比熱容,Ts為表面溫度,Ta為大氣溫度,ra空氣動力學阻力[9]。
3)海面蒸發熱量
海面蒸發熱量又稱為潛熱通量,與感熱通量不同,該通量為表面水分蒸發作用帶走的熱量,對于海面水分可以采用以下公式[8]:

其中,L為水的汽化潛熱,ρv為絕對濕度。
4)海面向下熱傳導熱量
海面向下熱傳導熱量的計算滿足傅里葉方程,對于海面薄層,按一維熱傳導計算,方程為[10]

其中,T為表面層溫度,λ為材質的導熱系數。特別地,表面層傳導的熱量G為

由于在研究海面溫度計算中,我們忽略洋流、內部漩渦和表面溫度差異等因素,認為海面在某一時刻具有相同的溫度,并且在一定深度海水溫度是保持一個定值[11],那么可以把研究對象縮小成具有一定深度的海水薄層,海洋表面滿足上方給出的表面熱平衡方程,海水薄層底部再往深處的海水則保持一個恒定的溫度,薄層內部則滿足一維傅里葉熱傳導方程,只存在層與層之間的熱量傳遞。
經過上述對海面模型的簡化和分析,將海水薄層沿著垂直方向上進行空間劃分,從海面到薄層底部分成若干層,Zi為每一層的厚度,如果不考慮層與層之間溫度變化速度的差異,可以把每層厚度設成相同的值,本文采用等比數列設置薄層厚度,越往深處薄層厚度越大。
除了進行空間劃分,還需要對海水薄層進行時間上的劃分,由于紅外仿真系統是模擬海面一天24小時的輻射變化,所以時間維度上從0到24小時劃分,步長的△t設成定值。
對于薄層內部,由于我們忽略了洋流等復雜因素的影響,海水只存在層與層之間的能量傳遞,即滿足一維傅里葉熱傳導方程,那么將熱傳導方程從空間和時間維度上進行離散劃分,轉為隱式差分方程[12]:

其中λi為每一層的導熱系數,由于忽略復雜因素的影響,每層海水的導熱系數相同;△Zi表示第i層的厚度,Tit表示在t時刻下的第i層海水的溫度;ρi表示第i層海水的密度;ci表示第i層海水的比熱容;這里每層的密度和比熱容都是相同值。
在公式中t時刻下每層溫度值和每一層厚度、海水導熱系數、密度、比熱容都是已知的,需要推導t+Δt下每一層的溫度,那么將上述差分公式變換成如下形式:

其中經過公式換算,得到A、B、C、D系數表達式如下所示,其中i大于1,且小于n:經過上述公式轉換,可以很方便地將方程式轉換成矩陣運算,如下所示:


由于這里討論的是海水薄層內部的能量交換,所以上述矩陣運算中,從第二層到第n-1層的系數是已知的,第一層和最后一層的A、B、C、D系數仍是未知量,還需要對上下邊界進行單獨推導分析。
海水表面不但要和第二層海水進行熱量交換,還和大氣直接接觸,涉及到大氣背景輻射、風速和相對濕度等因素的影響,所以一維傅里葉熱傳導方程不適合用于計算表面溫度。這里需要結合之前提到的表面熱平衡方程來解算表面溫度場數據,同樣采用隱式差分的方式,對該方程進行空間和時間維度上離散化,得到如下公式:


將一階微分近似后公式代入表面熱平衡方程中,得到如下公式:

和上小節類似,轉換成如下形式:
其中經過公式換算,得到A、B、C、D系數如下所示:
通過對表面熱平衡方程進行隱式差分轉換,得到矩陣運算中系數A、B、C、D的表達式。值得一提是這里系數A值為0,這是由于海面不存在上一層海水,無需與上層海水交換熱量,所以不需要上層海水溫度的分量。
對于海水薄層下邊界而言,除了要和倒數第二層海水進行熱量交換,它還要和薄層底部以下的海水進行熱量交換;這個區別在于,我們把薄層底部以下海水經過復雜的熱量交換,達到穩定的溫度值,而這部分輻射熱量計算較為復雜,需要單獨處理。這里采用方式是將底部能量交換總結為一個換熱系數αin,從而底部熱傳導方程如下所示[12]:

其中,αin為換熱系數;Tin為海水薄層底部以下海水的溫度值,采用經驗數據;Tn為海水薄層底部的溫度值。
同樣由于該熱傳導方程存在微分項,需要對其進行空間和時間維度上離散化,經過隱式差分轉換得到如下公式:

和上小節類似,轉換成如下形式:

其中經過公式換算,得到A、B、C、D系數如下所示:

通過對海水薄層上下邊界和內部熱量傳遞方程的解算,最終將熱傳導方程轉成矩陣運算,公式中系數矩陣和等號右邊的列向量為已知量,公式中第二項的列向量表示為下一Δt時刻每層海水的溫度值,為待求解。通過分析公式中第一項系數矩陣為三對角矩陣,其非零元素集中分布在主對角線和其相鄰兩個次對角線上,我們把這種矩陣公式稱為三對角方程組,可以通過追趕法來進行求解[13]。
在利用追趕法求解該三對角方程組過程中,由于方程式從空間和時間維度上進行離散化,意味著溫度場數據的求解是一個遞推的過程,需要擬定一個初始值才能解算出下一Δt時刻的海水薄層每層的溫度值。本文擬定以下起始條件:時間起點為凌晨4點,此時海面層溫度與大氣溫度相同,海水表層以下每層溫度按線性插值計算得到。設定該起始條件基于以下假設:經過整個夜晚的降溫和熱交換,海水表面的溫度應與大氣溫度相近。
求解過程如圖1所示,圖1左側為初始溫度場數據的設置,選定的時間為凌晨四點,海水表面溫度與大氣溫度相等,然后通過線性插值得到表層以下每層海水的溫度值,接著根據初始值按步長Δt進行迭代,直到解算出仿真時刻的海水表面溫度。
首先本文先進行海面溫度與大氣溫度比較分析。由于海水表面直接接觸著大氣,大氣溫度直接影響著海水溫度,使得大氣溫度和海水溫度有著相同的上升下降趨勢;接著將海面溫度仿真數據和實際測量數據進行對比驗證,但是由于海面溫度逐小時數據難以獲取,比如中國氣象網中需要較高權限才能下載全球海平面溫度資料數據,普通學生用戶無法得到,只能從國家海洋環境預報中心網站中獲取鄰近海域的最低最高溫度數據,所以本文將用仿真計算得到的海面溫度數據與同一海域的最高最低溫度數據進行比較,觀察仿真數據是否分布在最高最低溫度的范圍附近。

圖1 迭代計算海面溫度
本文采用網上氣象數據和海洋溫度數據的時間為2016年8月和2016年1月,分別代表著夏天和冬天季節,海域經緯度位置約為東經118.1°,北緯38.7°。其中海面溫度計算過程中所使用到的參數信息如表1所示。

表1 參數信息

圖2 夏天海溫和氣溫比較
以上兩幅圖分別為夏天和冬天季節下海面溫度仿真數據和氣象網大氣溫度數據的對比圖,從圖中可以看出,凌晨時段海水溫度比大氣溫度高,海水溫度變化比較緩慢,然后從日出時段后,由于大氣的比熱容較小,升溫速度較海面快,大概下午14點到15點大氣溫度達到了峰值,隨后海面溫度也跟著達到峰值,最后隨著日落兩者的溫度也逐漸降低。因此可以看出,海面溫度仿真數據和實測大氣溫度具有大致相同的上升和下降趨勢,說明本文海面溫度計算模型能較好反映大氣溫度和太陽輻射對海面溫度的影響。

圖3 冬天海溫和氣溫比較

圖4 夏天溫度數據對比

圖5 冬天溫度數據對比
從上面兩幅圖可以看出,雖然海面溫度仿真數據與實測海面溫度范圍有一定的誤差,這是因為本文忽略了洋流等復雜因素的影響,在整體上本文的仿真數據還是能很好地分布在最高最低值附近,較好反映實測海面溫度變化的范圍,從而表明本文海面溫度計算模型具有較好的正確性。
本文從海面的熱平衡方程出發,介紹了平衡方程中的各項:太陽輻射、天空的長波輻射、海水的蒸發潛熱、海水對大氣的熱傳導和海水向深處的熱傳導分量,在建立海面熱平衡方程基礎上,利用隱式差分方法將海面模型在空間和時間維度上進行離散化,接著將方程轉換成矩陣運算形式,最后利用追趕法迭代求出仿真時刻的海面溫度數據。并對海面溫度的仿真數據進行對比驗證,分別和歷史大氣溫度數據、海面最高最低溫度數據進行比較,實驗結果表明仿真數據滿足實際溫度變化趨勢,說明本文海面溫度計算模型具有較好實用性,能為海面背景紅外輻射模型提供相關輸入,同時可以為海洋環境的紅外檢測提供參考。