一種高精度并行主偏度分析算法及其在遙感圖像中的應(yīng)用

王大虎 劉 暢* 王 健 姚 鍇 張 振

①(中國科學(xué)院空天信息創(chuàng)新研究院 北京 100190)

②(中國科學(xué)院大學(xué)電子電氣與通信工程學(xué)院 北京 100094)

③(陸軍航空兵學(xué)院 北京 101123)

1 引言

隨著遙感技術(shù)的發(fā)展，人們可以獲得越來越多的高維大規(guī)模遙感圖像數(shù)據(jù)。在過去的幾十年中，許多學(xué)者使用遙感圖像數(shù)據(jù)在光譜融合[1]、目標(biāo)分類[2]、目標(biāo)檢測[3]、圖像去噪[4]、圖像壓縮[5]等研究領(lǐng)域做了大量研究。高維大規(guī)模數(shù)據(jù)所蘊(yùn)含的巨大的信息為后續(xù)遙感圖像處理帶來巨大挑戰(zhàn)，過高的數(shù)據(jù)維度不利于后續(xù)的存儲、傳輸以及運(yùn)算，而特征提取成為解決這一問題的關(guān)鍵技術(shù)。

近幾年研究人員提出了多種用于特征提取的算法，典型方法包括主成分分析方法(Principal Component Analysis, PCA)[6]、奇異值分解(Singu lar Value Decom position, SVD) 和獨(dú)立成分分析方法(Independen t Com ponen t Analysis, ICA)[7]。PCA是基于數(shù)據(jù)的2階相關(guān)統(tǒng)計特征找到一組正交向量集，這個向量集從最小二乘意義上對原始信號進(jìn)行表示。SVD則是在最大能量的意義上選擇較大的變量對原矩陣進(jìn)行分解，以此來提取數(shù)據(jù)矩陣的主要特征。上述方法主要針對數(shù)據(jù)的2階統(tǒng)計特征，并假設(shè)數(shù)據(jù)滿足高斯分布。然而，許多真實(shí)數(shù)據(jù)集的分布通常不滿足高斯分布。在這種情況下，基于高階統(tǒng)計量的ICA方法近年來受到了越來越多的關(guān)注。ICA是對數(shù)據(jù)進(jìn)行高階統(tǒng)計分析以尋求非高斯量最大化來分離各個獨(dú)立的信號源。FastICA是ICA的一種改進(jìn)算法，它可以選擇偏度、負(fù)熵、峭度等指標(biāo)作為非高斯度量[8]。雖然FastICA算法能達(dá)到3次收斂速度，優(yōu)于其他大多數(shù)常用的ICA算法[9]，但是其每次迭代都需要涉及所有的像素來尋找最優(yōu)的投影方向，這導(dǎo)致其在處理高維數(shù)據(jù)時耗時大。

為了解決這一問題，Geng等人[10]提出了一種主偏度分析方法(Principal Skew ness Analysis,PSA)，該方法的求解可以轉(zhuǎn)化為計算數(shù)據(jù)3階統(tǒng)計張量特征對，因此PSA可以看作PCA的3階泛化。當(dāng)選擇偏度作為非高斯指標(biāo)時，可等效為FastICA。在此基礎(chǔ)上，文獻(xiàn)[11,12]進(jìn)一步提出了動量PSA(Momentum Principal Skewness Analysis, MPSA)和主峰度分析方法(Principal Kurtosis Analysis,PKA)。超對稱張量并不具有天然的正交性[13,14]，然而PSA和MPSA算法在求解過程中均采用了正交補(bǔ)策略，因此它們得到的解(特征對)是相互正交的。這導(dǎo)致除了第1個特征對外，得到的其他特征對都偏離了協(xié)偏度張量的實(shí)際特征對。針對這個問題，Geng等人[15]又提出了一種新的非正交的PSA(Nonorthogonal Principal Skew ness Analysis,NPSA)算法。NPSA算法通過引入克羅內(nèi)克積提出了一種新的搜索策略，該策略可以在更大的空間內(nèi)搜索解，因此可以獲得比PSA算法更精確的解，但是NPSA算法得到的仍然是近似解。

