不同時間尺度多源時序數據的FEEMD分解比較研究

王 正, 邱士可, 曾 群, 呂言利, 王 超, 張起萍, 李雙權

(1.河南省科學院地理研究所, 鄭州 450052; 2.河南省遙感與GIS重點實驗室,鄭州 450052;3.華中師范大學城市與環境科學學院, 武漢 430079; 4.北京國遙新天地信息技術股份有限公司, 北京 100083)

眾多的海洋生態學和動力學研究表明,海表葉綠素a濃度的變化是浮游植物對多個物理環境因素微弱變動信號的非線性放大響應[1].因此,表征浮游植物豐度的海表葉綠素a濃度與多個環境因子間存在密切的內在關系[2-3].研究影響長時間序列海表葉綠素a濃度的環境因子[4]發現,太陽有效光合輻射(photosynthetically active radiation,PAR )和海面風場影響海洋上層環境;海表溫度(sea surface temperature, SST)、海面高度(sea surface height, SSH)、海表鹽度(sea surface salinity, SSS)、海表密度(sea surface density, SSD)、水體的分層結構(用混合層深度(mixed layer depth, MLD)來表示)表征海洋上層的熱力和動力狀況[5];基于海表風場的風速(wind speed)和風應力(wind stress)計算的風應力旋度和?？寺槲俾蕝⒘?表征浮游植物生長的溫度和營養鹽供給,與葉綠素濃度關系密切[6-12].這些環境參量都是影響海表葉綠素a濃度的關鍵因子.

中國南海東北部海域受到黑潮、閩浙沿岸流、呂宋西北部上升流、內波、季風、河流輸入等多個因素的影響,位置特殊,水文環境復雜,水體動態度較高,該區域的浮游植物和相關多個環境因子也呈現出高動態度特征.因此, 對該區域21 a(1998—2018年)的葉綠素a濃度及多個相關環境因子8 d和月尺度數據進行分解面臨較多挑戰,迫切需要一種適合該區域特征的數據分解方法,以分析葉綠素a濃度和多個環境因子蘊含的內在機制和規律.

目前,常見的時序數據分解方法有經驗正交函數分解法(empirical orthogonal function,EOF),傅里葉變換, 小波變換等[13-14],然而,使用這些方法需要滿足一定的前提條件[15-16].傅里葉變換需要假定信號按照固定的時間完全重復,EOF變換極度敏感于所選數據的時間與空間范圍;而小波變換需要先確定基函數,再以基函數為基數對時間序列數據進行分解.海洋葉綠素濃度的變化是對海洋物理環境因子變化的非線性多尺度響應,多個環境因子數據量大、影響因子多樣、變化比較劇烈,因此如果利用常用的方法進行時序數據分解具有一定的局限性,需要用非線性自適應的方法來進行數據分析.經驗模特分解法(empirical mode decomposition, EMD)以及基于EMD發展的集合經驗模態分解法(ensemble empirical mode decomposition, EEMD)、互補集合經驗模態分解CEEMD(complementary ensemble empirical mode decomposition, CEEMD)等系列方法能將頻率不同的不規則時間序列的信號按照不同頻率進行分解[17],是一類自適應的數據處理方法,非常適合非線性、非平穩態和多尺度時間序列數據的處理,已在數據預測、信號分解領域得到推廣, 尚未在遙感領域廣泛應用[18-21].

因此,本文選用南海東北部區域1998—2018年共21 a的8 d和月尺度: SST、 海表面高度異常數據(sea level anomaly, SLA)、 SSS、MLD、 wind stress、wind speed、PAR、SSD這些物理要素與葉綠素a濃度數據進行FEEMD分解和比較,以分析和探討其內在的機制關系.

1 數據及方法

1.1 研究區域

