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

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

關鍵詞： 鄰近內點法； 單調包含； 分裂方法； 高光譜解混； 自適應

中圖分類號： O221.2

文獻標志碼： 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）就是將傳感器獲取的混合像元的光譜分解為不同純地物的光譜以及估計其相應的比例。文獻［1］提出了一個高效的數值解法——鄰近內點算法（proximal interior point algorithm，PIPA）。該算法的外迭代是內點法，而內迭代是運用向前、向后分裂算法來求解子問題。

本文則將內迭代中的子問題轉化為三算子單調包含，并運用文獻［2］中的方法來求解。

1 問題描述

設s和r分別是獲得圖像的光譜帶數和像素數，q是材料（端元）的個數。我們使用式（1）線性模型來關聯數據、端元和豐度，

Y=SX+ω，（1）

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

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

更為具體的，我們依據文獻［1］將高光譜解混轉化為約束極小化問題，

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是Frobenius范數;Xi∈Rr，i∈｛1，…，q}代表豐度矩陣的第i行;

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

引入對數內點化函數，可將上述問題轉化為

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）

鄭 州 大 學 學 報 （理 學 版）第57卷 第2期董云達，等：求解高光譜解混的三算子分裂方法

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

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

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

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

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

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

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

其中：表示x的真實值。

2 算法

最近，文獻［3］在無窮維實Hilbert空間中提出求解帶有一個線性復合的三算子單調包含問題的自適應分裂算法，該算法是在文獻［2］的基礎之上，給出了鄰近因子自適應選取準則。引入如下記號，

n=qr，Q=In

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

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

H=diag（S，…，S）。

我們可以得到問題（3）的最優性條件為

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

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

文獻［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采用的是向前向后分裂方法，而在本文中，我們將每個子問題轉化為如下三算子單調包含問題

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

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

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

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

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

θ=1.01。

計算

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

然后分別計算得到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則停止，否則計算

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。

計算下一迭代點，

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

αkxk+ak-γk

（xk-k）+HTy，其中y是矩陣Y拉直后的向量（不要與迭代點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）‖，

注：該調比技巧在文獻［3］中給出，該調比準則是對文獻［4］中處理兩算子單調包含問題時提出的DR算法（相關的討論參見文獻［5］）相應部分準則的一個延伸。

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

文獻［1］針對該問題提出了一個鄰近內點法，該算法的外迭代同樣是使用內點法將問題轉化為（3），內迭代則是采用向前、向后分裂算法。記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，通過計算

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

Δ

φμl（xk））

得到候選點k（γ）。若γ=γk滿足某個Armijo型步長準則，則下一迭代點

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

3 數值實驗

本節運用算法1求解上述高光譜解混問題，并驗證其有效性。本實驗是由Matlab R2018b運行。

此外，我們還能將問題（3）轉化為

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上的示性函數。這種情況下我們記作算法2。文獻［1］中當變尺度矩陣

Ak恒為單位矩陣I時，我們記作算法3，當Ak取函數φμl在xk處的Hessian陣時，即

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的列向量。數值結果如表1。

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

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

4 結論

本文將高光譜解混問題轉化為求解三算子單調

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

參考文獻：

［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］ 張園園. 求解三算子單調包含的自適應分裂方法及其應用［D］. 鄭州：鄭州大學， 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］ 康倍倍， 董云達， 王亞麗. 關于凸極小化的Douglas-Rachford分裂方法的一個注［J］. 鄭州大學學報（工學版）， 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.