針對上述局限，本文提出了一種新的高精度并行主偏度分析(Parallel Principle of Skewness Analysis, PPSA)算法，PPSA算法把協(xié)偏度張量切片的特征向量作為初始值進(jìn)行迭代求解。與現(xiàn)有PSA算法施加約束再依次求解各特征對不同，PPSA算法不再施加約束，而是將各個特征對的初始值并行地、一次性地全部算出，有效地避免了施加約束帶來的誤差。相比以前只能得到近似解，PPSA算法可以準(zhǔn)確得到實(shí)際解。

2 PSA算法

3 PPSA算法

3.1 現(xiàn)有PSA算法局限

現(xiàn)有的 PSA 算法主要采用了串行方式依次求解各個特征對，為了避免收斂到同一特征對而施加約束，但是這也將導(dǎo)致除了第1個特征對外，其他特征對都偏離了實(shí)際特征對[15]。并且現(xiàn)有PSA算法初始化均采用隨機(jī)生成初始值，這意味著給定不同的初始化，可能收斂到不同的解(特征對)[18]。在這里舉一個簡單的例子來更直觀地說明這兩種現(xiàn)象?？紤]一個2階超對稱張量，其兩個切片分別為

這個簡單的例子直觀地說明了現(xiàn)有的PSA算法并不能準(zhǔn)確地得到協(xié)偏度張量的實(shí)際特征對，還證明了現(xiàn)有的PSA算法隨機(jī)生成初始值，導(dǎo)致每次運(yùn)行都會收斂到不同的特征對。并且上述問題會隨著張量維數(shù)的增加表現(xiàn)得更加明顯，在3階2維超對稱張量求解中，可能會收斂到兩個不同的特征對，在3階維超對稱張量求解中，因?yàn)槌跏贾档牟煌赡苁諗康絺€不同的特征對。因此如何設(shè)計出計算求解精確、穩(wěn)定的方法至關(guān)重要。

3.2 PPSA算法

針對上述問題，本文在PSA算法基礎(chǔ)上進(jìn)行改進(jìn)，提出一種高精度并行PSA算法。PPSA算法不再施加約束，而是將各個特征對的初始值并行地全部算出，有效地避免了施加約束和串行計算帶來的誤差。

求解協(xié)偏度張量的實(shí)際特征對可以理解為求解張量的局部極大值解[18]，只要保證初始值落在協(xié)偏度張量各個收斂區(qū)域，然后再使用不動點(diǎn)迭代法進(jìn)行求解，就能保證精確地求解到所有極大值解。因此求解協(xié)偏度張量的特征對問題可轉(zhuǎn)化為初始值選取的問題。常見的初始化方法有隨機(jī)初始化，基于正交約束的隨機(jī)初始化。隨機(jī)初始化的方法雖然可以精確求解，但是它們受初始值影響很大，并且是不可重復(fù)的。因此設(shè)計出一種初始值可重復(fù)，可以穩(wěn)定的精確的求解出協(xié)偏度張量的特征對的方法是很重要的。

圖1 不同初始值，PSA,MPSA算法求解結(jié)果

圖2 不同初始值，NPSA算法求解結(jié)果

PPSA算法充分考慮數(shù)據(jù)結(jié)構(gòu)，選用協(xié)偏度張量的全部切片的特征向量作為迭代的初始值。首先，此方法是可重復(fù)的。其次，此方法是并行的。最后，對比隨機(jī)初始化方法，此方法可以獲得穩(wěn)定精確的解的原因是我們相信它將提供在最優(yōu)解附近的初始猜測。

假設(shè)S ∈R[3,n]，如果λ是協(xié)偏度張量S的一個特征值，設(shè)u是S關(guān)于λ的單位特征向量，即u=[a1,a2,...,a n]T且uTu=1。求 解S×1u×3u=λu，其等價求解

由式(7)可知，求解協(xié)偏度張量的局部極值解等價多個主成分分析聯(lián)合求極大極小值。把協(xié)偏度張量的切片的特征向量作為迭代初始值可以理解為從單個主成分的極值點(diǎn)去逼近聯(lián)合主成分的極值點(diǎn)。因此對比隨機(jī)選取初始值，此初始化方法生成的初始值距離實(shí)際特征對更近。為了驗(yàn)證這一結(jié)論，本文在4.1節(jié)做了蒙特卡羅仿真實(shí)驗(yàn)。