南海位于中國大陸的南方,是太平洋西部海域,夏半年5—9 月盛行西南季風, 臺風頻繁;冬半年盛行東北季風,季風穩定且強度很大.同時,南海既處于太平洋、印度洋和青藏高原的交互區域,又位于印度洋和太平洋兩個沃克環流的交匯點上,還受黑潮入侵、閩浙沿岸流、珠江和湄公河等多條大河輸入的影響.多尺度強迫影響下南海處于氣候變化的敏感區域,各種尺度的海洋動力過程如中尺度渦旋和上升流頻發.南海東北部海區(圖1的NE區域)除了受南海多尺度多因素影響外,更位于太平洋和南海水體交換的主要通道,許多研究表明該區域常年有冷渦、上升流和大規模浮游植物葉綠素a濃度極值現象的發生[22-25].因此,本文在南海東北部海洋環境狀況極端復雜的NE區域開展葉綠素a濃度和多個密切相關環境因子間的長時序數據分解研究.

1.2 研究數據及來源

本文使用的L3級8 d和月尺度合成葉綠素a濃度ChlOC5數據自Global Color免費獲取 (表1), Global Color提供了4 km、9 km和25 km 的產品,包括 MERIS、MODIS、SeaWiFS 、VIIRS和OLCI以及這些多傳感器合成的產品,數據空間覆蓋范圍為全球,時間范圍從1997年9月1日—2023年7月(本文成稿時間).SST數據由美國國家氣候數據中心(National Climatic Data Center,NCDC)提供,該數據是經過最優插值法插值的AVHRR4級的產品,時間分辨率為每天,空間分辨率為0.25度,數據的時間尺度跨度自1981年9月至今.法國AVISO數據中心提供了融合多源衛星測高數據的SSH和海面高度異常數據(sea surface anomaly, SLA),該數據是TOPEX/Poseidon、ERS-1/2、ENVISAT Jason-1和Geosat Follow-On這4種測高衛星經交叉定標的網格化海平面高度融合數據.ASCAT風場數據為插值后的8 d和月尺度的全球風場數據,空間分辨率為0.25度×0.25度,參數包括海表面的wind speed和 wind stress.MLD數據是基于閾值法計算的, 使用了周和月尺度混合層再處理數據, 這兩種再處理數據的空間分辨率為25 km, 數據共有33層,本文所用的表層數據是從1998年1月1日到2018年12月31日的周和月尺度數據.SSD數據通過微波遙感技術間接反演獲得,再結合實測的數據和數值模擬數據同化而來的,SSD數據與SSS數據都下載自法國CMEMS. PAR數據下載自歐空局的Global Color計劃所提供的長時間序列的多傳感器融合數據,參數包括了生地化參數,海洋表層光學參數、大氣參數和表觀光學參數;本文使用的是從1998年到2018全年8 d和月尺度時間分辨率的融合數據.本研究所用數據及下載地址見表1.

表1 本文所用數據及下載地址(1998-01—2018-12)

1.3 EMD、 EEMD FEEMD數據分解法

經驗模態分解方法是一種具有良好自適應數據處理能力的方法,非常適合非線性、多尺度非平穩態的長時間序列數據的處理,該方法的本質是對數據序列代表的信號進行平穩化處理[26].其基本思想是將時間序列數據信號看作是頻率不同的不規則波, 并按不規則波的不同頻率將時序數據分解為多個單一頻率的波和剩余殘波的形式.其基本步驟為:

1) 求取時間序列信號x(t)的極值點,包括極小值點和極大值點;

2) 利用三次樣條函數構造時間序列數據極大值點和極小值點的上下包絡線,并計算其均值函數m1;

3) 驗證h1=x(t)-m1是否符合本征模態函數(intrinsic mode function, IMF)的條件,如符合則進行下一步計算,如不符合則重復步驟1和步驟2,直至符合條件為止,結果即為第一個IMF,imf1=h1k;

4) 得到原始時間序列數據扣除第一個IMF后的第一個殘留r1=x(t)-imf1.重復以上步驟,則可得到從高頻到低頻的各模態IMF;

