數(shù)據(jù)融合技術在海洋二號衛(wèi)星數(shù)據(jù)中的應用

齊亞琳 林明森

（1 中國空間技術研究院，北京 100094） （2 國家衛(wèi)星海洋應用中心，北京 100081）

1 引言

我國第一顆海洋動力環(huán)境探測衛(wèi)星——海洋二號（HY-2）于2011年8月成功發(fā)射，可執(zhí)行業(yè)務化監(jiān)測和全球海洋動力環(huán)境參數(shù)探測。在HY-2衛(wèi)星數(shù)據(jù)產(chǎn)品中，為最大限度地保存觀測數(shù)據(jù)的信息，并節(jié)省存儲空間，數(shù)據(jù)產(chǎn)品采用沿軌存儲的方式，即沿衛(wèi)星軌道方向存儲風場、海面溫度等海洋參數(shù)，并同時記錄對應的觀測時間等信息；并且HY-2衛(wèi)星搭載的微波散射計和掃描微波輻射計的刈幅寬度雖然大于1 400km，但一天沿軌數(shù)據(jù)仍不能完全覆蓋全球海域。而海洋環(huán)境預報等部門在實際使用衛(wèi)星遙感數(shù)據(jù)中，通常需要用到全球化網(wǎng)格化的標準數(shù)據(jù)，并要求數(shù)據(jù)具有較高的空間和時間分辨率。為了解決上述問題，數(shù)據(jù)融合被認為是非常有效的手段之一。通過數(shù)據(jù)融合可以優(yōu)化遙感信息資源利用，實現(xiàn)多源數(shù)據(jù)的優(yōu)勢互補。

本文以HY-2衛(wèi)星獲取的海面風場和海面溫度作為數(shù)據(jù)源，結(jié)合美國國家環(huán)境預報中心（National Centers for Environmental Prediction，NCEP）的再分析數(shù)據(jù)，分別進行了海面溫度和海面風場的數(shù)據(jù)融合，并得到了數(shù)據(jù)融合后的結(jié)果。通過定性分析，結(jié)果表明文中使用的融合方法對數(shù)據(jù)改善明顯，并且可用于探測中、大尺度的海洋現(xiàn)象，證明數(shù)據(jù)融合技術可為HY-2衛(wèi)星高分辨率數(shù)據(jù)產(chǎn)品的開發(fā)提供借鑒。

數(shù)據(jù)融合算法主要包括逐步訂正算法、混合分析算法、截斷濾波窗的變分分析算法、客觀分析算法、三維變分同化分析算法、卡爾曼濾波算法和最優(yōu)插值算法。逐步訂正算法由于僅使用當前的觀測數(shù)據(jù)，因此時間上的記憶性差；混合分析算法也沒有時間上的記憶性，且時空分辨率都很低，分析場被嚴重平滑；統(tǒng)計最優(yōu)估計算法和使用截斷濾波窗的變分分析算法均選用了包括紅外與微波在內(nèi)的多種衛(wèi)星遙感數(shù)據(jù)，空間分辨率高，但在有多種衛(wèi)星遙感數(shù)據(jù)時，僅選取其中一種數(shù)據(jù)則可能給融合結(jié)果帶入較大的噪聲；三維變分同化分析算法對計算機存儲要求低，計算速度快，三維變分同化分析和最優(yōu)插值算法分析一樣，融合結(jié)果受到數(shù)據(jù)分布稀疏的限制；最優(yōu)插值算法可看作是卡爾曼濾波法的簡化版，與卡爾曼濾波法相比，最優(yōu)插值算法有以下幾個方面的特點：在誤差的計算上，從卡爾曼濾波法的公式看，它的誤差是發(fā)展的，這是卡爾曼濾波法最吸引人的地方，但這一直也是制約著卡爾曼濾波法應用的致命之處，因為它的計算量很大，而最優(yōu)插值算法完全省略了這點。因此最優(yōu)插值算法減少了很大的計算量，使得其誤差變?yōu)殪o止的；在最優(yōu)插值算法的資料選擇方面，它通常只是選取分析點附近的資料來進行分析。這樣，大大減少了計算量；最優(yōu)插值算法通常只進行單變量分析（例如海面溫度）；時間空間距離加權(quán)插值的優(yōu)點是公式比較簡單，特別適用于結(jié)點散亂、不是網(wǎng)格點的問題。