為方便起見，總結(jié)PPSA算法如算法1所示。

需要注意的是，PPSA算法中步驟7中的循環(huán)終止條件包含最小允許誤差ε和最大循環(huán)次數(shù)K。本文ε設(shè)置為0.0001，K設(shè)置為1000。U是最終的非正交主偏度變換矩陣。

算法1 PPSA算法

同樣對上述3.1節(jié)的例子進(jìn)行計算，用PPSA算法重復(fù)上述操作，求解得到的結(jié)果為=[0.8812,-0.4727]T,=[0.3756,0.9268]T。

圖3 PPSA算法求解結(jié)果

4實(shí)驗(yàn)與復(fù)雜度分析

本文將PPSA算法應(yīng)用于盲圖像分離、SAR圖像去噪以及高光譜特征提取，并與幾種經(jīng)典方法開展了對比試驗(yàn)。所有的算法都是在1臺CPU AMD Ryzen 7 5800H,16 GB RAM,@3.20 GHz的筆記本電腦上完成的，所有程序在MATLAB R2021a上編程和實(shí)現(xiàn)。

4.1 蒙特卡羅仿真實(shí)驗(yàn)

為了檢驗(yàn)本文所提出的初始化算法的效果，本文采用蒙特卡羅方法進(jìn)行驗(yàn)證。實(shí)驗(yàn)采用有n個極大值特征對的3階n維超對稱張量進(jìn)行求解，n分別取2～9。采用3種不同的初始化方法進(jìn)行迭代求解，進(jìn)行10萬次蒙特卡羅仿真后，對比求解正確的結(jié)果，其結(jié)果如表1所示。

從實(shí)驗(yàn)結(jié)果上看，基于協(xié)偏度張量的單個切片的特征向量作為初始值的方法好于基于正交約束的隨機(jī)初始化方法以及隨機(jī)初始值的初始化方法。但是隨著協(xié)偏度階數(shù)的增加，此初始化得到全部解的概率會有所下降。為了改善這種情況，可以適當(dāng)增加初始值的數(shù)量，選用全部切片的特征向量作為迭代初始值重復(fù)上述實(shí)驗(yàn)，其結(jié)果如表2所示。

從實(shí)驗(yàn)結(jié)果上看，隨著初值點(diǎn)數(shù)的增加，各種初始化方法得到所有特征對的概率都有明顯提升。尤其是基于全部切片的特征向量作為初始值的方法，其在所有的實(shí)驗(yàn)中準(zhǔn)確率均達(dá)到了99.9%。因此，對比隨機(jī)初始化方法，選用協(xié)偏度張量的全部切片的特征向量作為迭代的初始值距離實(shí)際特征對更近。

4.2 盲圖像分離實(shí)驗(yàn)

本文首先將PPSA算法應(yīng)用于盲圖像分離(B lind Im age Separation,BIS)問題上，BIS的目的是估計混合矩陣，記為B(或其逆矩陣，即分解矩陣，記為U)。本文將PPSA算法與Fast ICA[8],PSA[10],MPSA[11],NPSA[15]和MSDP[18]等算法進(jìn)行比較，評價各算法的分離性能。在本實(shí)驗(yàn)中，本文選取n幅大小為256像素×256像素的灰度圖像作為源圖像，其中n分別取值2～6。由于已知真實(shí)的源圖像，因此可以將其與各算法獲得的結(jié)果進(jìn)行比較。由于篇幅的限制，本文沒有顯示所有實(shí)驗(yàn)結(jié)果，只選擇n=3的情況進(jìn)行展示，結(jié)果如圖4所示。為了確保結(jié)論的可靠性，本文還進(jìn)行了其他兩種組合實(shí)驗(yàn)，在每種組合中本文隨機(jī)選擇3種不同的源圖像進(jìn)行實(shí)驗(yàn)。