5) 原始時間序列數據信號扣除各IMF,直至剩余的信號最后變為單調信號Res(趨勢),即只存在一個極值點的情況為止.則原始時間序列的信號可表示為:

(1)

其中,每個本征模態都必須滿足數據序列中極值點的個數與過零點的數目相等或相差不超過一個點, 并且時序數據中由局部最小值和局部極大值所構成的上下包絡線的均值為零[26].

然而,雖然EMD方法適用于非線性、非穩定態時間序列數據的分析,但在實際時序數據處理中由于原始信號極值點不均勻分布和極值突變會使得數據丟失一部分的時間或頻率尺度,導致模態混疊問題.模態混疊是指相近的特征時間尺度分布在不同的IMF中, 或者單個IMF包含差異極大的特征時間尺度, 導致相鄰的2個IMF相互影響造成波形混疊難以辨認.當混疊模態出現時,分解的本征模函數(IMF)是沒有實際物理意義的.為了解決模態混疊的問題,吳召華教授等在EMD的基礎上,發展出了集合經驗模態分解法[27].其基本思想是:將頻率平均分布的高斯白噪聲加入到時間序列數據中,并將其看作是信號,則信號具有了在不同尺度上的連續性,可以減小在EMD中存在的模態混疊現象.EEMD分解的步驟與EMD分解的步驟類似:

1) 將均值為0,頻譜正態均勻分布的高斯白噪聲ni(t)多次加入到初始序列數據信號s(t)中,加入的白噪聲的標準差為原始信號標準差的0.1～0.4倍,即

xi(t)=s(t)+ni(t),

(2)

xi(t)為第i次加入的白噪聲的信號;

2) 對序列數據信號xi(t)進行集合經驗模態分解,即可得到分解后的各模態分量imfij和趨勢項res;

3) 通過n次重復步驟1)和2)對相應的各IMF分量取均值,以消除高斯白噪聲對序列數據信號的影響,結果得到的IMF為

(3)

若重復次數越多則n越大,則對應的白噪聲的各個IMF和越趨近于0.因此,可將EEMD分解的結果展示為:

(4)

EEMD通過加入高斯白噪聲的方式,消除EMD中的模態混疊現象,比EMD的分解更科學準確.然而,EEMD中也存在殘余輔助噪聲的影響(圖2c).Yeh等通過將隨機高斯白噪聲以正負成對的方式加入到時序信號中,發展出了互補集合經驗模態分解方法(complementary ensemble empirical mode decomposition, CEEMD)[28],CEEMD就是FEEMD的核心算法.實驗結果表明,FEEMD能有效消除時序數據信號中殘余的輔助噪聲,不僅保留了EEMD在處理非平穩信號方面的優點,又完善了EMD的模態混疊問題,且計算結果表明FEEMD的去噪和抗干擾能力優于EMD和EEMD(圖2d).CEEMD分解的步驟為:

圖2 基于8 d尺度數據的EMD、EEMD和FEEMD分解結果的對比

1) 將n組正負成對的輔助白噪聲加入到原始數據序列信號中.白噪聲的取值范圍為0.1～0.4,即白噪聲的標準差為原始信號的0.1～0.4倍.則可獲得一個包含2n個數據序列信號的集合:

(5)

其中,加入正負成對白噪聲后的時序數據信號為M1和M2,S是原始的時序數據信號,N為加入的輔助白噪聲;

2) 按照EMD分解的步驟依次對數據集合中的每個時序數據信號分解,則每個信號可得一組IMF分量,因此,imfij即為第i個分量和第j個本征模態函數的分量;

3)imfj即為多組分量取均值后得到的分解后結果:

(6)

其中,imfj為FEEMD分解后獲得的第j個IMF分量.FEEMD分解后,產生的各個imfj根據頻率由高到低依次排列,高頻和低頻噪聲分別出現在靠前和靠后的幾個模態分量中,通過顯著模態檢驗可判斷分解出的是噪聲還是具有實際物理意義的信號,分解至最后剩余的信號即為該時序數據的總體變化趨勢.

