多源海洋遙感葉綠素數據融合

崔建勇，劉曉東，岳增友，李連偉

(1.中國石油大學(華東) 地球科學與技術學院，山東 青島 266580；2.海洋科學與技術國家實驗室 海洋礦產資源評價與探測技術功能實驗室，山東 青島 266071)

0 引言

地球表面70%的面積被海洋所覆蓋，海洋為人類生活提供了豐富的魚類資源。隨人類社會技術的不斷進步，特別是遙感技術的快速發展，人類與海洋已融合為一個生命共同體。1978年美國宇航局(National Aeronautics and Space Administration，NASA)發射了Nimbus-7(雨云)號衛星，揭開了利用遙感技術探測海洋的新篇章。之后，NASA和其他各國衛星發射機構相繼發射了多顆海洋衛星，實現了大區域、高頻率、高精度的海洋觀測，為科研人員研究海洋提供了海量遙感數據，對研究海洋生態環境演變、生物地球化學循環、氣候變化、漁業捕撈等具有重要的意義[1-2]，也為探測海洋葉綠素α提供了堅實的數據基礎。

不同的海洋衛星遙感傳感器具有不同的空間分辨率和時間分辨率，光譜和覆蓋率也各有不同，故單個海洋葉綠素α濃度數據源存在數據覆蓋率、分辨率和數據可利用率等方面的缺陷。為彌補該缺陷，科研人員提出了多源數據融合技術。多源數據融合技術可實現多遙感數據信息互補，可擴展單個海洋遙感的空間覆蓋率和時間分辨率[3]，增加整體的信息量，提高數據產品的時空連續性、一致性和可靠性。

多源數據融合技術最早被廣泛應用于海表溫度、衛星高度計測高數據等方面，隨海洋遙感衛星的發展，該技術被用于海色數據的融合。1997年，NASA正式提出了SIMBIOS計劃，同時歐空局(European Space Agency，ESA)提出了GLOBcolour計劃[4]，相關計劃大多使用平均法、GSM(Garver-Siegel-Maritorena)生物光學模型等方法，將多顆海色傳感器的數據進行融合，建立了高質量、高精度、高覆蓋率的海色數據集。李新星等[1]指出，在2005年GSM 方法提出后，Maritorena等人在2010年將該方法用于SeaWIFS、MODIS、MERIS的數據融合。國內學者在數據融合方面稍落后于國外學者，較早的研究是張靈凱等[5-6]用小波分析方法對SeaWIFS和MODIS的葉綠素α數據進行融合。

通過研究分析，國內外相關科研機構生產的海洋葉綠素α融合產品存在精度低、覆蓋率低、時間跨度短等問題。本文基于前人研究成果，利用SeaWIFS、Terra-MODIS、Aqua-MODIS、MERIS和VIIRS共5個傳感器的葉綠素α數據，實現了小波變換與Kalman濾波技術相結合的融合方法，構建了1998—2017年全球海洋表面的葉綠素α數據集，并與實測數據進行了對比分析。

1 數據獲取與預處理

1.1 多源遙感數據

本文使用的遙感數據為SeaWIFS、Terra-MODIS、Aqua-MODIS、MERIS和VIIRS共5個傳感器的葉綠素α濃度數據，時間范圍為1998—2017年，數據產品為Level 3標準網格數據，數據格式為Netcdf(network common data form)和HDF4(hierarchical data format)，空間分辨率為4 km和9 km，如表1所示。相關數據均可從美國NASA的oceancolor網站免費下載。

表1 5個傳感器葉綠素α濃度數據參數

因5顆衛星發射時間不同，因此各個衛星數據的時間覆蓋范圍也不同。考慮到融合產品的質量和精度，需要對同一時間段內的數據進行融合。

1.2 實測數據

實測葉綠素α濃度數據用于檢驗融合產品的精度和可靠性。實測數據共獲取了5 710個觀測文件。觀測文件記錄了站點編號、日期、時間、緯度、經度、水深、葉綠素濃度等，數據量共1.53 GB。下載網址為http：//seabass.gsfc.nasa.gov。

1.3 數據預處理

為方便后期融合算法的流程化計算，需要將5個傳感器的數據轉換成相同的格式和相同的空間分辨率。SeaWIFS、Terra-MODIS、Aqua-MODIS和VIIRS傳感器的葉綠素α濃度數據格式為Netcdf格式，而MEIRS傳感器的葉綠素α濃度數據是HDF4格式。空間分辨率方面也存在差異，SeaWIFS傳感器的葉綠素α濃度數據的空間分辨率是9 km，其他4個傳感器的葉綠素α濃度數據空間分辨率是4 km。