從實(shí)驗(yàn)結(jié)果來看，上述6種算法均可以從混合圖像中分離出源圖像，為了定量評估上述6種算法的性能，本文使用5個評估指標(biāo)對上述6種算法獲得的分離結(jié)果進(jìn)行評價。5個評價指標(biāo)分別為符號間干擾(InterSymbol-Interference,ISI)[15]、總均方誤差(Total M ean Square Error,TMSE)[15]、相關(guān)系數(shù)指數(shù)(Correlation Coefficient,記為ρ)[15]、峰值信噪比(Peak Signal to Noise Ratio,PSNR)[15]和運(yùn)行時間(Time,T)。值得注意的是除PPSA算法與MSDP算法外，其余4種算法獲得的結(jié)果均具有隨機(jī)性[18]，因此本文取10次運(yùn)行的平均值作為最終結(jié)果。

表3詳細(xì)地列出了3種不同組合下6種算法的符號間干擾、總均方誤差、相關(guān)系數(shù)、峰值信噪比和運(yùn)行時間。對于ISI和TM SE指標(biāo)，指數(shù)值越小，代表該方法的性能越好；ρ和PSNR越大，代表該方法的性能更好。通過對表3進(jìn)行分析可以得到如下結(jié)論：

表3 FastICA,PSA, MPSA,NPSA,MSDP和PPSA算法評估結(jié)果

(1)MSDP算法可以準(zhǔn)確得到張量的特征對[18]，通過觀察MSDP算法與PPSA算法的評價指標(biāo)，可以發(fā)現(xiàn)PPSA算法與MSDP算法獲得的結(jié)果一致，這說明PPSA算法也可以準(zhǔn)確得到張量的實(shí)際特征對，驗(yàn)證了本文所提方法的有效性。同時，PPSA算法的運(yùn)行時間遠(yuǎn)小于MSDP算法的運(yùn)行時間，這說明PPSA算法在精確求解的同時還可以做到快速求解。

(2)PPSA算法獲得的ISI與TMSE指標(biāo)始終小于其余4種方法的，ISI與TMSE指標(biāo)越小，圖像分離效果越好，這說明PPSA算法在BIS應(yīng)用中可以獲得更好的整體效果。

(3)PPSA算法獲得的相關(guān)系數(shù)與峰值信噪比也在這6種算法中處于較高水平，相關(guān)系數(shù)值和峰值信噪比越大，分離圖像質(zhì)量越好，這說明PPSA算法性能穩(wěn)定。

最后，結(jié)合上述5個指標(biāo)的比較結(jié)果，說明PPSA算法在BIS應(yīng)用中具有更準(zhǔn)確和更魯棒的性能。

4.3 SAR圖像去噪實(shí)驗(yàn)

本文將PPSA算法用于SAR圖像去噪問題上，實(shí)驗(yàn)使用的是高分三號拍攝的一幅美國舊金山圖像，圖像大小為448像素×364像素，極化方式為HV極化。對源圖像施加噪聲來模擬低信噪比圖像，通常認(rèn)為SAR圖像的斑點(diǎn)噪聲服從伽馬分布或者高斯分布，因此分別在源圖上施加均值為0，方差為100的高斯噪聲(記為Egs=0，=100)和均值為2，方差為80的伽馬噪聲(記為Egm=2，=80)來仿真低信噪比的混合圖像。將原始圖像按逐行的方式轉(zhuǎn)化為相應(yīng)的3個1維信號，產(chǎn)生隨機(jī)混合矩陣rand(3,3)。1維信號與隨機(jī)混合矩陣混合后得到混合圖像，再使用FastICA,PSA,NPSA,PPSA、均值濾波和中值濾波等算法對混合圖像進(jìn)行去噪處理，結(jié)果如圖5所示。

因?yàn)镕astICA,PSA和NPSA的處理結(jié)果受隨機(jī)初值影響，所以本文選擇它們10次實(shí)驗(yàn)的最佳的處理結(jié)果進(jìn)行對比以及展示。同樣，本文也選擇出均值濾波以及中值濾波最佳的結(jié)果進(jìn)行對比展示，均值濾波與中值濾波均選擇5×5大小的濾波器。

通過觀察圖5，發(fā)現(xiàn) PSA算法、NPSA算法、PPSA算法的去噪效果優(yōu)于FastICA算法、中值濾波和均值濾波。為了客觀評價分離濾波結(jié)果，本文分別計算圖像去噪后的峰值信噪比(PSNR)[15]和平均絕對誤差(MAE)。對像素為M×N的圖像，其MAE定義為