2 結果

2.1 三種數據分解方法的比較

為比較EMD、EEMD和FEEMD的數據分解方法結果,基于南海北部NE區域8 d合成的21 a時間序列葉綠素a濃度數據進行分解.圖2即為對該時間序列數據(圖2a)進行EMD分解(圖2b)、EEMD分解(圖2c)、FEEMD分解(圖2d)的結果.其中,由于21 a的8 d尺度數據共有966期,則可分解出10個模態變量,其中第1個變量為加噪聲的原始信號,第2～8個分量表示的是從高頻到低頻的各個模態的變量,第9個模態表示的是整個周期時間序列詳細的變化過程,第10個分量表示21 a總體的變化趨勢.

圖2中的b、c、d都能反映出數據的總體趨勢(R),且三者分解出趨勢是一致的,在低頻的第9模態(C9)和第10模態(R),三者是相似的.然而, 圖2b和2c都存在不同程度的模態混疊現象,該現象在C2和C3高頻模態中尤為明顯.其中,C2和C3模態的波形高度相似,此外在EMD分解的結果中(圖2b),模態C4和模態C5在同一模態中出現了其他模態才出現的頻率,這種現象在EEMD分解的結果中有所改善(圖2c).FEEMD的結果(圖2d)則避免了這些問題,同時FEEMD分解的運行效率相較于EMD和EEMD提升了10倍以上.因此,本研究選用FEEMD來進行數據分解和后續處理是合理的.

2.2 研究區域的時序數據的FEEMD分解

研究表明南海東北部區域(NE區域)的葉綠素a濃度存在顯著的周期性和季節性[3].為了對比不同時間尺度數據的分解結果,本文對比分析了NE區域月尺度數據(NE-monthly)來與8 d尺度數據(NE-8d)的分解結果.此外,為了分析影響NE區域顯著周期性現象的關鍵因子,本文對影響海洋葉綠素a濃度的環境因子SST、SSH、MLD、PAR、SSS、SSD、 wind stress和wind speed的8 d和月尺度數據也進行了FEEMD分解.需注意的是,分解出的模態個數由數據的長度來確定,數據長度越長,則分解模態數目越多,可以通過log2n來確定模態數,其中的n為數據的長度.因此,8 d尺度21 a共966期數據可以分解出10個模態,月尺度21 a共有252期數據,則可分解為8個模態.

FEEMD對8 d尺度數據加0.4的噪聲,集合訓練200次;共分解出10個模態,其中第1個模態(C1)為加噪聲的原始信號,模態的性質與噪聲類似,故各變量的第一模態(C1)不參與討論.第9個模態表示的是整個周期時間序列詳細變化過程,第10個模態表示21 a總體變化趨勢,因此在顯著性檢驗時,僅檢驗C2～C8這7個模態(圖4).本文分解并展示了NE-8 d(圖3)和NE-monthly(圖4)的結果.各環境因子C1模態的性質與噪聲類似,因此8 d和月尺度的各變量的C1模態均不參與討論.

2.3 FEEMD數據分解結果的顯著性檢驗

使用IMF平均周期與能量關系分布特征顯著性檢驗來判別FEEMD分解出的各IMF是噪聲還是具有實際物理意義的量.如果分解出的各個模態的結果是在95%的置信線以上,則為具有真實物理意義的參量,在95%的置信線以下,則認為是噪聲.

圖5為NE-8d的FEEMD分解出的各個模態的顯著性檢驗結果,圖5中除了ChlOC5的C8模態和PAR的C6模態未通過顯著性檢驗外,其它各參數的各模態均通過顯著性檢驗,為具有實際物理意義的信號.圖6 NE-monthly為基于月尺度數據分解出的各個模態的顯著性檢驗的結果,SST、MLD、PAR、SSD和wind speed為僅C2模態通過顯著性檢驗的參量;ChlOC5和wind stress為C2和

