多源海洋遙感葉綠素?cái)?shù)據(jù)融合

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

(1.中國(guó)石油大學(xué)(華東) 地球科學(xué)與技術(shù)學(xué)院，山東 青島 266580；2.海洋科學(xué)與技術(shù)國(guó)家實(shí)驗(yàn)室 海洋礦產(chǎn)資源評(píng)價(jià)與探測(cè)技術(shù)功能實(shí)驗(yàn)室，山東 青島 266071)

0 引言

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

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

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

通過(guò)研究分析，國(guó)內(nèi)外相關(guān)科研機(jī)構(gòu)生產(chǎn)的海洋葉綠素α融合產(chǎn)品存在精度低、覆蓋率低、時(shí)間跨度短等問(wèn)題。本文基于前人研究成果，利用SeaWIFS、Terra-MODIS、Aqua-MODIS、MERIS和VIIRS共5個(gè)傳感器的葉綠素α數(shù)據(jù)，實(shí)現(xiàn)了小波變換與Kalman濾波技術(shù)相結(jié)合的融合方法，構(gòu)建了1998—2017年全球海洋表面的葉綠素α數(shù)據(jù)集，并與實(shí)測(cè)數(shù)據(jù)進(jìn)行了對(duì)比分析。

1 數(shù)據(jù)獲取與預(yù)處理

1.1 多源遙感數(shù)據(jù)

本文使用的遙感數(shù)據(jù)為SeaWIFS、Terra-MODIS、Aqua-MODIS、MERIS和VIIRS共5個(gè)傳感器的葉綠素α濃度數(shù)據(jù)，時(shí)間范圍為1998—2017年，數(shù)據(jù)產(chǎn)品為L(zhǎng)evel 3標(biāo)準(zhǔn)網(wǎng)格數(shù)據(jù)，數(shù)據(jù)格式為Netcdf(network common data form)和HDF4(hierarchical data format)，空間分辨率為4 km和9 km，如表1所示。相關(guān)數(shù)據(jù)均可從美國(guó)NASA的oceancolor網(wǎng)站免費(fèi)下載。

表1 5個(gè)傳感器葉綠素α濃度數(shù)據(jù)參數(shù)

因5顆衛(wèi)星發(fā)射時(shí)間不同，因此各個(gè)衛(wèi)星數(shù)據(jù)的時(shí)間覆蓋范圍也不同?？紤]到融合產(chǎn)品的質(zhì)量和精度，需要對(duì)同一時(shí)間段內(nèi)的數(shù)據(jù)進(jìn)行融合。

1.2 實(shí)測(cè)數(shù)據(jù)

實(shí)測(cè)葉綠素α濃度數(shù)據(jù)用于檢驗(yàn)融合產(chǎn)品的精度和可靠性。實(shí)測(cè)數(shù)據(jù)共獲取了5 710個(gè)觀測(cè)文件。觀測(cè)文件記錄了站點(diǎn)編號(hào)、日期、時(shí)間、緯度、經(jīng)度、水深、葉綠素濃度等，數(shù)據(jù)量共1.53 GB。下載網(wǎng)址為http：//seabass.gsfc.nasa.gov。

1.3 數(shù)據(jù)預(yù)處理

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

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

2 數(shù)據(jù)融合方法研究

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

經(jīng)數(shù)據(jù)預(yù)處理，可得到進(jìn)行融合處理的標(biāo)準(zhǔn)圖像。根據(jù)時(shí)間范圍，將其劃分為6個(gè)時(shí)間段進(jìn)行融合，如表2所示。

表2 融合時(shí)間分配表

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

2.1 融合權(quán)重確定

不同數(shù)據(jù)源的數(shù)據(jù)精度不同，融合時(shí)會(huì)賦予不同的權(quán)重值。融合系數(shù)計(jì)算原理是：數(shù)據(jù)源對(duì)應(yīng)的誤差越小權(quán)重越大；一致性測(cè)度數(shù)值越大，權(quán)重越大。結(jié)合這2個(gè)信息，計(jì)算一致性測(cè)度，來(lái)確定低頻系數(shù)矩陣的融合權(quán)重值。

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

(1)

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

(2)

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

(3)

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

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

(4)

(5)

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

(6)

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

(7)

2.2 小波變換融合與Kalman濾波技術(shù)

對(duì)數(shù)據(jù)源計(jì)算權(quán)重后，利用小波變換對(duì)數(shù)據(jù)進(jìn)行加權(quán)融合。融合算法對(duì)2個(gè)或多個(gè)數(shù)據(jù)源進(jìn)行小波變換。首先采用小波多分辨率分解和Mallat快速算法，將原始圖像分解成近似圖像和細(xì)節(jié)圖像，其分別代表了圖像的不同結(jié)構(gòu)；然后在各層的特征域上進(jìn)行有針對(duì)性的融合。因?yàn)楸容^容易提取原始圖像的結(jié)構(gòu)信息和細(xì)節(jié)信息，所以融合效果較好，并且由于小波變換具有完善的重建能力，故保證了信號(hào)在分解重建過(guò)程中，很少有信息損失和信息冗余產(chǎn)生[14]。

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

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