式中，h(x,y)和g(x,y)分別表示原始圖像與待評價圖像的第 個像素的灰度值。為了保證結(jié)論的可靠性，本文還選擇不同參數(shù)的伽馬噪聲與高斯噪聲進(jìn)行了2次實(shí)驗(yàn)，計算結(jié)果如表4所示。

表4 不同算法的平均峰值信噪比和平均絕對誤差(計算10次運(yùn)行的平均結(jié)果)

MAE越小，算法性能越好。由表4可看出，對于相同強(qiáng)度的高斯噪聲和伽馬噪聲,均值濾波法和中值濾波法的去噪效果相對較差，而FastICA方法相對較好，PSA方法在所有組合中均獲得了不錯的效果，尤其是PPSA方法，能顯著地去除噪聲，提高復(fù)原圖像的質(zhì)量。結(jié)合上述兩個指標(biāo)的比較結(jié)果，說明PPSA算法在SAR圖像去噪處理中具有更強(qiáng)的去噪能力。

4.4 高光譜特征提取實(shí)驗(yàn)

本節(jié)使用高光譜數(shù)據(jù)集Cuprite來評價PPSA算法的性能。該數(shù)據(jù)集是由機(jī)載可見紅外成像光譜儀于1997年6月在內(nèi)華達(dá)州銅礦礦區(qū)獲得的。將數(shù)據(jù)嵌入到可視化圖像環(huán)境軟件中，圖像包含50個波段，每個波段有350像素×400像素，偽彩色圖像如圖6所示。

圖6 高光譜Cup rite數(shù)據(jù)彩色圖像

圖7 4種算法運(yùn)行時間對比圖

本文直接選擇數(shù)據(jù)集Cup rite全部的50個波段進(jìn)行特征提取。當(dāng)選用全部50個波段進(jìn)行特征提取時，MSDP和FastICA算法的時間消耗是巨大的[18]，因此本文使用PSA,MPSA,NPSA和PPSA算法進(jìn)行特征提取。與4.2節(jié)的情況不同，在實(shí)際數(shù)據(jù)處理中，通常不知道獨(dú)立分量的數(shù)量，所以PSA,MPSA和NPSA算法的獨(dú)立分量的數(shù)量總是手動設(shè)置為頻帶的數(shù)量。因此用上述3種算法進(jìn)行求解，它們最終全都獲得50個獨(dú)立分量。相比之下，PPSA算法是參數(shù)自適應(yīng)的，它不需要事先設(shè)置獨(dú)立分量的數(shù)量，它可以自動確定獨(dú)立分量的數(shù)目，用PPSA算法進(jìn)行求解，最終獲得了6個獨(dú)立分量。由于不知道真實(shí)數(shù)據(jù)的參考獨(dú)立分量，因此4.2節(jié)仿真實(shí)驗(yàn)部分使用的評價指標(biāo)不可用。因?yàn)槠南拗疲疚牟灰灰伙@示PSA,MPSA和NPSA這3種算法獲得的結(jié)果。本文僅顯示PPSA算法獲得的結(jié)果，如圖8所示。

圖8 PPSA算法對高光譜圖像Cup rite的偏度比較

其余3種算法雖然均獲得了50個獨(dú)立分量，但是所提取到50個特征對中大多是沒有信息的冗余結(jié)果，并且提取到的特征對除了第1個外其余都為近似解。把PSA,MPSA和NPSA這3種算法得到的50個獨(dú)立分量分別作為初始值進(jìn)行固定點(diǎn)迭代求解，最終全都收斂為6個獨(dú)立分量，這6個獨(dú)立分量與PPSA算法得到的結(jié)果完全相同，說明PPSA算法在對高維數(shù)據(jù)進(jìn)行特征提取時具有更好的自適應(yīng)性和準(zhǔn)確性。

此外，本文還對上述4種算法的運(yùn)行速度進(jìn)行比較。PPSA算法可以并行地求解各特征對，因此PPSA算法可以使用Parfor并行加速計算。本文將使用Parfor并行加速計算的PPSA算法命名為FastPPSA。

圖7繪制了不同數(shù)量波段下4種算法的運(yùn)行時間曲線，根據(jù)Cuprite數(shù)據(jù)各波段的偏度值，從小到大排列依次選取n個 波段進(jìn)行特征提取，n的范圍為2～50。