為解決上述2項問題，采用ENVI-IDL語言和Python語言編寫數據預處理程序，分別對葉綠素α濃度數據進行處理。處理完成后，5個傳感器的葉綠素α濃度數據的格式均為TIFF，空間分辨率均為4 km，地理空間投影一致。

2 數據融合方法研究

融合算法以小波變換為基礎，并結合Kalman濾波技術。小波變換算法中正交或雙正交多分辨率小波融合方法已用于多傳感器圖像信息的融合，該方法最大限度地保留了原多光譜圖像的光譜信息[7-8]，并且能將圖像分解為具有不同空間分辨率、頻率特性和方向特性的子圖像序列。其分頻功能相當于一對共軛鏡像正交濾波器組，把圖像解構為1個低頻帶子圖和3個高頻帶子圖(包括水平、垂直和對角線3個方向)[9-11]。低頻部分反映的是圖像的整體視覺信息，高頻部分反映的是圖像的細節特征，具有多尺度、多分辨率、方向性和無冗余數據等優點[12-13]。Kalman濾波是一種基于協方差遞歸思想的預測器，能夠尋求最優的預測結果，效率較高，將其應用于對高頻信息的處理，可以生成質量更好的高頻信息。

經數據預處理，可得到進行融合處理的標準圖像。根據時間范圍，將其劃分為6個時間段進行融合，如表2所示。

表2 融合時間分配表

從表2可知，1998—1999年遙感數據只有SeaWIFS一種數據源，無法采用本文提出的融合方法進行處理，可采用查找表法進行融合；2000—2017年均有多個數據源，根據數據源個數的不同，研制不同的小波變換融合算法。

2.1 融合權重確定

不同數據源的數據精度不同，融合時會賦予不同的權重值。融合系數計算原理是：數據源對應的誤差越小權重越大；一致性測度數值越大，權重越大。結合這2個信息，計算一致性測度，來確定低頻系數矩陣的融合權重值。

假設有n個傳感器，這些傳感器獲取同一觀測范圍的地物屬性。設傳感器i在K時刻的采集值為Xi(k)(k=1，2，…，s；k表示像素編號；s為觀測范圍像素個數)，傳感器j在K時刻的采集值為Xj(k)，時間為K時，傳感器i與j在k位置測量值的支撐力度矩陣如式(1)所示。

(1)

時間為K時，像素位置k的各傳感器觀測結果總支持度如式(2)所示。

(2)

總支持度矩陣中每一行的支持度之和，代表了對應傳感器對最終觀測結果的支持度。因此，時間為K時，像素位置k各傳感器觀測結果支持度如式(3)所示。

(3)

式中：一致性測度ri(k)僅反映了在k像素位置，傳感器i的測量值與所有傳感器測量值的相近水平。單個像素ri(k)的值，并不代表在整個觀測區域上傳感器的可靠性，應考慮到整個觀測范圍上傳感器節點的可靠性。

第K時刻第i個傳感器觀測支持度的均值和方差分別如式(4)、式(5)所示。

(4)

(5)

若某傳感器觀測范圍內像素平均值較大，且灰度的波動較小，則該傳感器較穩定，具有較高的支持度，權重較大。因此，可用信噪比描述此傳感器的支持度，信噪比的值為平均值除以信號波動的方差。時間為K時，傳感器i的權重wi(K)如式(6)所示。

(6)

歸一化后值如式(8)所示。

(7)

2.2 小波變換融合與Kalman濾波技術

對數據源計算權重后，利用小波變換對數據進行加權融合。融合算法對2個或多個數據源進行小波變換。首先采用小波多分辨率分解和Mallat快速算法，將原始圖像分解成近似圖像和細節圖像，其分別代表了圖像的不同結構；然后在各層的特征域上進行有針對性的融合。因為比較容易提取原始圖像的結構信息和細節信息，所以融合效果較好，并且由于小波變換具有完善的重建能力，故保證了信號在分解重建過程中，很少有信息損失和信息冗余產生[14]。

小波變換需要指定小波基函數，本文采用Haar小波基。Haar函數是小波分析中最早用到的一個具有緊支撐正交小波函數。經過小波變換后，圖像被分解為4個維度的子圖像，由1個低頻信息、3個高頻信息組成，圖像分辨率降為原來圖像的一半。首先對每一個數據源小波變換的低頻信息按權重進行加權處理，權重值根據式(7)計算得到，經過加權處理后的信息直接作為小波反變換的低頻信息參與融合。接下來對高頻信息進行Kalman濾波處理。

首先對n(n為傳感器數量)個數據源的高頻信息進行Kalman濾波計算，生成3×n對高頻信息；然后按照取最大值的原則，對n組高頻信息進行篩選，最終由一個加權運算得到的低頻信息和3個新的高頻信息進行小波反變換，實現多源圖像加權小波融合過程，算法對每幅圖像都作上述處理；最后得到相應的8 d融合產品，再根據8 d融合產品進行月季年產品的合成。

