基于風帶的海表面亮溫增量統計模型*

李長軍， 王穎芝， 趙 紅

(中國海洋大學數學科學學院，山東 青島 266100)

在世界各大洋中，海表面鹽度(Sea Surface Salinity, SSS)的變化是了解海洋在海洋生態系統和地球氣候系統中扮演何種角色的一個關鍵參數。迄今為止，Argo剖面浮標實現了SSS領域最好的時空覆蓋，提供分辨率約為3°×3°、每10天的SSS測量，然而，這個分辨率對于解決在SSS顯著變化的重要時空尺度上仍然太粗糙[1]。

SMOS衛星于2009年11月2日發射，是世界上最早能夠同時對土壤濕度和海水鹽度變化進行觀測的衛星，運轉期間，將提供分辨率為200×200 km2的10～30 d平均的全球鹽度數據，精度可達到0.1，但需要亮溫的測量精度控制在0.02～0.07 K范圍內[2]。在反演過程中，0.5K的亮溫誤差會導致1psu的鹽度反演誤差，且亮溫誤差和鹽度反演誤差成正相關，所以鹽度反演精度對亮溫精度有極為嚴格的要求[3]。L-波段海表面亮溫的測量值依賴于SSS的同時，也依賴于海表面粗糙度、海表面溫度、泡沫、入射角等其他因素，海表面粗糙度對亮溫的測量影響尤為明顯[4]。事實上，現實中的海洋是不平靜的，海洋表面因為風的存在變成由波浪構成的粗糙表面，海表面亮溫由兩部分組成，一部分是平靜海面輻射的亮溫，另一部分是由海表面粗糙度導致的亮溫增量[5-6]。在L-波段亮溫反演SSS的過程中，海表誤差的主要原因是風和泡沫誘導的海表面粗糙度對亮溫的影響。

為了提高SSS的反演精度，特別是利用正演模型來減少海面粗糙度引起的誤差，許多學者已經進行了大量的工作。在諸多亮溫增量模型中，歐洲空間局選擇3種正演模型對由海表面粗糙度導致的亮溫增量進行計算：模型1是將海表面假想成很多小波浪疊加在大波浪之上，海表面粗糙度按照截止波長λ0分為大尺度和小尺度的雙尺度模型[7]；模型2是利用微擾法(SPA)和小斜率近似(SSA)展開發射率理論的直接發射率模型[8]；模型3為半經驗模型，SMOS衛星采取查表法(Look Up
Table)來描述作為若干個參數函數的亮溫增量，該模型的輸入參數分別是：入射角θ，風速WSn，風向φ_WS和有效波高SWH[9]。這三種模型實現SMOS衛星鹽度反演得到的三種L2鹽度產品數據分別記作SSS1，SSS2，SSS3，并且這3種鹽度產品數據比較接近[10]。新版本V622的輸出產品為SSS1，將SSS1作為SMOS L2鹽度產品進行下面的分析與比較。通過本文構建的新亮溫增量模型獲得的海表鹽度產品記作SSS4。

一般來說，海水鹽度的變化較小，但會因海域所處的緯度位置不同有較大差異，主要受緯度、河流、洋流等因素的影響。SMOS衛星利用同一個亮溫增量模型在低緯度海域和開闊海域反演的鹽度產品精度較高，但在高緯度海域和近岸海域反演的鹽度產品精度非常低[11]。考慮到SSS的時空分布與全球大洋環流密切相關，而洋流的形成與全球風帶有著密切的聯系。因此，本文根據緯度劃分的風帶，依據風對海面粗糙度影響，將全球海域按照風帶劃分成東南信風帶(0°～30°S)、盛行西風帶(30°S～60°S)、極地東風帶(60°S～90°S)、東北信風帶(0°～30°N)、盛行西風帶(30°N～60°N)、極地東風帶(60°N～90°N)6個海域，分別建立亮溫增量模型，再計算各個海域的SSS。基于Argo實測鹽度數據、SMOS衛星數據，我們可計算出受海表面粗糙度影響的亮溫增量數據。在本文中，白冠覆蓋率(Whitecap coverage)、有效波高(Significant wave height,以下簡稱SWH)、海表面溫度(Sea surface temperature,以下簡稱SST)、降雨率(Rainfall rate)和蒸發量(Evaporation capacity)這5個變量的一階、二階、交叉項作為輸入變量，然后利用LASSO統計方法篩選出與該海域亮溫增量相關性較大的影響變量，在每個海域構建一個全新的二次曲線回歸亮溫增量模型。最后，將新模型嵌入SMOS衛星反演算法的程序計算中(L2OS軟件)，實現對軟件的修正。該研究提高了粗糙海面條件下亮溫數據的準確性，進一步提高了SSS的反演精度。

