何偉怡 李向春 明德烈

（1.華中科技大學(xué)自動化學(xué)院 武漢 430074）（2.火箭軍研究院 北京 100094）

1 引言

海面紅外輻射特性與本身的物理特性和環(huán)境因素密切相關(guān)，其中海面在紅外輻射特性分析下，涉及到的物理特性主要包括海面溫度、海面發(fā)射率和反射率［1］。海面溫度和發(fā)射率影響著自身紅外輻射能量，海面反射率則影響其對外界太陽輻射和天空背景輻射的反射比例［2］，兩者最終影響海面在紅外波段下輻射到探測器的能量值，所以在對海面紅外輻射特性的研究中，這三個參數(shù)就顯得至關(guān)重要。本文主要介紹海面溫度模型的建立和求解。

海面海水溫度與測量點(diǎn)所處經(jīng)緯度、季節(jié)天氣均有緊密聯(lián)系，因而具有多變性。目前，根據(jù)不同的要求，測量海水溫度主要有下面幾種方法［3］：一是使用投入式測溫儀進(jìn)行測量，這種方法主要適用于對特定海域固定深度的多年溫度變化檢測；二是采用遙感技術(shù)對海水溫度進(jìn)行測量，此方法主要用于對大面積海域或者全球海域溫度的檢測。目前，包括中國在內(nèi)的許多國家已經(jīng)發(fā)射了多顆極軌氣象衛(wèi)星，為全球海域全方位、多角度觀測提供了便利。

然而，由于海洋極端天氣條件對船只及人員安全性的威脅，以及特定海域局勢軍事敏感性的限制，使用投入式測溫儀對特定海域進(jìn)行海水溫度的測量就受到較大的限制?；谶b感技術(shù)的海溫測量對于反應(yīng)海水表層24小時內(nèi)的溫度變化精確度以及時效性有待考量。所以，前述測量方法適用于氣象學(xué)研究，但無法為紅外實(shí)時仿真所使用。

針對海面紅外仿真實(shí)時性高、對海水表層溫度變化精度要求高的特點(diǎn)，需要建立專門的模型計(jì)算海水表皮在一天之中隨時間推移的變化，其中一種常用的方法是查閱當(dāng)?shù)睾Ｓ虻膶?shí)測資料（海面溫度、大氣溫度、相對濕度和風(fēng)速等參數(shù)），然后對實(shí)測數(shù)據(jù)進(jìn)行分析，研究水溫、氣溫以及其它氣象參數(shù)之間的數(shù)學(xué)關(guān)系，推導(dǎo)包含以上參數(shù)的回歸方程和回歸曲線，得到一個計(jì)算海面溫度的經(jīng)驗(yàn)公式［4］。這種方法的優(yōu)點(diǎn)在于計(jì)算方便，能滿足一些實(shí)時性要求較高的系統(tǒng)需求，但是缺點(diǎn)也很明顯，因?yàn)榻?jīng)驗(yàn)公式是通過當(dāng)?shù)貙?shí)測數(shù)據(jù)的分析推導(dǎo)而來，存在一定的誤差，并且公式的局限性較大，如果要計(jì)算另一片海域的海面溫度，就需要重新分析實(shí)測數(shù)據(jù)和推導(dǎo)回歸方程。

計(jì)算海面溫度的另一種方法就是對海面溫度進(jìn)行數(shù)學(xué)建模，即利用熱平衡理論對海面建立熱平衡方程，把海面溫度當(dāng)作熱平衡方程的未知量，通過求解平衡方程中其它分項(xiàng)的值，最后再計(jì)算得到海面溫度。

本文將采用熱平衡方程來計(jì)算海面溫度，并針對海面熱平衡方程中的海面熱傳導(dǎo)分量的計(jì)算方法進(jìn)行研究。在求解方程過程中，首先通過隱式差分方法對方程進(jìn)行空間和時間維度上離散化，接著將方程轉(zhuǎn)換成矩陣運(yùn)算形式，最后利用追趕法求解海面溫度。

2 海面溫度場景模型建立

根據(jù)能量守恒定律，對于整個海水表面薄層來說，任意時刻，表面從外界吸收的能量與表面向外界散失的能量滿足以下熱平衡方程［5］：

其中各項(xiàng)含義如下：

R為凈入射輻射通量，即表面從外界（太陽、天空等）獲取的凈輻射通量，方向?yàn)閺纳舷蛳轮赶虮砻鏋檎?。不考慮外部熱源的情況下，其值等于太陽入射通量和天空入射通量與表面自發(fā)輻射散失的輻射通量之差：

