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

齊亞琳 林明森

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

1 引言

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

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

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

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

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

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

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

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

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

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

式中：H為插值算子。

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

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

權重系數矩陣K的表達形式為

式中：B為背景場誤差協方差矩陣，R為觀測誤差協方差矩陣。欲計算權重系數K，必須先估算背景場誤差協方差矩陣B和觀測誤差協方差矩陣R。

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

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

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

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

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

2．1．3 矩陣的求解

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

對式（7）兩邊再同時進行矩陣轉置，得

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

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

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

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

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

從圖2可以發(fā)現，融合后的全球海面溫度，不再存在由于軌道覆蓋范圍不夠導致的軌道間空隙，而且全球海面溫度的空間分辨率不降低。

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

圖1 HY-2衛(wèi)星掃描微波輻射計數據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)星掃描微波輻射計數據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)星掃描微波輻射計數據西北太平洋區(qū)域SST 分布Fig．3 SST of the west north Pacific from the fused HY-2radiometer data

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

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

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

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

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

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

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

1）求變異函數

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

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

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

2）利用變異函數求權重系數

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

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

其中

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

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

圖4和圖5 分別是同時間內的HY-2 衛(wèi)星微波散射計風場和美國國家環(huán)境預報中心（NCEP）再分析數據。該軌道的HY-2 衛(wèi)星微波數據不能覆蓋東海和南海，而NCEP 數據可覆蓋全球范圍內的海域。故在實際應用中需將高分辨率的HY-2風場數據和不受任何軌道等影響的NCEP數據進行融合。圖6是圖4 和圖5 所示的HY-2 微波散射計風場數據和NCEP數據的融合結果。圖6 顯示的融合產品中，風場的空間分辨率達到與HY-2風場空間分辨率相同，并且覆蓋東海和南海等海域。HY-2衛(wèi)星微波散射計一天最多可覆蓋全球90%的海域，而融合后的產品可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)星遙感數據不能全覆蓋的不足，同時提高了NCEP 數據的空間分辨率，可為兩種數據更好地使用，尤其是為業(yè)務化海洋監(jiān)測應用提供便利。

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

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

4 結束語

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

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

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

（References）

［1］奚萌．基于最優(yōu)插值算法的紅外和微波遙感海面溫表溫度數據融合［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］馬寨璞．海洋流場數據同化方法與應用的研究［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)方程的產生與改進方法［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，Freilich 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