總的來說，最優(yōu)插值算法由于其計算量少，融合的“性價比”高，因此可以在單變量數(shù)據(jù)融合——海面溫度融合業(yè)務上得到廣泛的應用［1］?；谝陨峡紤]，本文選用最優(yōu)插值算法對海面溫度數(shù)據(jù)進行融合，并且可利用該方法進行工程化和業(yè)務化。

根據(jù)散射計數(shù)據(jù)和插值方法本身的研究發(fā)現(xiàn)，時間空間權(quán)重插值法最適合結(jié)點散亂、非網(wǎng)格點的散射計數(shù)據(jù)融合，并且實現(xiàn)起來也是比較容易，適宜工程化和業(yè)務化操作。

2 海洋二號衛(wèi)星海面溫度數(shù)據(jù)融合

2．1 海洋二號衛(wèi)星海面溫度數(shù)據(jù)融合方法

2．1．1 最優(yōu)插值算法

在最優(yōu)插值算法中，空間網(wǎng)格點上的分析值是由網(wǎng)格點的背景場加上修訂值而確定的，其修訂值由周圍各觀測點的觀測值與背景場值的偏差加權(quán)求得，其權(quán)重系數(shù)（即最優(yōu)插值系數(shù)）不是任意選擇的，應該使得網(wǎng)格點分析值的誤差達到最小。對于海面溫度，最優(yōu)插值的表達式為［2］

式中：vam代表海面溫度在空間網(wǎng)格點上的分析值；vem代表海面溫度在空間網(wǎng)格點上的背景場值；K為權(quán)重系數(shù)矩陣；v0s代表海面溫度在觀測點的觀測值；ves代表海面溫度在觀測點的背景場值。當觀測點和網(wǎng)格點不重合時

式中：H為插值算子。

合并式（1）和式（2）得vam＝vem＋K（v0s-Hvem）。如果背景場和觀測值的空間分辨率不相同，則背景場和觀測值所屬的權(quán)重不同，也即H不同。所以本研究分別用H1，H2作為雙線性插值算子，分別對應背景場和觀測點，通過該算子分別將網(wǎng)格點的背景場值和觀測值插值到分析點上，轉(zhuǎn)換為式（3）。

經(jīng)過數(shù)據(jù)預處理后，可直接進行式（3）的運算，式（3）中H1、H2及K分別為雙線性插值算子和權(quán)重系數(shù)矩陣。

權(quán)重系數(shù)矩陣K的表達形式為

式中：B為背景場誤差協(xié)方差矩陣，R為觀測誤差協(xié)方差矩陣。欲計算權(quán)重系數(shù)K，必須先估算背景場誤差協(xié)方差矩陣B和觀測誤差協(xié)方差矩陣R。

2．1．2 背景場誤差及觀測誤差協(xié)方差矩陣計算

對B的確定基于以下兩個假設條件：（1）B是定常的；（2）背景場誤差的水平相關性滿足水平距離的增加相關性呈指數(shù)遞減的規(guī)律。則

式中：ρ為預報背景場水平相關系數(shù)矩陣，D為背景場誤差組成的對角線矩陣。對于觀測誤差協(xié)方差矩陣R采用與背景場誤差協(xié)方差矩陣B相同的處理方法。

對于相關系數(shù)矩陣ρ，鑒于以往的經(jīng)驗和觀測站點橫向比縱向稀疏分布的特點［3-5］，對相關系數(shù)矩陣做了相應的改進，使其更符合相關尺度的物理規(guī)律。相關系數(shù)矩陣ρ定義為

式（6）的相關性范圍呈橢圓形狀。式中：a為經(jīng)度方向的相關尺度；b為緯度方向的相關尺度；，分別為經(jīng)度、緯度方向上的距離。在本文中a取200km，b取150km。

2．1．3 矩陣的求解

由于觀測點的不規(guī)則分布，在求解權(quán)重系數(shù)矩陣K時，（HBHT＋R）-1的求解過程通常出現(xiàn)誤差，因此，應在實際求解時應盡可能避免求逆過程。本文采用求解方程的方法來求解權(quán)重系數(shù)矩陣K。將式（4）進行改寫，兩邊同時右乘（HBHT＋R），得

對式（7）兩邊再同時進行矩陣轉(zhuǎn)置，得