G為海面向下熱傳導(dǎo)熱量，即地面表面朝下向更深層熱傳導(dǎo)散失的熱量，方向?yàn)閺谋砻嫦蛳聻檎?/p>

H為海面與大氣交換熱量，即表面與大氣發(fā)生熱交換散失的能量，方向?yàn)閺谋砻嫦蛏蠟檎?/p>

LE為海面蒸發(fā)熱量，即表面由于水分蒸發(fā)作用而散失的熱量，方向?yàn)閺谋砻嫦蛏蠟檎?/p>

除了G分量，方程剩余分量與傳統(tǒng)方法中的熱平衡方程一致。傳統(tǒng)方法計(jì)算得到的G分量表達(dá)式較為簡單，不能正確反應(yīng)出海水表面向下傳導(dǎo)的熱量，得到的結(jié)果誤差較大。

本文采用一維傅里葉熱傳導(dǎo)方程來計(jì)算海面向下傳導(dǎo)的熱量，求解方程過程中，對方程進(jìn)行空間和時間維度上進(jìn)行離散化，通過設(shè)定初始值迭代計(jì)算得到逐小時海面溫度。得到海面溫度結(jié)果較傳統(tǒng)方法更準(zhǔn)確，因此計(jì)算復(fù)雜度也較大，可以通過離線計(jì)算和查表方式來解決實(shí)時性問題。

接下來具體介紹熱平衡方程中每個分項(xiàng)的具體計(jì)算方式。

1）凈入射輻射熱量

凈入射輻射通量包括太陽入射熱量、天空入射熱量和海面自發(fā)輻射熱量，其中太陽輻射［6］為短波輻射，到達(dá)地球大氣層外的太陽輻射可以通過日地距離修正太陽常數(shù)精確計(jì)算得到；天空輻射［7］分為大氣輻射和云層輻射，均為長波輻射；物體通過輻射向外散發(fā)能量，由于地面物體的溫度不高，自發(fā)輻射也主要是長波輻射［7］。以上三種輻射能量的計(jì)算比較常用，這里就不詳細(xì)展開。

2）海面與大氣交換熱量

海面與大氣之間交換的熱量也稱感熱通量，熱量交換的速度和方向與表面與近地層空氣的溫差以及風(fēng)速等因素相關(guān)。通常采用空氣動力學(xué)公式進(jìn)行計(jì)算［8］：

其中，ρa(bǔ)為濕空氣的密度，Cp為空氣的定壓比熱容，Ts為表面溫度，Ta為大氣溫度，ra空氣動力學(xué)阻力［9］。

3）海面蒸發(fā)熱量

海面蒸發(fā)熱量又稱為潛熱通量，與感熱通量不同，該通量為表面水分蒸發(fā)作用帶走的熱量，對于海面水分可以采用以下公式［8］：

其中，L為水的汽化潛熱，ρv為絕對濕度。

4）海面向下熱傳導(dǎo)熱量

海面向下熱傳導(dǎo)熱量的計(jì)算滿足傅里葉方程，對于海面薄層，按一維熱傳導(dǎo)計(jì)算，方程為［10］

其中，T為表面層溫度，λ為材質(zhì)的導(dǎo)熱系數(shù)。特別地，表面層傳導(dǎo)的熱量G為

2.1 海面溫度場模型分析

由于在研究海面溫度計(jì)算中，我們忽略洋流、內(nèi)部漩渦和表面溫度差異等因素，認(rèn)為海面在某一時刻具有相同的溫度，并且在一定深度海水溫度是保持一個定值［11］，那么可以把研究對象縮小成具有一定深度的海水薄層，海洋表面滿足上方給出的表面熱平衡方程，海水薄層底部再往深處的海水則保持一個恒定的溫度，薄層內(nèi)部則滿足一維傅里葉熱傳導(dǎo)方程，只存在層與層之間的熱量傳遞。

經(jīng)過上述對海面模型的簡化和分析，將海水薄層沿著垂直方向上進(jìn)行空間劃分，從海面到薄層底部分成若干層，Zi為每一層的厚度，如果不考慮層與層之間溫度變化速度的差異，可以把每層厚度設(shè)成相同的值，本文采用等比數(shù)列設(shè)置薄層厚度，越往深處薄層厚度越大。

除了進(jìn)行空間劃分，還需要對海水薄層進(jìn)行時間上的劃分，由于紅外仿真系統(tǒng)是模擬海面一天24小時的輻射變化，所以時間維度上從0到24小時劃分，步長的△t設(shè)成定值。