1 數據和方法

本次研究使用3個數據集，分別是Argo浮標數據集、SMOS數據集和輔助數據集。通過文獻[12-14]了解到，海表面粗糙度受海表面溫度、有效波高、白冠覆蓋率、降雨率和蒸發量的影響。新模型將選取這5個參量作為海表面粗糙度的影響變量，相關數據由歐洲中期天氣預報中心(ECMWF)提供。研究將集中在2015年5月的全球海洋，并在2015年6月進行驗證(由于數據的限制，以一個月為例，其他時間的數據結果將在后續繼續完成) 。

1.1 Argo實測數據

Argo(Array for Real-time Geostrophic Oceanography)數據由法國CORIOLIS數據中心(http://www.coriolis.eu.org)提供，作為本文的實測數據。Argo項目是一個實時的海洋觀測項目，Argo浮標可以測量海水溫度、鹽度等海洋參數，具有自定位能力。取Argo測量的有效的第一層(-5 m)鹽度數據作為海表面鹽度的實際值，選擇從2015年5月1日～2015年5月31日，在非極地水域均勻分布的13 783組數據進行建模，每個Argo數據組包括日、小時、緯度、經度和鹽度。圖1顯示了每個可用的Argo點的位置。

圖1 2015年5月份Argo測量點的位置

1.2 SMOS數據

SMOS衛星數據來自歐洲空間局(ESA)，包括L1C亮溫產品(數據版本為V620，https://smos-ds-02.eo.esa.int/oads/access/)、L2海水鹽度產品(數據版本為V622)。SMOS衛星搭載的MIRAS儀器(Microwave Imaging Radiometer with Aperture Synthesis)能夠探測地球表面L-波段(1.4GHz)的微波輻射，重建一個二維的、擁有多個入射角的亮溫視野圖像。L2海水鹽度產品由L2OS軟件生成，在SSS反演過程中，系統每2周進行一次OTT(Ocean Target Transformation)修正，旨在修正MIRAS實測亮溫與模型模擬亮溫之間總是存在的系統偏差，但季節偏差和緯度偏差依然存在[1]。

1.3 輔助數據

歐洲中期天氣預報中心(The European Centre for Medium-Range Weather Forecasts，以下簡稱ECMWF，https://www.ecmwf.int/)集合了全球大部分的地面氣象觀測數據以及90%以上的環境衛星數據，可提供多種類型的數據用于地球地表及其周圍大氣環境的研究[15]。本文使用的ECMWF數據包括海表面溫度、有效波高、白冠覆蓋率、降雨率和蒸發量，數據版本為V318。

SMOS亮溫數據與ECMWF輔助數據使用同一種等區域網格系統—ISEA 4-9。在1.2和1.3節中，利用2015年5月8日～2015年5月31日的所有上行和下行數據，建立亮溫增量模型，并2015年6月1日～2015年6月30日的所有上行和下行數據驗證模型。

1.4 數據匹配

由于Argo數據、SMOS衛星數據和輔助數據存在時空上的差異，因此在構建模型之前，需對南北半球的東南(北)信風帶、盛行西風帶、極地東風帶6個海域的訓練數據進行時間、空間上的匹配，包含以下幾個步驟[16]：首先，將2015年5月的Argo數據按照經緯度信息與ISEA 4-9全球離散格網匹配，匹配遵循距離最近原則，最終獲得12 111組全球有效數據。然后，將2015年5月的SMOS衛星數據、ECMWF數據按照網格點編號、最大時間間隔為1天與Argo數據進行時空匹配。考慮到不同入射角的問題，最后，得到南半球數據各為63 694、100 897和14 859組(H極化)，63 621、100 834和14 853組(V極化)；北半球數據各為80 693、76 887和5 916組(H極化)，80 732、76 964和5 899組(V極化)。這些數據和由1.5節計算得到的亮溫增量數據為本次研究的訓練數據集，是構建6個海域亮溫增量模型的數據基礎。

1.5 亮溫增量

海表面亮溫可以分解為兩部分：平靜海面輻射的亮溫和由海表面粗糙導致的亮溫增量。海表面亮溫在p極化(H水平極化，V垂直極化)時可以表示成以下形式[9]:

Tb,p(θ,SST,SSS,Prough)=Tb,p,flat(θ,SST,SSS)+ΔTb,p,rough(θ,SST,SSS,Prough) 。

(1)

式中：Tb,p,flat是平靜海面輻射的亮溫；ΔTb,p,rough是由海表面粗糙度導致的亮溫增量；Prough為粗糙度參數；θ為入射角。

平靜海面的亮溫被定義為：

Tb,flat(θ)=e(θ)·SST。

(2)

e(θ)是L-波段的海表面發射率，攜帶了關于鹽度的主要信息，與反射率Γ存在e=1-Γ的等式關系。在平靜海面中，反射率可以通過菲涅爾反射定律直接計算，并提供一個精確的介電常數模型，對于每個極化方式下的反射系數R，均定義為海水介電常數ε和入射角θ的函數：

(3)

(4)

因此，平靜海面的Tb,h,flat和Tb,v,flat可表示為：

Tb,h,flat(θ,SST,SSS)=(1-Rh(θ,ε))·SST,

(5)

Tb,v,flat(θ,SST,SSS)=(1-Rv(θ,ε))·SST。

(6)

海水的復介電常數ε依賴于溫度和鹽度，在微波段的任何頻率都可依據Debye表達式計算[17]：

(7)

式中：i是虛數單位；ε∞是無線高頻介電常數；εs是靜態介電常數；ε0= 8.854×10-12[F/m]為真空下的介電常數；τ是豫馳時間；σ是離子電導率，ω=8.7965×109為電磁波的角頻率。

基于1.4節匹配后的數據，可以計算出平靜海面下，與Argo SSS和其他參數數據相關，各個入射角θ下的Tb,h,flat和Tb,v,flat，亮溫增量可以表示為：

ΔTb,p,rough(θ,SST,SSS,Prough)=Tb,p(θ,SST,SSS,Prough)-Tb,p,flat(θ,SST,SSS) 。

(8)

Tb,p是由SMOS衛星通過反演迭代算法提供的亮溫值，在這里認為是準確的。儀器校正后仍然存在的亮溫偏差通過OTT進行校正，但是目前業務上OTT算法仍然受正演模型缺陷的影響，表現在當天頂角和入射角不同時，誤差分布將有所不同，隨著地表條件參數的改變偏差訂正會有0.5K左右的變化[18]。

到目前為止，已經得到了所有的變量數據作為訓練數據集，亮溫增量是因變量(輸出變量)，而海表面溫度、有效波高、白冠覆蓋率、降雨率和蒸發量是自變量(輸入變量)。然后，采用LASSO統計方法和二次曲線回歸模型建立了亮溫增量的統計模型。

1.6 LASSO統計方法

Tibshirani[19]在1996年提出了LASSO(The Least Absolute Shrinkage and Selection Operator)算法,通過該方法可以做到變量選擇，即去除冗余變量，得到更重要的解釋變量。它的具體表達式為：

(9)

通過選取非零系數的變量，得到一個可解釋的模型。然后將所選的變量用于多元線性回歸(MLR)中作為輸入項，計算由海面粗糙度引起的亮溫增量。所有匹配的數據點都將在訓練模型中使用。

1.7 二次曲線回歸模型

以二維自變量為例，假設因變量Y和自變量X1，X2之間有如下關系，稱為二次曲線回歸模型[20]：

(10)

其中，E(ε)=0且Var(ε)=σ2。

(11)

這里“經驗”兩字表示這個回歸方程是基于前面的n次觀測數據(x1i，x2i，yi)，i=1,…,n獲得的。

2 結果

在SMOS衛星的亮溫觀測中，同一個測量點會包括多個入射角θ對應的亮溫結果，但其他輔助數據保持不變。本次研究僅對入射角為θ0的ΔTb,v,rough和ΔTb,h,rough建立統計模型，其他角度的亮溫增量通過以下公式進行轉換[21]:

(12)

式中：θ0為42.5°入射角；ρp(θ0)是不同極化方式下入射角為θ0的菲涅爾反射率。

2.1 模型建立

利用1.4和1.5節獲得的訓練數據集，選擇白冠覆蓋率、有效波高、海表面溫度、降雨率和蒸發量的一階、二階和交叉項，共20個變量，作為待選變量。然后利用LASSO統計方法進行變量選擇得到解釋亮溫增量的重要變量，再建立二次曲線亮溫增量模型。最終篩選出各個海域在42.5°入射角，不同極化方式下5個相同的變量，分別是：SWH和SST的一階、二階和交叉項，表明這些變量在全球海洋中的重要作用。最后，使用選出的變量構造訓練模型，并代入表1所示的估計回歸系數。

ΔTb,p,rough(θ,SST,SSS,Prough)=A+B·XSWH+C·XSST+D·XSWH2+E·XSST2+F·XSWH·SST。

(13)

式中：XSWH2和XSST2是 SWH和SST的二階項，XSWH·SST是SWH和SST的交叉項。

表 1 全球6個海域的亮溫增量模型擬合系數

Note：①Area;②Polarization;③Southeast trade-wind zone;④Prevailing westerlies;⑤Polar easterlies;⑥Northeast trade-wind zone;⑦Prevailing westerlies;⑧Polar easterlies.

雖然在不同的海域通過LASSO選出了相同的影響變量，但這些變量的系數在不同的區域是不同的，表明每個變量對亮溫增量的影響是區域相關的。因此，在全球海洋中使用同一模型來反演SSS是不合適的。

模型建成后，將新模型替換L2OS軟件中現有的亮溫增量模型，并將反演得到的SSS產品命名為“SMOS SSS4”。接下來，將在訓練數據集和驗證數據集上測試結果。

2.2 訓練數據集誤差分析

依據2015年5月的Argo浮標所在網格點，輸出各個海域新模型反演鹽度數據SSS4，原模型反演鹽度數據SSS1與該時期、該區域的Argo浮標SSS數據進行匹配，匹配結果：東南信風帶(0°～30°S)1242組，盛行西風帶(30°S～60°S)2028組，極地東風帶(60°S～90°S)284組，東北信風帶(0°～30°N)1466組，盛行西風帶(30°N～60°N)1280組，極地東風帶(60°N～90°N)285組。經過計算，新模型在全球海域的絕對平均誤差(MAE)和均方根誤差(RMSE)分別是0.956 3和1.519 6，而SMOS SSS1產品的分別為1.711 7和2.528 8。各個海域的誤差如表2所示，在全球海域和各個海域，新模型的鹽度反演精度均優于原模型，同時繪制了SSS4和SSS1的誤差百分比分布圖(見圖2，其中，model error對應SSS4誤差，SMOS error對應SSS1誤差)。

表2 訓練數據集上的SSS4、SSS1與Argo實測SSS的誤差對比

Note：①Area;②Error;③Southeast trade-wind zone;④Prevailing westerlies;⑤Polar easterlies;⑥Northeast trade-wind zone;⑦Prevailing westerlies;⑧Polar easterlies.

2.3 驗證數據集誤差分析

為了對模型進行驗證，依據2015年6月Argo浮標所在網格點，全球共11 297個，輸出新舊模型反演的鹽度，與Argo SSS匹配，匹配結果：東南信風帶(0°～30°S)1 741組，盛行西風帶(30°S～60°S)2 438組，極地東風帶(60°S～90°S)139組，東北信風帶(0°～30°N)2 351組，盛行西風帶(30°N～60°N)1 461組，極地東風帶(60°N～90°N)204組。圖3是SSS4和SSS1的誤差分布圖。圖中30°N～60°N以及60°～90°(NS)中有偏離Argo數值較大的點，通過查看數據發現，30°N～60°N范圍內陸地較多，偏離Argo數值較大的點大多分布在海陸交界處，受RFI影響較大，故反演誤差大。而60°～90°(NS)為高緯度區域，Argo浮標少，用于建模的數據點少，導致模型自身的誤差偏大，故反演誤差大。經過計算，新模型和原模型在全球海域的MAE、RMSE分別為0.971 1和1.591 6，1.691 9和2.482 2，按風帶劃分的各個海域的鹽度反演誤差如表3所示，進一步說明新模型生成的鹽度產品比SMOS L2鹽度產品具有更高的精度，從實驗的角度肯定了按風帶劃分全球海域反演鹽度的正確性。

圖2 訓練數據集上的SSS4和 SSS1誤差百分比直方圖

圖3 驗證數據集上的SSS1和 SSS4 誤差分布圖

Note：①Area;②Error;③Southeast trade-wind zone;④Prevailing westerlies;⑤Polar easterlies;⑥Northeast trade-wind zone;⑦Prevailing westerlies;⑧Polar easterlies.

取2015年6月1～10日的全球SSS數據進行平均，利用Surfer11軟件繪制新模型的十日平均SSS圖(見圖4右)，從ESA網站上下載SMOS衛星L3十日平均海水鹽度產品數據，得到SMOS衛星L3十日平均SSS圖(見圖4左)，同時繪制了SMOS L3 SSS與新模型SSS4的十日平均差值圖(見圖5)。可以看出2幅圖的SSS變化趨勢基本一致，圖中新模型鹽度值異于L3產品的海域大多靠近陸地，以及分布在高緯度海域，受RFI、數據點少模型誤差大等因素的影響。地中海區域蒸發旺盛，周圍河流淡水注入少，海水鹽度高，新模型的反演精度高于原模型，所以新模型鹽度值高于原模型是合理的。圖中顯示新模型的十日平均SSS圖比SMOS衛星L3十日平均SSS圖的鹽度覆蓋范圍更為廣闊，更加完整地描繪了全球鹽度的分布情況，具有更高的實際應用價值。

圖4 SMOS L3 SSS與新模型SSS4的十日 (2015年6月1～10日)平均圖

圖5 SMOS L3 SSS與新模型SSS4的十日平均差值圖

3 結語

本文首次提出將全球海域按照風帶劃分成6個部分，利用LASSO統計方法挑選出包括SWH和SST一次項、二次項、交叉項在內的5個影響變量，在各個海域重新建立了有關亮溫增量的二次曲線回歸模型。在驗證數據集上將新模型反演的SSS4產品與Argo實測SSS進行比較：南半球的東南信風帶(0°～30°S)、盛行西風帶(30°S～60°S)、極地東風帶(60°S～90°S)以及北半球的東北信風帶(0°～30°N)、盛行西風帶(30°N～60°N)、極地東風帶(60°N～90°N)6個海域新模型對亮溫修正后反演的SSS4，MAE分別為0.76、 0.88、 0.93、 0.92、1.28和1.21。相對應的SMOS衛星L2日平均鹽度產品數據(SSS1)的MAE分別為0.98、 1.61、 2.82、 1.50、 2.35和3.13。結果表明本次研究得到的亮溫增量模型更適用于全球海域SSS研究，可以輸出更高精度的SMOS衛星SSS產品。通過比較誤差的直方圖和分布圖發現，新模型的鹽度反演誤差更加集中，且SSS變化趨勢基本一致。

從整體上看，南半球海域的鹽度反演精度高于北半球海域，因為北半球陸地覆蓋多于南半球，受陸地射頻(RFI)干擾大，因此模型誤差較大。由于實驗數據的不易獲得以及完成一次實驗耗費的時間長達1個月這兩方面限制，模型的驗證數據集僅停留在6月，與5月的訓練數據集在時間上間隔較短，難以斷定時間對模型預測的干擾程度，可繼續利用7月、8月，甚至更大時間尺度的驗證數據集評估模型，以獲得模型預測的最大時間范圍。同時，在實驗中未考慮不同風帶邊界處是否存在不連續現象，也未對邊界數據進行單獨處理。在未來的研究中，可以通過考慮更多的影響因素，適當的結合相關的亮溫物理機制，建立更為復雜的數學模型，更加準確地刻畫亮溫增量，得到更高精度的SSS產品。