求解高光譜解混的三算子分裂方法

摘要： 對于高光譜解混問題，提出了一個相應(yīng)的三算子單調(diào)包含和求解該問題的一個分裂方法，其中鄰近因子可以自適應(yīng)選取。該方法還可用于求解更一般的帶有線性復(fù)合的三算子單調(diào)包含問題。數(shù)值實(shí)驗(yàn)表明，該算法的性能遠(yuǎn)遠(yuǎn)超過最近提出的鄰近內(nèi)點(diǎn)方法，并且與其變尺度版本相當(dāng)。

關(guān)鍵詞： 鄰近內(nèi)點(diǎn)法； 單調(diào)包含； 分裂方法； 高光譜解混； 自適應(yīng)

中圖分類號： O221.2

文獻(xiàn)標(biāo)志碼： A

文章編號： 1671-6841（2025）02-0085-04

DOI： 10.13705/j.issn.1671-6841.2023172

A Three-operator Splitting Method for Solving Hyperspectral Unmixing

DONG Yunda， ZHANG Yuanyuan， LI Yiyi

（School of Mathematics and Statistics， Zhengzhou University， Zhengzhou 450001， China）

Abstract： We propose A three-operator monotone inclusion for hyperspectral unmixing was proposed and a corresponding splitting algorithm was given in which proximal factors could be selected adaptively. The algorithm could be also used to solve more general three-operator monotone inclusion with one linear composition. Numerical experiments showed that the algorithm outperformed the recently proposed proximal interior point algorithm and was comparable to its variable metric version.

Key words： proximal interior point algorithm; monotone inclusion; splitting method; hyperspectral unmixing; self-adaption

0 引言

高光譜解混（hyperspectral unmixing）就是將傳感器獲取的混合像元的光譜分解為不同純地物的光譜以及估計(jì)其相應(yīng)的比例。文獻(xiàn)［1］提出了一個高效的數(shù)值解法——鄰近內(nèi)點(diǎn)算法（proximal interior point algorithm，PIPA）。該算法的外迭代是內(nèi)點(diǎn)法，而內(nèi)迭代是運(yùn)用向前、向后分裂算法來求解子問題。

本文則將內(nèi)迭代中的子問題轉(zhuǎn)化為三算子單調(diào)包含，并運(yùn)用文獻(xiàn)［2］中的方法來求解。

1 問題描述

設(shè)s和r分別是獲得圖像的光譜帶數(shù)和像素?cái)?shù)，q是材料（端元）的個數(shù)。我們使用式（1）線性模型來關(guān)聯(lián)數(shù)據(jù)、端元和豐度，

Y=SX+ω，（1）

其中：Y=R^s×r表示測得的高光譜立方體；S∈R^s×q是一個已知的數(shù)據(jù)庫，它的每一列都對應(yīng)著預(yù)計(jì)會出現(xiàn)在場景中的一種材料（或端元）的光譜；

X∈R^q×r表示豐度矩陣，它用來描述每種材料（或端元）在每個像素中的比例或豐度；ω代表加性白噪聲。

更為具體的，我們依據(jù)文獻(xiàn)［1］將高光譜解混轉(zhuǎn)化為約束極小化問題，

minX∈R^q×r

12‖Y-SX‖2F+σ∑qi=1‖（WXi）d‖1，

s.t. ∑qi=1Xi，j≤1，j∈｛1，…，r}，

Xi，j≥0，i∈｛1，…，q}，j∈｛1，…，r}，（2）

其中：‖·‖F(xiàn)是Frobenius范數(shù);Xi∈Rr，i∈｛1，…，q}代表豐度矩陣的第i行;

W∈R^r×r表示小波分解算子;‖（·）d‖1是細(xì)節(jié)小波系數(shù)的1-范數(shù);σ≥0是正則化參數(shù)。

引入對數(shù)內(nèi)點(diǎn)化函數(shù)，可將上述問題轉(zhuǎn)化為

minx∈R^qr

12‖y-diag（S，…，S）x‖22+

σ∑qi=1‖（WPix）d‖1-

μ∑qi=1∑rj=1

lnxi，j+

∑rj=1

ln1-∑qi=1xi，j，（3）

鄭 州 大 學(xué) 學(xué) 報(bào) （理 學(xué) 版）第57卷 第2期董云達(dá)，等：求解高光譜解混的三算子分裂方法

其中：μ>0是障礙因子;

diag（S，…，S）表示塊對角矩陣;

y∈R^sr、x∈R^qr分別表示矩陣Y、X拉直后的向量;

Pix=Xi、Pi∈R^r×n是抽取矩陣;σ=0.01。

本文中使用的是真實(shí)的高光譜數(shù)據(jù)集Urban數(shù)據(jù)集，s取162，q取6，r取256×256，小波分解算子選取的是超過2分辨率的正交db4小波分解矩陣。

此外，我們給出信噪比的定義

SNR=20lg（‖‖/‖x-‖），

其中：表示x的真實(shí)值。