圖6 NE-Monthly FEEMD分解的各個模態的顯著性檢驗

C3模態通過顯著性檢驗的參量; SLA為C2和C5模態顯著的參量, SSS為C2、C4、C5和C6通過顯著性檢驗的參量.通過對比8 d尺度和月尺度數據FEEMD分解結果的顯著性檢驗可知,8 d尺度數據分解出比月尺度更豐富的、具有實際物理意義的顯著模態.

3 不同時間尺度數據的趨勢、周期對比及分析

已有的長時序數據時空變化研究中,常用月尺度數據進行分析[29-31].8 d尺度數據相較于月尺度數據,有何異同, 能否發現月尺度數據未能發現的現象或規律,是本節也是本文關注的問題之一.

3.1 NE區域8 d與月尺度總體趨勢的比較

圖7和圖8分別為8 d和月尺度葉綠素a濃度與相關環境因子數據在NE區域21 a的變化趨

圖7 基于8 d尺度數據分解出的NE區域各個因子21 a總體趨勢

圖8 基于月尺度數據分解出的21 aNE區域的總體趨勢

勢.從圖7可知在NE區域21 a的葉綠素a濃度的總體趨勢是下降的,與之對應的SST呈上升趨勢,MLD呈先上升后下降的趨勢,NE區域的PAR也呈現下降趨勢,SLA呈先上升后略有下降的趨勢,SSD呈下降趨勢,SSS在1998-2010年呈現下降趨勢,2010-2018年呈現上升趨勢,21 a的wind stress, wind speed 也呈下降的趨勢.對比8 d(圖7)和月尺度(圖8)數據的變化趨勢,發現兩者趨勢一致,表明FEEMD方法具有較好穩定性.

3.2 月尺度與8 d尺度數據分解出的周期對比

為進一步比較8 d尺度與月尺度數據的FEEMD分解結果,利用經過顯著性檢驗的模態來計算顯著模態周期和方差貢獻.計算周期的方法主要有兩種,基于瞬時頻率和基于零點數的方案.在利用瞬時頻率方案計算平均周期時,容易出現異常值.研究結果顯示,出現異常值是因為出現了負的瞬時頻率,零點數方案則可以有效的避免此種情況的發生,因此采用零點數方案來計算顯著模態的平均周期.

圖9a為8 d尺度的葉綠素a濃度數據分解出的顯著模態平均周期,C7模態為平均63個月的周期,C6模態為31.5個月的周期,C4和C5模態為準年周期,C3為季節周期/半年周期,C2為月尺度周期.綜合圖9中的所有因子可知,基于8 d尺度的數據可分解出準年周期的模態(C4模態),也可分解出年際(C5/C6模態)和年代際周期的模態(C8/C7模態),值得注意的是,基于8 d尺度的數據可以分解出季節尺度(C3模態)甚至是月尺度(C2模態)的周期.圖10為各因子顯著模態的方差貢獻,由圖10 (a～i)可知,第四模態的方差貢獻最大,對于ChlOC5因子C2～C4模態累計貢獻了所有模態的總方差的83.4%左右(圖10a),SST的C2～C4模態累積貢獻了總方差的97%左右(圖10b),MLD的C2～C4模態累計貢獻了總方差的89.2%左右(圖10c), PAR的C2～C4模態累計貢獻了總方差的95.6%左右(圖10d); SLA (圖10e), SSD (圖10f), SSS (圖10g), wind stress (圖10h) 和wind speed (圖10i)這些因子的C2～C4模態累計的方差貢獻分別為79%,96%,63%,85%,92.5%左右.

圖9 區域NE-8 d各環境因子的顯著模態的周期

圖10 區域NE-8 d尺度各環境因子顯著模態的方差貢獻/%