2.2 海水薄層內(nèi)部溫度場

對于薄層內(nèi)部，由于我們忽略了洋流等復(fù)雜因素的影響，海水只存在層與層之間的能量傳遞，即滿足一維傅里葉熱傳導(dǎo)方程，那么將熱傳導(dǎo)方程從空間和時間維度上進(jìn)行離散劃分，轉(zhuǎn)為隱式差分方程［12］：

其中λi為每一層的導(dǎo)熱系數(shù)，由于忽略復(fù)雜因素的影響，每層海水的導(dǎo)熱系數(shù)相同；△Zi表示第i層的厚度，Tit表示在t時刻下的第i層海水的溫度；ρi表示第i層海水的密度；ci表示第i層海水的比熱容；這里每層的密度和比熱容都是相同值。

在公式中t時刻下每層溫度值和每一層厚度、海水導(dǎo)熱系數(shù)、密度、比熱容都是已知的，需要推導(dǎo)t+Δt下每一層的溫度，那么將上述差分公式變換成如下形式：

其中經(jīng)過公式換算，得到A、B、C、D系數(shù)表達(dá)式如下所示，其中i大于1，且小于n：經(jīng)過上述公式轉(zhuǎn)換，可以很方便地將方程式轉(zhuǎn)換成矩陣運(yùn)算，如下所示：

由于這里討論的是海水薄層內(nèi)部的能量交換，所以上述矩陣運(yùn)算中，從第二層到第n-1層的系數(shù)是已知的，第一層和最后一層的A、B、C、D系數(shù)仍是未知量，還需要對上下邊界進(jìn)行單獨(dú)推導(dǎo)分析。

2.3 海水薄層上邊界溫度場

海水表面不但要和第二層海水進(jìn)行熱量交換，還和大氣直接接觸，涉及到大氣背景輻射、風(fēng)速和相對濕度等因素的影響，所以一維傅里葉熱傳導(dǎo)方程不適合用于計(jì)算表面溫度。這里需要結(jié)合之前提到的表面熱平衡方程來解算表面溫度場數(shù)據(jù)，同樣采用隱式差分的方式，對該方程進(jìn)行空間和時間維度上離散化，得到如下公式：

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

和上小節(jié)類似，轉(zhuǎn)換成如下形式：

其中經(jīng)過公式換算，得到A、B、C、D系數(shù)如下所示：

通過對表面熱平衡方程進(jìn)行隱式差分轉(zhuǎn)換，得到矩陣運(yùn)算中系數(shù)A、B、C、D的表達(dá)式。值得一提是這里系數(shù)A值為0，這是由于海面不存在上一層海水，無需與上層海水交換熱量，所以不需要上層海水溫度的分量。

2.4 海水薄層下邊界溫度場

對于海水薄層下邊界而言，除了要和倒數(shù)第二層海水進(jìn)行熱量交換，它還要和薄層底部以下的海水進(jìn)行熱量交換；這個區(qū)別在于，我們把薄層底部以下海水經(jīng)過復(fù)雜的熱量交換，達(dá)到穩(wěn)定的溫度值，而這部分輻射熱量計(jì)算較為復(fù)雜，需要單獨(dú)處理。這里采用方式是將底部能量交換總結(jié)為一個換熱系數(shù)αin，從而底部熱傳導(dǎo)方程如下所示［12］：

其中，αin為換熱系數(shù)；Tin為海水薄層底部以下海水的溫度值，采用經(jīng)驗(yàn)數(shù)據(jù)；Tn為海水薄層底部的溫度值。

同樣由于該熱傳導(dǎo)方程存在微分項(xiàng)，需要對其進(jìn)行空間和時間維度上離散化，經(jīng)過隱式差分轉(zhuǎn)換得到如下公式：

和上小節(jié)類似，轉(zhuǎn)換成如下形式：

其中經(jīng)過公式換算，得到A、B、C、D系數(shù)如下所示：

2.5 海水薄層溫度場計(jì)算過程

通過對海水薄層上下邊界和內(nèi)部熱量傳遞方程的解算，最終將熱傳導(dǎo)方程轉(zhuǎn)成矩陣運(yùn)算，公式中系數(shù)矩陣和等號右邊的列向量為已知量，公式中第二項(xiàng)的列向量表示為下一Δt時刻每層海水的溫度值，為待求解。通過分析公式中第一項(xiàng)系數(shù)矩陣為三對角矩陣，其非零元素集中分布在主對角線和其相鄰兩個次對角線上，我們把這種矩陣公式稱為三對角方程組，可以通過追趕法來進(jìn)行求解［13］。