因此，將對K的求解轉(zhuǎn)化為對KT的求解，且避免了對矩陣的求逆。但是因為稀疏矩陣H數(shù)據(jù)元素分布的不規(guī)則性（由于觀測數(shù)據(jù)分布不規(guī)則導致），得到的轉(zhuǎn)置矩陣方程組可能是病態(tài)的。對于病態(tài)方程組的求解，通常有三種方法求解：第一種方法是采用病態(tài)迭代算法，即用全選主元高斯法求解，得到方程組的一組近似解，將近似解帶回原方程，得到剩余向量方程，然后求解剩余向量方程，直到滿足精度為止。第二種方法是采用QR 分解與Household轉(zhuǎn)換相結(jié)合，轉(zhuǎn)化為求解線性方程組的最小二乘問題，得到線性方程組的最小二乘解。第三種方法是采用奇異值分解（SVD），得到線性方程組的最小二乘解。本文采用奇異值分解方法來處理可能出現(xiàn)的病態(tài)方程。將K分解為3個矩陣

式中：U為下三角矩陣；S為對角線矩陣；VT為上三角矩陣。SVD 法可以有效地消除不正常解，并滿足方程求解精度。

為了解決計算機內(nèi)存不足和觀測數(shù)據(jù)數(shù)量稀少的問題，采用逐點逐次插值過程。其基本思路是以點為基本計算單位，尋求相關半徑內(nèi)的相關點，并計算其相關性，求得該點的融合值，同時在空間上逐行逐列推進。

2．2 海洋二號衛(wèi)星海面溫度數(shù)據(jù)融合結(jié)果

圖1所示為HY-2 衛(wèi)星掃描微波輻射計2011年12月3日升軌（白天）和降軌（夜間）的逐日全球海面溫度（SST）分布，圖2為使用最優(yōu)插值算法制作的HY-2衛(wèi)星掃描微波輻射計全覆蓋的2011年12月3日白天和夜間全球SST 分布。

從圖2可以發(fā)現(xiàn)，融合后的全球海面溫度，不再存在由于軌道覆蓋范圍不夠?qū)е碌能壍篱g空隙，而且全球海面溫度的空間分辨率不降低。

融合前，HY-2衛(wèi)星掃描微波輻射計測量的海面溫度不能完全反應全球海面溫度的分布情況，特別是不能反應出如黑潮這樣的西邊界洋流（圖1）。經(jīng)數(shù)據(jù)融合后，全球海面溫度分布非常明顯地分為熱帶、亞熱帶、寒帶等，并且能反應出諸如黑潮等高溫海流區(qū)域，如圖3所示。

圖1 HY-2衛(wèi)星掃描微波輻射計數(shù)據(jù)2011年12月3日的全球SST 分布Fig．1 SST of the world ocean from HY-2satellite radiometer data on December 3，2011

圖2 融合后HY-2衛(wèi)星掃描微波輻射計數(shù)據(jù)2011年12月3日全覆蓋的全球SST 分布Fig．2 SST of the world ocean from the fused HY-2radiometer data of December 3，2011

圖3 融合后的HY-2衛(wèi)星掃描微波輻射計數(shù)據(jù)西北太平洋區(qū)域SST 分布Fig．3 SST of the west north Pacific from the fused HY-2radiometer data

大洋中的中尺度渦旋有冷渦和暖渦之分，顧名思義冷渦在海面表現(xiàn)為海面溫度相對臨近海域低、暖渦在海表面表現(xiàn)為海面溫度較臨近海域高。融合后的海面溫度可用來探測全球海洋的中尺度渦旋，而融合前不具備這樣的能力。除此之外，融合后的海面溫度可用來研究厄爾尼諾與南方濤動（El Ni?o-Southern Oscillation，ENSO）等海洋災害，而融合前的數(shù)據(jù)由于不完整，不具有真正意義上的使用價值。

綜上所述，針對HY-2衛(wèi)星掃描微波輻射計獲得的海面溫度數(shù)據(jù)，本文使用最優(yōu)插值算法進行單變量融合，融合結(jié)果在不降低精度的基礎上提高了空間分辨率。完全可以在現(xiàn)有的基礎上，融合國內(nèi)外同類或紅外衛(wèi)星獲取的海面溫度數(shù)據(jù)，以提高空間分辨率和精度，并實現(xiàn)業(yè)務化和工程化。

3 海洋二號衛(wèi)星海面風場數(shù)據(jù)融合

3．1 海洋二號衛(wèi)星海面風場數(shù)據(jù)融合方法

采用時間空間權(quán)重插值法對多源衛(wèi)星遙感海面風場進行融合。該方法實質(zhì)上是一種比較簡單的時