綜合圖9和圖10來看,各個模態中,方差貢獻最大的是C4模態,C4模態主要體現的是年尺度的變化趨勢,C4模態的方差貢獻在57%以上,最高可達90%左右.值得注意的是,除了方差貢獻最大的C4模態,方差貢獻第二和第三的模態周期一般為C3和C2模態的周期(圖10),個別的因子除外.C3和C2模態分別對應于季節周期和月尺度的周期.圖9所有因子的C2模態反映的周期從1.8～2.1個月不等,但主要是1.8個月左右的短時間頻率的周期,圖9的C3模態反映的是從3.7～5.86個月的季節尺度和半年尺度的周期.

表2是NE區域基于月尺度的各個因子的數據,在FEEMD分解后,經過顯著性檢驗的模態所計算的平均周期,方差貢獻和累積方差表.可以看出,NE區域的月尺度數據FEEMD分解結果的能量分布主要集中在C2模態上,NE區域的ChlOC5-monthly的C2模態方差貢獻達到70.5%,顯著的C2和C3模態的方差貢獻達到了93%,ChlOC5-monthly主要為年尺度周期.此外,SST-monthly,MLD-monthly, PAR-monthly, SSD-monthly, wind speed-monthly通過顯著性檢驗的模態為C2,周期分別為12、12、11.7、12.3、11.5個月,C2模態的方差貢獻分別為94.7%、86%、86.4、95.2%、78.4%.除了這些環境因子外,還有其他環境因子有兩個及以上模態通過顯著性檢驗,SLA-monthly的C2和C5模態通過了顯著性檢驗,C2模態反映了顯著的年周期變化情況,C5模態反映了顯著的168個月的年代際的變化情況,C2模態的方差貢獻為81.6%,C2和C5模態累積的方差貢獻為87.9% .SSS-Monthly的顯著模態有C2、C4、C5和C6模態,周期分別為12.3、45.8、100.8和168個月,顯著模態的方差貢獻最大的是C2模態,為58.2%,幾個顯著模態累積的方差貢獻為87.3%.wind stress-monthly的顯著模態C2和C3的周期可看作準年周期,累計方差貢獻達到92.1%.NE-monthly的ChlOC5和各環境因子分解結果可以反映出該區域顯著的年周期變化情況以及年際的變化情況.

綜合來看,NE區域的8 d和月尺度數據的FEEMD分解結果都可將該區域最顯著的年尺度周期信息表征出來,此外兩者都分解出年際和年代際的周期.然而兩者又有不同,8 d尺度數據的結果,可分解出顯著的半年尺度的變化周期、季節變化的周期和月尺度變化周期,這是月尺度的數據所無法做到的.另外,基于FEEMD分解的8 d尺度數據能分解出從短到長的不同周期,而月數據的FEEMD分解結果主要以準年周期為主.

4 討論

本文利用FEEMD方法分解了月尺度和8 d尺度數據,通過比較分析兩者具有實際物理意義的顯著模態、模態周期、模態方差貢獻等發現,8 d尺度數據可分解出較多的顯著模態,較詳細的不同時間長度的周期(既有年代際、多年周期,也有年際和年周期,更有半年周期、季節周期和1.8個月周期);與之對比明顯的是基于月尺度數據分解的周期主要是年周期,8 d尺度數據可以發現月尺度數據發現不了的周期和規律.研究發現FEEMD在復雜環境下具有較強的數據分解能力.

NE區域處于氣候變化的敏感區,該區域葉綠素a濃度表征的浮游植物和相關環境因子共同處于海洋和大氣動力學過程影響之下,受到多尺度因素和物理過程的影響,使得該區域的葉綠素a濃度和相關環境因子具有顯著的年內、年際、代際特征和趨勢.將不同時間尺度的葉綠素a濃度及環境因子信號正確地歸因,有助于了解不同季節、ENSO、年代際、長期變化背景下各因子的特征、影響及相互聯系.利用FEEMD方法能成功的將多尺度信息從干擾較多的混合信息中分解出來,較適用于復雜環境下的數據分解,能為研究不同尺度各因子之間的驅動關系提供借鑒.