2 算法

最近，文獻(xiàn)［3］在無窮維實(shí)Hilbert空間中提出求解帶有一個線性復(fù)合的三算子單調(diào)包含問題的自適應(yīng)分裂算法，該算法是在文獻(xiàn)［2］的基礎(chǔ)之上，給出了鄰近因子自適應(yīng)選取準(zhǔn)則。引入如下記號，

n=qr，Q=In

D，q=（0，…，0，-1，…，-1）T，

D=diag（a1，…，a1），a1=-（1，…，1）∈Rq，

H=diag（S，…，S）。

我們可以得到問題（3）的最優(yōu)性條件為

0∈σ（‖（WPix）d‖1）+HT（Hx-y）+QTB（Qx-q），其中：

Bz=-μdiag［z^-11，…，z^-1p］。

文獻(xiàn)［1］中PIPA在第l步的子問題為

minx∈R^qr

12‖y-diag（S，…，S）x‖22+

σ∑qi=1‖（WPix）d‖1-μl

∑qi=1∑rj=1

lnxi，j+∑rj=1

ln1-∑qi=1xi，j，

其中：μl>0是障礙因子。PIPA采用的是向前向后分裂方法，而在本文中，我們將每個子問題轉(zhuǎn)化為如下三算子單調(diào)包含問題

0∈σ（‖（WPix）d‖1）+HT

（Rx-y）+QTB（Qx-q），

其中：Bz=-μldiag［z^-11，…，z^-1p］。我們參考文獻(xiàn)［3］提出的自適應(yīng)選取鄰近因子的三算子分裂算法運(yùn)用到求解上述子問題中可以表述為：

記m=65 536，選取初始點(diǎn)：x0=ones（6*m，1）/7，

y0=zeros（7*m，1），u0=zeros（7*m，1）， =0.45，α0=1，μ0=1，

θ=1.01。

計(jì)算

β=100×（‖Q‖1‖Q‖∞/4αk+1/4），

然后分別計(jì)算得到k、k、k，k=uk-（yk-Qxk+q）/β，（αkI+0.01‖（·）d‖1）（k）瘙綍αkxk-ak-QTk，ki-μl1ki=

yki+^ki，i=1，…，（q+1）r。

如果xk=k，uk=k則停止，否則計(jì)算

k=αk‖xk-k‖2+

‖yk-k‖2+

〈k-Qk+q，uk-k〉，

φk=‖xk-k‖2+‖yk-k‖2+

‖k-Qk+q‖2，

γk=kφk。

計(jì)算下一迭代點(diǎn)，

（αkI+RTH）（x^k+1）=

αkxk+ak-γk

（xk-k）+HTy，其中y是矩陣Y拉直后的向量（不要與迭代點(diǎn)yk混淆），

y^k+1=yk-γk

（yk-k），

u^k+1=uk-γk

（k-Qk+q），

a^k+1=αk（xk-x^k+1）+ak-γk（xk-k）。

最后更新鄰近因子αk+1

αk+1=（1+τk）αk，若σk≤0.5，

（1-τk）αk，若σk≥2，

αk，其他，

其中：τk=0.1^k+1; σk選取為

σk=‖αk（x^k+1-xk）‖+

‖（y^k+1-yk）‖+‖u^k+1-uk‖‖HTH（x^k+1-xk）‖，

注：該調(diào)比技巧在文獻(xiàn)［3］中給出，該調(diào)比準(zhǔn)則是對文獻(xiàn)［4］中處理兩算子單調(diào)包含問題時(shí)提出的DR算法（相關(guān)的討論參見文獻(xiàn)［5］）相應(yīng)部分準(zhǔn)則的一個延伸。

根據(jù)文獻(xiàn)［2］中的引理1.14、1.15、假設(shè)2.1和定理2.1，可以證明上述算法1的弱收斂性，詳見文獻(xiàn)［3］。

文獻(xiàn)［1］針對該問題提出了一個鄰近內(nèi)點(diǎn)法，該算法的外迭代同樣是使用內(nèi)點(diǎn)法將問題轉(zhuǎn)化為（3），內(nèi)迭代則是采用向前、向后分裂算法。記f，φ為

f（x）=0.01∑qi=1‖（WPix）d‖1，

φμl（x）=12‖y-diag（S，…，S）x‖22-

μl∑qi=1∑rj=1

lnxi，j+∑rj=1

ln1-∑qi=1xi，j，

選擇合適的變尺度矩陣Ak，通過計(jì)算

k（γ）=（I+γA^-1kf）^-1（xk-γA^-1k

Δ

φμl（xk））

得到候選點(diǎn)k（γ）。若γ=γk滿足某個Armijo型步長準(zhǔn)則，則下一迭代點(diǎn)

x^k+1=k（γk）。在實(shí)際操作中，變尺度矩陣Ak一般取I或是函數(shù)φμl在xk處的Hessian陣及其近似。其他相關(guān)參數(shù)選取詳見文獻(xiàn)［1］。