3 結(jié)果與討論

3.1 全球葉綠素α濃度融合產(chǎn)品的分析

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

1)均值分析。均值可以反映整體變化情況。通過(guò)分析全球葉綠素α融合產(chǎn)品均值(圖1)，可以發(fā)現(xiàn)融合產(chǎn)品的均值變化有一定的規(guī)律。每年的圖像都是46幅，圖中剛好有5個(gè)波峰，說(shuō)明一年的數(shù)據(jù)中葉綠素值總會(huì)達(dá)到峰值，峰值出現(xiàn)的時(shí)間在一年的中間，在年初和年末的產(chǎn)品中均值較低。3組融合產(chǎn)品的均值曲線變化規(guī)律相似，峰谷時(shí)間大致吻合，說(shuō)明融合產(chǎn)品在時(shí)間尺度上有很大的一致性和可靠性。

圖1 多年融合產(chǎn)品均值對(duì)比

2)方差分析。通過(guò)分析全球葉綠素α融合產(chǎn)品方差(圖2)，可知融合產(chǎn)品的方差曲線擬合性很好，曲線變化趨勢(shì)一致，說(shuō)明融合產(chǎn)品質(zhì)量穩(wěn)定，沒(méi)有缺失數(shù)據(jù)存在，其中方差最小的是第1組數(shù)據(jù)，說(shuō)明該組數(shù)據(jù)穩(wěn)定性最好。

圖2 多年融合產(chǎn)品方差對(duì)比

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

圖3 多年融合產(chǎn)品信息量對(duì)比

3.2 全球葉綠素α濃度融合產(chǎn)品的驗(yàn)證

1)融合產(chǎn)品與實(shí)測(cè)值的對(duì)比分析。從實(shí)測(cè)數(shù)據(jù)集中選擇每8 d最大值和對(duì)應(yīng)的坐標(biāo)，從融合產(chǎn)品中對(duì)應(yīng)選擇相同坐標(biāo)和相同時(shí)間的葉綠素值，將2組數(shù)據(jù)繪制成圖(圖4)。從圖4可知，二者的點(diǎn)匹配度較高，超過(guò)60%的點(diǎn)擬合度很好。

圖4 融合產(chǎn)品與實(shí)測(cè)值相關(guān)性

2)融合產(chǎn)品與已有產(chǎn)品的對(duì)比分析。

①可利用率對(duì)比分析。將2008年的融合產(chǎn)品與歐空局對(duì)應(yīng)時(shí)間的GSM葉綠素產(chǎn)品進(jìn)行可利用率分析，如圖5所示?？芍?，融合產(chǎn)品的可利用率高于GSM產(chǎn)品，以編號(hào)9產(chǎn)品為例，可利用率差為37.46%，其中22號(hào)產(chǎn)品的可利用率差最小，為28.20%，故本文的融合產(chǎn)品可利用率較已有產(chǎn)品有很大提高，為數(shù)據(jù)的更深層次應(yīng)用提供了充分的保障。

圖5 融合產(chǎn)品與GSM產(chǎn)品可利用率對(duì)比

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

圖6 融合產(chǎn)品和GSM產(chǎn)品、實(shí)測(cè)值對(duì)比

選取2008年的融合產(chǎn)品、GSM產(chǎn)品與實(shí)測(cè)最大值對(duì)比，共匹配了24對(duì)點(diǎn)，分別以實(shí)測(cè)-融合、實(shí)測(cè)-GSM作為橫縱坐標(biāo)值畫出相關(guān)性圖，如圖7所示。根據(jù)數(shù)據(jù)的相關(guān)性，可知融合產(chǎn)品與實(shí)測(cè)值的相關(guān)性為0.792 2，GSM與實(shí)測(cè)值的相關(guān)性為0.349 4，證明本文融合產(chǎn)品質(zhì)量較高，與實(shí)測(cè)點(diǎn)匹配較好。

圖7 融合產(chǎn)品、GSM產(chǎn)品與實(shí)測(cè)值相關(guān)性對(duì)比

4 結(jié)束語(yǔ)

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

本文中因有部分圖像數(shù)據(jù)存在缺值現(xiàn)象，并沒(méi)有在融合之前進(jìn)行補(bǔ)值處理，因此影響了融合產(chǎn)品的質(zhì)量，對(duì)融合產(chǎn)品的覆蓋率和精度有一定的影響。另外，融合產(chǎn)品只與歐空局的GSM產(chǎn)品進(jìn)行了對(duì)比，沒(méi)有在實(shí)際項(xiàng)目中應(yīng)用，對(duì)實(shí)際工作的貢獻(xiàn)有待后續(xù)考察。