式中：Uest表示最終獲得的估計值；下標k（k＝1，2，…，n）表示衛(wèi)星觀測數(shù)據(jù)點，Wk表示在k點的權(quán)重，Uk表示在k點的觀測值，（xk，yk，tk）表示k點的時空坐標；下標o表示網(wǎng)格點，（xo，yo，to）表示其時空坐標；n表示插值點個數(shù)；R為距離影響半徑，T為時間影響半徑。

本文選定的時間空間權(quán)重插值法為普通Kriging插值法進行插值。Kriging 插值法是法國科學家Matheron［6］提出的，是對空間分布的數(shù)據(jù)求線性最優(yōu)、無偏內(nèi)插估計的一種方法，并在獲得預測結(jié)果同時可以獲得預測誤差。其基本步驟如下：

1）求變異函數(shù)

變異函數(shù)r（h）反映了區(qū)域化變量的空間自相關性。

式中：h為觀測點之間的空間間隔距離；N（h）為距離等于h的點對數(shù)；Z（xi）為處于點xi處變量的觀測值；為與點xi偏離h處變量的實測值。在實際計算中，也可以將所有觀測點的相對距離劃分為若干等級，計算每個等級內(nèi)觀測點的個數(shù)，然后對每個等級內(nèi)的所有點數(shù)取距離的平均值及r（h）的平均值。將計算所得的（h，r（h））點連接后就可以得到實驗變異函數(shù)，再以最小二乘法計算出理論變異函數(shù)及其參數(shù)。理論變異函數(shù)模型有線性模型、指數(shù)模型、球面模型、高斯模型等。通過比較選取的高斯模型作為理論變異函數(shù)。

式中：C0稱為塊金常數(shù)；C0＋C為基臺值；，C為拱高。

2）利用變異函數(shù)求權(quán)重系數(shù)

假設為未觀測點，xi（i＝1，2，…，n）為其周圍的觀測點，Z表示計算的區(qū)域化變量。則對處某個區(qū)域化變量的估計值為

式中：λi表示權(quán)重，可由理論變異函數(shù)求得，為了保證是無偏估計，要求＝1。式（13）也可以改寫為矩陣形式

其中

式中：r（x1，xn）為點x1和xn之間的變異量，可由理論變異函數(shù)根據(jù)兩點間距離求出。同時該點估計值誤差的平方可由下式求得δ2（x0）＝BTW-1B。

3．2 海洋二號衛(wèi)星海面風場數(shù)據(jù)融合結(jié)果

圖4和圖5 分別是同時間內(nèi)的HY-2 衛(wèi)星微波散射計風場和美國國家環(huán)境預報中心（NCEP）再分析數(shù)據(jù)。該軌道的HY-2 衛(wèi)星微波數(shù)據(jù)不能覆蓋東海和南海，而NCEP 數(shù)據(jù)可覆蓋全球范圍內(nèi)的海域。故在實際應用中需將高分辨率的HY-2風場數(shù)據(jù)和不受任何軌道等影響的NCEP數(shù)據(jù)進行融合。圖6是圖4 和圖5 所示的HY-2 微波散射計風場數(shù)據(jù)和NCEP數(shù)據(jù)的融合結(jié)果。圖6 顯示的融合產(chǎn)品中，風場的空間分辨率達到與HY-2風場空間分辨率相同，并且覆蓋東海和南海等海域。HY-2衛(wèi)星微波散射計一天最多可覆蓋全球90%的海域，而融合后的產(chǎn)品可100%覆蓋全球海域。

圖4 用于融合的海洋二號散射計風場Fig．4 HY-2satellite scatterometer wind field to be fused

圖5 用于融合的NCEP風場Fig．5 NCEP wind field to be fused

圖6 西北太融合海面風場Fig．6 Fused wind field in the west north Pacific

從圖6可以看出，融合后的海面風場不僅彌補了高分辨率的HY-2衛(wèi)星遙感數(shù)據(jù)不能全覆蓋的不足，同時提高了NCEP 數(shù)據(jù)的空間分辨率，可為兩種數(shù)據(jù)更好地使用，尤其是為業(yè)務化海洋監(jiān)測應用提供便利。

另外，從圖6可以看出，融合后的海面風場在東海黑潮上空加速，風向平行于東海黑潮主軸，這與已有的研究結(jié)果［7-10］一致。雖然海面溫度對海面風速影響機制目前尚不完全清楚，但融合后的海面風場可為研究這一機制提供非常好的數(shù)據(jù)源。