3 數(shù)值實(shí)驗(yàn)

本節(jié)運(yùn)用算法1求解上述高光譜解混問題，并驗(yàn)證其有效性。本實(shí)驗(yàn)是由Matlab R2018b運(yùn)行。

此外，我們還能將問題（3）轉(zhuǎn)化為

0∈diag（S，…，S）T（diag（S，…，S）x-y）+

0.01‖（·）d‖1+QTδC（Qx-q），

其中：C=R^（q+1）r+;

δC表示集合C上的示性函數(shù)。這種情況下我們記作算法2。文獻(xiàn)［1］中當(dāng)變尺度矩陣

Ak恒為單位矩陣I時(shí)，我們記作算法3，當(dāng)Ak取函數(shù)φμl在xk處的Hessian陣時(shí)，即

Ak=diag［STS，…，STS］+μldiag

1x211，…，1x2qr+

μl∑rt=1（eTxt-1）^-2eteTt，

記作算法4，其中：e=（1，…，1）T∈Rq，xt表示矩陣X的列向量，

et∈R^qr表示從第（t-1）q+1行到第tq行元素為1，其余全為0的列向量。數(shù)值結(jié)果如表1。

表1展示的是四種算法分別用來處理高光譜解混問題獲得的關(guān)于表中五種物質(zhì)：草、樹、屋頂、金屬、泥土的信噪比，根據(jù)信噪比的定義我們可知，信噪比值越高說明算法的數(shù)值性能越好。通過觀察表中數(shù)據(jù)，我們可以發(fā)現(xiàn)，算法1的數(shù)值效果遠(yuǎn)遠(yuǎn)超過算法3，并且與算法4的結(jié)果相差無幾，這說明在不引入變尺度矩陣的情況下，我們的分裂算法效果更好。算法2的數(shù)值結(jié)果也優(yōu)于算法3，但不如算法1和算法4，這說明對于求解該類問題來說，內(nèi)點(diǎn)法要優(yōu)于簡單引入示性函數(shù)。

圖1、圖2分別展示了算法1、算法2、算法4關(guān)于屋頂、泥土這兩種物質(zhì)的豐度圖和真實(shí)圖像的對比圖。通過對比不難發(fā)現(xiàn)，算法1和算法4效果優(yōu)于算法2，并且二者效果差別不大。圖3、圖4分別 展示了相同迭代步數(shù)下精度和信噪比的比較，其中精度ε=lg（‖xk-x*‖/‖x*‖）。通過上圖也能發(fā)現(xiàn)，算法1和算法4能達(dá)到更Dd5sqM58DBf/XQqPUgnoZg==高的精度和信噪比。

4 結(jié)論

本文將高光譜解混問題轉(zhuǎn)化為求解三算子單調(diào)

包含。不同于引入示性函數(shù)，本文第二小節(jié)中的處理方式更容易發(fā)揮三算子分裂算法的性能，這是本文的主要貢獻(xiàn)之一。綜上所述，在不引入變尺度矩陣的情況下，算法1的數(shù)值效果遠(yuǎn)超算法3，并且能達(dá)到和引入變尺度矩陣Ak的算法4近乎相同的效果，但是運(yùn)行時(shí)間相對久一些，因此后續(xù)我們可以考慮將變尺度矩陣Ak引入到我們的算法1之中，期待可以獲得更好的數(shù)值結(jié)果。此外，參考文獻(xiàn)［3］可知，算法1的應(yīng)用范圍更為廣泛，參數(shù)選取更加簡單實(shí)用。

參考文獻(xiàn)：

［1］ CHOUZENOUX E， CORBINEAU M C， PESQUET J C. A proximal interior point algorithm with applications to image processing［J］.Journal of mathematical imaging and vision， 2020， 62（6/7）： 919-940.

［2］ DONG Y D. A new splitting method for systems of monotone inclusions in Hilbert spaces［J］. Mathematics and computers in simulation， 2023， 203： 518-537.

［3］ 張園園. 求解三算子單調(diào)包含的自適應(yīng)分裂方法及其應(yīng)用［D］. 鄭州：鄭州大學(xué)， 2023.

ZHANG Y Y. Adaptive splitting method for solving three-operator monotone inclusions with its applications［D］. Zhengzhou：Zhengzhou University， 2023.

［4］ DONG Y D， FISCHER A. A family of operator splitting methods revisited［J］. Nonlinear analysis： theory， methods & applications， 2010， 72（11）： 4307-4315.

［5］ 康倍倍， 董云達(dá)， 王亞麗. 關(guān)于凸極小化的Douglas-Rachford分裂方法的一個注［J］. 鄭州大學(xué)學(xué)報(bào)（工學(xué)版）， 2017， 38（4）： 94-96.

KANG B B， DONG Y D， WANG Y L. A note on Douglas-rachford splitting method for convex minimization［J］. Journal of Zhengzhou university （engineering science）， 2017， 38（4）： 94-96.