通過觀察圖7可以看出在處理相同數(shù)量的數(shù)據(jù)樣本時(樣本數(shù)大于20時)，MPSA算法運(yùn)行時間小于FastPPSA算法運(yùn)行時間小于PSA算法運(yùn)行時間小于NPSA算法運(yùn)行時間小于PPSA算法運(yùn)行時間。同時，還可以看出在處理高階數(shù)據(jù)時FastPPSA算法運(yùn)行時間遠(yuǎn)小于PPSA算法運(yùn)行時間。然而在處理低階數(shù)據(jù)時FastPPSA運(yùn)行時間要略大于PPSA算法的運(yùn)行時間。這是因?yàn)楫?dāng)計算機(jī)處理的數(shù)據(jù)量較小時，Parfor并行加速在初始化過程中需要消耗大量的時間，這部分時間在總的運(yùn)行時間里占了很大的一部分。因此，當(dāng)所要處理圖像規(guī)模較小時，利用Parfor并行運(yùn)算的加速效果并不明顯。然而，隨著張量維數(shù)的逐漸增大，耗費(fèi)在并行計算額外開銷上的時間占算法執(zhí)行總時間的比例越來越小，所以并行加速效果明顯。FastPPSA采用Parfor并行計算導(dǎo)致其在處理高維數(shù)據(jù)時運(yùn)算速度很快，處理50個樣本數(shù)據(jù)，其運(yùn)行速度為3.1028 s,NPSA運(yùn)行速度為17.2797 s,PPSA算法運(yùn)行時間為28.1447 s。

結(jié)合上述實(shí)驗(yàn)結(jié)果與分析可知，本文所提出的PPSA算法不僅能夠準(zhǔn)確地對高光譜數(shù)據(jù)進(jìn)行特征提取，而且還可以做到快速提取。這說明PPSA算法在高光譜數(shù)據(jù)特征提取上具有更好的準(zhǔn)確性、自適應(yīng)性和魯棒性。

4.5 復(fù)雜度分析

本節(jié)從理論上比較PSA算法和PPSA算法的計算復(fù)雜度。由于這兩種方法基本流程都一樣，只是在求解張量特征對時主循環(huán)的次數(shù)不同，因此只需考慮主循環(huán)求解的計算復(fù)雜性。對于L×M大小的數(shù)據(jù)，PSA算法只需要循環(huán)L次，假設(shè)每次循環(huán)需要迭代K次，其中K是平均迭代次數(shù)，則PSA算法需要進(jìn)行K L次操作。同樣對于L×M大小的數(shù)據(jù)，PPSA算法則需要循環(huán)L2次，同樣也假設(shè)每次循環(huán)需要迭代K次，因此PPSA算法最終需要K L2次操作。雖然本文PPSA算法的計算復(fù)雜度略高于PSA算法，但PPSA算法相對于PSA算法的主要優(yōu)勢在于PPSA算法是并行的，所以它可以采用并行計算來提升算法運(yùn)行速度。由前面的高光譜運(yùn)行速度對比實(shí)驗(yàn)可知，在處理高維數(shù)據(jù)時，采用并行加速計算后的PPSA算法運(yùn)行速度更快。

5 結(jié)束語

本文針對現(xiàn)有PSA算法不能準(zhǔn)確得到實(shí)際解這一問題，在原有的PSA算法的基礎(chǔ)上進(jìn)行改進(jìn)，提出了一種高精度并行PSA算法。與現(xiàn)有PSA算法串行求解特征對不同，本文所提算法不再施加任何約束，而是選擇張量全部切片的特征向量作為初始值進(jìn)行自由迭代求解，有效地避免了施加約束帶來的誤差。通過將PPSA算法與現(xiàn)有PSA算法進(jìn)行實(shí)驗(yàn)對比，驗(yàn)證了PPSA算法在盲圖像分離，SAR圖像去噪以及高光譜數(shù)據(jù)特征提取上具有更好的準(zhǔn)確性和魯棒性。PSA,NPSA以及PPSA都關(guān)注數(shù)據(jù)集的3階偏度。對于其他數(shù)據(jù)集，偏度可能并不總是最佳選擇。在這種情況下，可以選用其他指標(biāo)，如4階峰度或其他高階統(tǒng)計量替代。