在利用追趕法求解該三對角方程組過程中，由于方程式從空間和時間維度上進(jìn)行離散化，意味著溫度場數(shù)據(jù)的求解是一個遞推的過程，需要擬定一個初始值才能解算出下一Δt時刻的海水薄層每層的溫度值。本文擬定以下起始條件：時間起點(diǎn)為凌晨4點(diǎn)，此時海面層溫度與大氣溫度相同，海水表層以下每層溫度按線性插值計(jì)算得到。設(shè)定該起始條件基于以下假設(shè)：經(jīng)過整個夜晚的降溫和熱交換，海水表面的溫度應(yīng)與大氣溫度相近。

求解過程如圖1所示，圖1左側(cè)為初始溫度場數(shù)據(jù)的設(shè)置，選定的時間為凌晨四點(diǎn)，海水表面溫度與大氣溫度相等，然后通過線性插值得到表層以下每層海水的溫度值，接著根據(jù)初始值按步長Δt進(jìn)行迭代，直到解算出仿真時刻的海水表面溫度。

3 實(shí)驗(yàn)結(jié)果分析

首先本文先進(jìn)行海面溫度與大氣溫度比較分析。由于海水表面直接接觸著大氣，大氣溫度直接影響著海水溫度，使得大氣溫度和海水溫度有著相同的上升下降趨勢；接著將海面溫度仿真數(shù)據(jù)和實(shí)際測量數(shù)據(jù)進(jìn)行對比驗(yàn)證，但是由于海面溫度逐小時數(shù)據(jù)難以獲取，比如中國氣象網(wǎng)中需要較高權(quán)限才能下載全球海平面溫度資料數(shù)據(jù)，普通學(xué)生用戶無法得到，只能從國家海洋環(huán)境預(yù)報(bào)中心網(wǎng)站中獲取鄰近海域的最低最高溫度數(shù)據(jù)，所以本文將用仿真計(jì)算得到的海面溫度數(shù)據(jù)與同一海域的最高最低溫度數(shù)據(jù)進(jìn)行比較，觀察仿真數(shù)據(jù)是否分布在最高最低溫度的范圍附近。

圖1 迭代計(jì)算海面溫度

本文采用網(wǎng)上氣象數(shù)據(jù)和海洋溫度數(shù)據(jù)的時間為2016年8月和2016年1月，分別代表著夏天和冬天季節(jié)，海域經(jīng)緯度位置約為東經(jīng)118.1°，北緯38.7°。其中海面溫度計(jì)算過程中所使用到的參數(shù)信息如表1所示。

表1 參數(shù)信息

圖2 夏天海溫和氣溫比較

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

圖3 冬天海溫和氣溫比較

圖4 夏天溫度數(shù)據(jù)對比

圖5 冬天溫度數(shù)據(jù)對比

從上面兩幅圖可以看出，雖然海面溫度仿真數(shù)據(jù)與實(shí)測海面溫度范圍有一定的誤差，這是因?yàn)楸疚暮雎粤搜罅鞯葟?fù)雜因素的影響，在整體上本文的仿真數(shù)據(jù)還是能很好地分布在最高最低值附近，較好反映實(shí)測海面溫度變化的范圍，從而表明本文海面溫度計(jì)算模型具有較好的正確性。

4 結(jié)語

本文從海面的熱平衡方程出發(fā)，介紹了平衡方程中的各項(xiàng)：太陽輻射、天空的長波輻射、海水的蒸發(fā)潛熱、海水對大氣的熱傳導(dǎo)和海水向深處的熱傳導(dǎo)分量，在建立海面熱平衡方程基礎(chǔ)上，利用隱式差分方法將海面模型在空間和時間維度上進(jìn)行離散化，接著將方程轉(zhuǎn)換成矩陣運(yùn)算形式，最后利用追趕法迭代求出仿真時刻的海面溫度數(shù)據(jù)。并對海面溫度的仿真數(shù)據(jù)進(jìn)行對比驗(yàn)證，分別和歷史大氣溫度數(shù)據(jù)、海面最高最低溫度數(shù)據(jù)進(jìn)行比較，實(shí)驗(yàn)結(jié)果表明仿真數(shù)據(jù)滿足實(shí)際溫度變化趨勢，說明本文海面溫度計(jì)算模型具有較好實(shí)用性，能為海面背景紅外輻射模型提供相關(guān)輸入，同時可以為海洋環(huán)境的紅外檢測提供參考。