3 結果與討論

3.1 全球葉綠素α濃度融合產品的分析

全球葉綠素α濃度融合產品可從數據均值、方差、信息量3個方面進行分析評價。因融合產品時間跨度大，為能充分反映不同時間段內融合數據的可利用率，本文將數據分3組進行對比分析：第1組數據的時間段為1998—2002年；第2組數據的時間段為2005—2009年；第3組數據的時間段為2012—2016年。

1)均值分析。均值可以反映整體變化情況。通過分析全球葉綠素α融合產品均值(圖1)，可以發現融合產品的均值變化有一定的規律。每年的圖像都是46幅，圖中剛好有5個波峰，說明一年的數據中葉綠素值總會達到峰值，峰值出現的時間在一年的中間，在年初和年末的產品中均值較低。3組融合產品的均值曲線變化規律相似，峰谷時間大致吻合，說明融合產品在時間尺度上有很大的一致性和可靠性。

圖1 多年融合產品均值對比

2)方差分析。通過分析全球葉綠素α融合產品方差(圖2)，可知融合產品的方差曲線擬合性很好，曲線變化趨勢一致，說明融合產品質量穩定，沒有缺失數據存在，其中方差最小的是第1組數據，說明該組數據穩定性最好。

圖2 多年融合產品方差對比

3)信息量分析。信息量反映的是圖像本身所含信息的豐富程度。對融合產品信息量進行分析(圖3)，3組數據中信息量的曲線變化趨勢一致，說明在全年或者長時間的數據中，信息量的變化是一致的，在特定的時間信息量會增大，而在其他時間就會減少。信息量最大的是第3組數據，其次是第2組數據，最后是第1組數據。

圖3 多年融合產品信息量對比

3.2 全球葉綠素α濃度融合產品的驗證

1)融合產品與實測值的對比分析。從實測數據集中選擇每8 d最大值和對應的坐標，從融合產品中對應選擇相同坐標和相同時間的葉綠素值，將2組數據繪制成圖(圖4)。從圖4可知，二者的點匹配度較高，超過60%的點擬合度很好。

圖4 融合產品與實測值相關性

2)融合產品與已有產品的對比分析。

①可利用率對比分析。將2008年的融合產品與歐空局對應時間的GSM葉綠素產品進行可利用率分析，如圖5所示。可知，融合產品的可利用率高于GSM產品，以編號9產品為例，可利用率差為37.46%，其中22號產品的可利用率差最小，為28.20%，故本文的融合產品可利用率較已有產品有很大提高，為數據的更深層次應用提供了充分的保障。

圖5 融合產品與GSM產品可利用率對比

②精度對比分析。將本文融合產品與2008年的實測值、歐空局的GSM產品進行對比分析，如圖6所示。可知，本文融合產品、GSM產品與實測值擬合性均較好，但有部分點葉綠素濃度值相差較大，融合產品與GSM產品的葉綠素濃度值均高于實測值。

圖6 融合產品和GSM產品、實測值對比

選取2008年的融合產品、GSM產品與實測最大值對比，共匹配了24對點，分別以實測-融合、實測-GSM作為橫縱坐標值畫出相關性圖，如圖7所示。根據數據的相關性，可知融合產品與實測值的相關性為0.792 2，GSM與實測值的相關性為0.349 4，證明本文融合產品質量較高，與實測點匹配較好。

圖7 融合產品、GSM產品與實測值相關性對比

4 結束語

由于在軌運行的海色衛星大多是極軌衛星，葉綠素α濃度產品受到空間覆蓋率和空間分辨率不同程度的影響，從而限制了衛星產品的實際使用。本文提出了基于小波變換和Kalman濾波技術的融合算法，將SeaWIFS、Terra-MODIS、Aqua-MODIS、MERIS和VIIRS共5個傳感器的葉綠素α濃度數據進行融合，完成了融合產品的均值、方差和信息量的分析，并完成了融合產品與實測值、歐空局的GSM產品的對比分析。通過對比分析可知，本文的融合產品在空間覆蓋率、數據可利用率、葉綠素α濃度值的精度方面均優于歐空局的GSM產品。本文構建的小波變換方法與Kalman濾波技術相結合的融合算法有效提高了融合產品的質量，在未來海色數據的應用方面具有較大的潛力。

本文中因有部分圖像數據存在缺值現象，并沒有在融合之前進行補值處理，因此影響了融合產品的質量，對融合產品的覆蓋率和精度有一定的影響。另外，融合產品只與歐空局的GSM產品進行了對比，沒有在實際項目中應用，對實際工作的貢獻有待后續考察。