目前國際上同類衛(wèi)星微波散射計數(shù)據(jù)非常稀少，但在不久的將來同類衛(wèi)星數(shù)據(jù)會逐漸增多。同時，目前國際上有較為成熟的風場再分析數(shù)據(jù)，如NCEP，歐洲中期天氣預報中心（European Centre for Medium-Range Weather Forecasts，ECMWF）的數(shù)據(jù)等，可以用來和HY-2衛(wèi)星散射計數(shù)據(jù)進行融合，并實現(xiàn)業(yè)務化和工程化。

4 結(jié)束語

本文針對HY-2衛(wèi)星獲取的海面溫度和海面風場數(shù)據(jù)，利用不同的融合方法對HY-2衛(wèi)星海面溫度和海面風場數(shù)據(jù)進行融合。通過定性分析，融合后的海面溫度和海面風場數(shù)據(jù)在空間分辨率不降低的前提下，彌補了自身的不足。

本文的研究結(jié)果證明數(shù)據(jù)融合方法應用到HY-2衛(wèi)星數(shù)據(jù)處理系統(tǒng)中，可以較好地解決衛(wèi)星觀測數(shù)據(jù)覆蓋范圍不足的問題。這為今后衛(wèi)星數(shù)據(jù)的高效使用和海洋環(huán)境預報部門更有效地使用衛(wèi)星數(shù)據(jù)，提供了一種技術途徑。

從本文的研究結(jié)果可以發(fā)現(xiàn)，HY-2 衛(wèi)星數(shù)據(jù)經(jīng)融合后得到的海面溫度、海面風場數(shù)據(jù)，可為全球海氣相互作用研究提供同時間內(nèi)相同區(qū)域的數(shù)據(jù)源，對進一步找到海面溫度和海面風場的定量物理關系奠定基礎。同時，HY-2衛(wèi)星海面溫度和海面風場數(shù)據(jù)融合的成功實現(xiàn)和廣泛應用，可為HY-2衛(wèi)星其它數(shù)據(jù)的高效使用提供技術基礎，并實現(xiàn)HY-2衛(wèi)星與國內(nèi)外同類衛(wèi)星數(shù)據(jù)融合的業(yè)務化和工程化。

（References）

［1］奚萌．基于最優(yōu)插值算法的紅外和微波遙感海面溫表溫度數(shù)據(jù)融合［D］．北京：國家海洋環(huán)境預報中心，2011

Xi Meng．Merging infrared radiometer and microwave radiometer sea surface temperature data based on the optimum interpolation［D］．Beijing：National Marine Environmental Forecasting Center，2011（in Chinese）

［2］Reynolds R W，Smith T M．Improved global sea surface temperature analyses using optimum interpolation［J］．J．Climate，1994，7：929-948

［3］馮天祥，李世宏．矩陣QR 分解［J］．西南民族學院學報（自然科學版），2001，27（4）：418-422

Feng Tianxiang，Li Shihong．QR factorization of the metrix［J］．Journal of Southwest University for Nationalities（Natural Science Edition），2001，27（4）：418-422（in Chinese）

［4］馬寨璞．海洋流場數(shù)據(jù)同化方法與應用的研究［D］．杭州：浙江大學，2002

Ma Zhaipu．Study of data assimilation method and its application to sea flow field［D］．Hangzhou：Doctoral dissertation of Zhejiang university，2002（in Chinese）

［5］張識．力法中病態(tài)方程的產(chǎn)生與改進方法［J］．力學與實踐，1989，11（4）：57-59

Zhang Shi．Generation and improving method of the morbid equation in the force method［J］．Mechanics in Engineering，1989，11（4）：57-59（in Chinese）

［6］Matheron G．Principles of geostatistics［J］．Economic Geology，1963，58：1246-1266

［7］Pan Jiayi，Yan XiaoHai，Zheng Quanan，et al．Observation of western boundary current atmospheric convergence zones using scatterometer winds［J］．Geophys．Res．Lett．，2002，29（17）：1832

［8］Chelton D B，Schlax M G，F(xiàn)reilich M H，et al．Satellite measurements reveal persistent small-scale features in ocean winds［J］．Science，2004，303（5660）：978-983

［9］Chelton D B，Xie S P．Coupled ocean-atmosphere interaction at oceanic mesoscales［J］．Oceanography，2010，23（4）：52-69

［10］Xu Haiming，Xu Mimi，Xie Shangping，et al．Deep atmospheric response to the spring Kuroshio over the East China Sea［J］．J．Climate，2011，24：4959-4972