CN106778530A - 一种基于双线性混合模型的高光谱图像非线性解混方法 - Google Patents

Abstract

本发明属于遥感图像处理技术领域，具体为一种基于双线性混合模型的高光谱图像非线性解混方法。本发明根据双线性混合模型的几何特性，通过将数据中的非线性混合成分表示为一个融合了共同非线性效应端点的线性贡献，使复杂的非线性解混转化为简单的线性解混问题，进而结合线性解混算法迭代估计正确的丰度。本发明从高光谱观测像元的混合模型出发，结合其几何与物理意义，不仅能有效地弥补线性解混的不足，而且能较好地克服共线性效应带来的不利影响。具有对噪声及端元数目较好的鲁棒性，可以作为一种解决高光谱遥感图像非线性解混的有效手段。在基于高光谱遥感图像的高精度解混以及地面目标的检测和识别方面具有重要的应用价值。

Description

一种基于双线性混合模型的高光谱图像非线性解混方法

技术领域

本发明属于遥感图像处理技术领域，具体涉及一种高光谱图像非线性解混方法。

背景技术

遥感技术是本世纪六十年代发展起来的新兴综合技术，与空间、电子光学、计算机、地理学等科学技术紧密相关，是研究地球资源环境的最有力的技术手段之一。高光谱遥感是将成像技术与光谱技术相结合的多维信息获取技术。其图像具有光谱分辨率高、图谱合一的特点，为地物信息的提取和分析提供了极其丰富的信息。然而由于空间分辨率较低以及地物分布的复杂多样性，图像中的像元大多是由数种不同的纯物质光谱按照一定方式构成的混合像元。混合像元的存在使得传统的像元级高光谱图像的分类、目标检测等应用受到限制[1]。为解决该问题，使分析和应用在亚像元级进行，就需要实现光谱解混，将这些像元分解为其组成成分的物质光谱(端元)及在像元对应区域中的比例(丰度)[2]。光谱解混可以大大提高高光谱遥感对亚像元级目标的探测能力，已被广泛应用于地质科学、环境科学、精细农业以及军事等领域中。

以往经常采用线性混合模型(Linear Mixture Model,LMM)来描述混合像元现象，认为像元是由端元依各自丰度线性组合而成。虽然LMM模型简单且具有一定物理意义，但是因为难以准确描述如沙地、矿物混合区域的紧密混合以及植被覆盖与城市地区的多层次混合等非线性混效应，会给结果带来较大误差。因此，为满足更高的应用要求，各种非线性混合模型及算法被提出来用于弥补线性解混方法的不足[3],[4]。其中，包括常用于紧密混合场景的Hapke模型[5]和仅考虑两两端元间二次散射效应的双线性混合模型(BilinearMixture Model,BMM)。BMM主要由线性混合与非线性混合两部分构成，根据后者的不同可分为Fan模型(Fan Model,FM)[6]，广义双线性模型(Generalized Bilinear Model,GBM)[7]，多项式后验非线性模型(Polynomial Post-Nonlinear Model,PPNM)[8]等。由于BMM具有较好的物理意义且模型相对简单而常作为目前非线性解混算法研究的基础。端元已知时，主要的相关丰度估计算法包括一阶泰勒展开的FCLS(Fully Constrained Least Squares)，梯度下降方法(Gradient Descent Algorithm,GDA)和Bayesian方法等[3],[4]。但是它们存在运算过于复杂且常会陷入局部极小的问题，而端元与虚拟端元间强相关性产生的共线性效应也会导 致结果过拟合并对噪声更为敏感[9]。因此，BMM解混算法需要得到进一步改进，以在提高算法的效率同时获取更高精度的非线性解混结果。

下面介绍与本发明相关的一些概念

双线性混合模型

对于高光谱数据X＝[x₁,x₂,…,x_m]∈R^n×m，其每列x_j∈R^n×1(j＝1,2,...m)都对应一个具有n个波段的像元向量，共存在m个这样的像元。以A＝[a₁,a₂,...,a_r]∈R^n×r表示端元矩阵(r是端元数)，s_ij是端元a_i在像元x_j中的丰度，ε_j是模型误差。LMM假设下的像元是端元按丰度的线性组合，为满足物理意义，丰度需满足非负与“和为一”的约束条件[2]：

BMM则在LMM的基础上增加了两两端元间的二次散射混合项，对于FM，GBM和PPNM来说，像元分别表示为如式(13)(14)(15)所示：

其中，a_i.*a_k＝[a_i,1a_k,1,a_i,2a_k,2,...,a_i,na_k,n]^T为端元间的Hadamard乘积，也称为虚拟端元，用于描述光线在两种端元物质间发生的二次散射非线***互作用。在FM中，该二次散射的强度等于对应端元丰度间的乘积，当像元中不存在某端元时也就不存在与该端元相关的二次散射。但是因为FM是个相对严格的模型，GBM便在FM的基础上向两两端元间的交互中引入了一个非线性参数γ_i,k,j∈[0,1]使模型变得更加灵活，且当对于 ,…,r-1,k＝i+…1,有rγ,_i,k,j＝0时GBM就是LMM，同样若γ_i,k,j＝1则GBM就等价为FM。PPNM相比前两种模型多考虑了端元与自身的二次散射，并用一个参数b_j调控像元中的非线性影响。

发明内容

本发明的目的在于提出一种精度高且运算复杂度低的高光谱图像非线性解混方法。

本发明提出的基于双线性混合模型高光谱图像非线性解混算法，具体步骤为：首先通过几何方法确定一个集中所有二次散射效应的非线性端点，使双线性混合数据近似投影为 其线性混合部分，从而转换为简单的线性解混问题。由于无需在优化中考虑随端元数目增加而大幅增加的虚拟端元和部分参数并直接求解丰度，在减弱共线性对求解精度影响的同时降低了算法复杂度。然后根据投影的近似关系，利用传统的线性解混算法迭代地修正投影位置并估计得到最后的丰度。与其他优秀的同类方法相比，本发明对噪声及端元数目具有较好的鲁棒性，不仅提高了丰度估计精度而且减少了运算时间。具体内容介绍如下：

一、双线性混合模型几何特性的分析

高光谱图像中的每个像元x_j都能被视为是n维空间中的一个点，在LMM假设下，由于丰度需满足非负与“和为一”的约束条件，端元集[a₁,a₂,...,a_r]构成了r-1维子空间(r是端元数)中一个包裹所有数据点的最小单形体Simplex_L。例如三个端元就构成了二维平面上的一个三角形。另一方面，从数学意义来说，BMM中的线性混合部分就是LMM，而且丰度同样需满足非负与“和为一”这两个约束条件。这样BMM的非线性混合项可视为是对LMM增加的一个扰动，使原本位于单形体中的数据向外部偏移，在r维子空间中形成类似鼓包状的分布。由于不同非线性系数的作用，各模型的数据分布存在差异。

由于r维子空间需由其中r个线性无关向量张成，可认为在该子空间中的一个位于r-1维子单形Simplex_L外部的端点与已知的r个端元共同确定了BMM数据点在该空间中的分布。也就是说FM，GBM和PPNM对应的数据点实际上都位于由真实端元和一个额外端点所构成的r维仿射包Affine_BL中，且Simplex_L为Affine_BL中的一个r-1维子单形。任意通过非线性端点与BMM数据点的直线都将与Simplex_L相交，投影点近似为BMM数据点的线性混合部分。

二、非线性端点p的确定

非线性端点p决定着BMM数据的投影位置，对最终解混结果起到关键的影响。实际上，每r-1个端元分别与它们所形成的对应BMM数据点确定了r维子空间中r个r-1维超平面，这些平面分别表征缺少一个端元时的二次散射效应。而点p正是这些超平面的交点并与端元构成最小仿射包，是所有二次散射作用的集中体现。点p对数据的线性贡献就是它们的非线性混合成分。

在确定这些平面时，先根据BMM和已知的端元计算各平面上的非线性混合中点ω₁,...,ω_r-1,ω_r(即各端元丰度相同)，对于FM和GBM如式(5)而PPNM则按式(6)计算：

然后用主成分分析(Principle Component Analysis,PCA)算法将数据降到r维，顶点集{ω₁,a₂,...,a_r},{a₁,ω₂,...,a_r}，…，{a₁,a₂,...,ω_r}便确定该空间中r个r-1维超平面H₁…H_r。以单形体体积为0表示共面，非线性端点p∈R^r×1是适定线性方程组(7)的唯一解：

其中U^T为数据前r个特征向量构成的矩阵，而方程组(7)的解存在且唯一，即这些平面在r维空间中总有唯一的交点p。在确定端点p后，可利用最小二乘法求解式(8)的最优解即可得到像元x_j在中的重心坐标h_j：

同时x_j在中的投影坐标s′_j为：

三、投影位置的修正与丰度估计

s′_j确定的投影位置将与BMM像元的线性混合部分所对应位置间存在一个较小的偏差，但它与像元真实丰度s_j对应的BMM非线性混合部分间存在近似的线性关系。利用该关系可以通过反复修正数据的投影位置，使其不断地向真实的线性混合位置逼近，而同时投影坐标s′_j也将逐步收敛于真实丰度。

求得投影坐标s′_j后，根据各BMM计算其非线性混合部分接着求解优化问题(10)得到使像元到其所在直线距离最短的系数λ_j：

其中，对于FM和GBM有PPNM是 求得式(10)的最小二乘解继而可将投影位置修正为修正后的点y_j将比最初的投影更逼近真实线性混合位置。最后就可用传统的线性丰度估计算法求解新投影y_j的丰度，这里采用的是FCLS[10]算法，并以此反复迭代最后收敛于正确丰度。而在求得丰度后，BMM中的非线性参数用线性规划方法直接求得。

根据上述内容，本发明采用的算法中采用的具体步骤归纳如下：

输入：高光谱数据X∈R^n×m和端元矩阵A∈R^n×r；

输出：丰度矩阵S∈R^r×m；

步骤1确定非线性端点p：

1.1：FM和GBM按式(5)而PPNM按式(6)分别计算各自对应的r个非线性混合中点ω₁,...,ω_r-1,ω_r(即各端元丰度值相等)；

1.2：利用PCA算法将数据降到r维，确定该空间中r个r-1维超平面H₁…H_r；

1.3：求解线性方程组(7)，非线性端点p即为其唯一解。

步骤2修正投影与估计丰度：

2.1：用最小二乘法计算(8)的最优解并通过式(9)得到近似丰度s′_j；

2.2：对于每一个观测像元x_j,j＝1,...,m执行如下操作：

2.2a)初始化迭代次数t＝0和丰度

2.2b)计算对应投影点的BMM非线性混合向量

2.2c)求解最优化问题(10)得系数并计算新投影位置

2.2d)用线性丰度估计算法求解的丰度s^(t+1)；

2.2e)如达到最大迭代次数或满足收敛精度则输出丰度矩阵S＝S^(t+1)，否则回到步骤2.2b)

步骤3：输出结果S＝(s₁,s₂,…,s_m)并用线性规划求得非线性参数γ与b

本发明中：

步骤2.2c)中的系数为最小二乘解：

步骤2.2d)中采用FCLS算法求解s^(t+1)。

步骤2.2e)中的最大迭代次数为200。

步骤2.2e)中，通过判断每次迭代与前一次迭代之间的相对变化是否小于10^-8，来判断是否收敛。

本发明的有益效果在于：其从高光谱观测像元的混合模型出发，结合其几何与物理意义，不仅能有效地弥补线性解混的不足，而且能较好地克服共线性效应带来的不利影响，可以作为一种解决高光谱遥感图像非线性解混的有效手段。在基于高光谱遥感图像的高精度解混以及地面目标的检测和识别方面具有重要的应用价值。

仿真和真实高光谱遥感图像的实验表明，与同类算法相比，本发明方法中采取的算法具有更好的丰度估计精度以及较低的运算复杂度，具有对噪声及端元数目较好的鲁棒性，这对于处理复杂的高光谱数据来说具有重要的实际意义。

附图说明

图1模型数据的几何分布。其中，(a)LMM，(b)FM，(c)GBM，(d)PPNM。

图2方法原理示意图。其中，(a)非线性端点的确定，(b)FM数据投影与投影偏差。

图3投影偏差的修正。

图4 USGS光谱库中5条端元光谱。

图5算法GAEB-FCLS收敛性分析。

图6真实高光谱遥感图像。其中，(a)AVIRIS Cuprite及子区域(b)HYDICEWashingtonDC及子区域。

图7 AVIRIS数据各算法估计的丰度图。其中，从左到右每列为(a)FCLS，(b)KFCLS，(c)Fan-FCLS，(d)PPNM-GDA，(e)GBM-GDA，(f)GAEB-FCLS(FM)，(g)GAEB-FCLS(GBM)，(h)GAEB-FCLS(PPNM)；从上到下每行分别为的蒙脱石，沙漠地表和明矾石三种物质的丰度。

图8 HYDICE数据算法GAEB-FCLS的丰度图。其中，从上到下每行依次是分FM,GBM和PPNM下的估计结果；每列分别对应于5种端元的丰度：(a)水体，(b)屋顶，(c)树木，(d)道路，(e)草地。

具体实施方式

下面，分别用模拟数据和实际遥感图像数据为例说明本发明的具体的实施方式。

本发明中采用的基于双线性混合模型的高光谱遥感图像非线性解混算法用GEAB-FCLS表示。

1、模拟数据实验

在本节中将GAEB-FCLS算法与线性丰度估计算法FCLS[10]，基于高斯核的数据驱动非线性解混方法KFCLS[12](核参数通过交叉验证法在0.01-300间取得)以及FM，GBM和PPNM三个模型对应的传统求解算法：Fan-FCLS[6]，GBM-GDA[7]以及PPNM-GDA[8]进行性能比较。并利用丰度的均方根误差RMSE(Root Mean Square Error)和图像的重构误差RE(Reconstructed Error)两个因子来评价各算法丰度估计的精度：

其中，s_i和x_i是真实丰度和数据而和是估计的丰度和按模型重构的数据，m是像元数，n是波段数，r是端元数。

首先，从美国地质调查局(USGS)光谱库中选取了10种不同地物的光谱作为端元生成模拟数据，图4列出了其中包括植被、沙地等地物在内的五条端元光谱。该光谱库的数据覆盖了0.38μm-2.5μm波长区间内的224个光谱波段，光谱分辨率为10nm。数据的生成方式如下：首先随机生成服从Dirichlet分布[11]的数据丰度，然后根据三种双线性混合模型FM，GBM和PPNM分别构造三类BMM数据。其中GBM数据的非线性参数γ在[0,1]中随机取值，而PPNM数据中参数b取值范围是[-0.3,0.3]。每类数据都设置了三组实验，分别通过改变噪声和端元数目并比较不同像元数目下运行时间来全面地评价算法的性能。每个实验都在相同条件下独立运行10次。

实验1噪声的影响分析：在该实验中分析比较了各算法的噪声鲁棒性。采用像元数为2000，图4中5个端元构成的三类BMM数据，均分别加入20dB，30dB，40dB，50dB和60dB的高斯白噪声并考虑无噪声的情况。表1，表2和表3依次列出了FM，GBM和PPNM三种模型下各算法对于不同信噪比噪声数据的丰度估计精度。总体来看，随着数据信噪比的降低，所有算法对三类数据的解混结果精度也逐渐变差，而本算法相比其它算法在所有情况下丰度估计精度总是最好的。线性算法FCLS的结果都要差于其它四种非线性解混算法，KFCLS算法虽然重构误差都小于FCLS但RMSE较差，表现出过拟合的问题。由于PPNM模型相对FM和GBM具有更好的鲁棒性[8]，因此PPNM-GDA算法在三类数据的结果都优于Fan-FCLS和GBM-GDA两种BMM解混算法。虽然20dB噪声时PPNM-GDA的RE较小，但是丰度的RMSE却差于本文的算法，也表现出一定的过拟合。该实验表明GAEB-FCLS相比其它算法对受不同程度噪声影响的数据具有优越性。

表1 FM不同噪声的数据下的各算法结果比较(×10^-2)

表2 GBM不同噪声的数据下的各算法结果比较(×10^-2)

表3 PPNM不同噪声的数据下的各算法结果比较(×10^-2)

实验2端元数目的影响分析：在该实验中分析比较了各算法在不同端元数目条件下解混结果的精度。数据的像元数为2000，加入50dB的高斯白噪声，端元数目分别取3，5和8。表4，表5和表6比较了FM，GBM和PPNM三种模型下各算法对于不同端元数目数据的解混结果。从丰度的RMSE来看，参与对比的算法结果基本上都随着端元数据数目的增加而降低，但KFCLS结果的RE却表现出了显著的随端元数目增加而减小的趋势，进一步显示出它的过拟合问题。对于三类BMM数据，算法GAEB-FCLS解混结果的RMSE和RE都比其它算法要小，依旧表现出最好的解混性能。

表4 FM不同端元数目的数据下的各算法结果比较(×10^-2)

表5 GBM不同端元数目的数据下的各算法结果比较(×10^-2)

表6 PPNM不同端元数目的数据下的各算法结果比较(×10^-2)

实验3像元数目与算法收敛性分析：本实验首先分析了所提出方法的收敛性，然后以不同像元数目和端元数目测试各算法并比较运行时间。图5描绘了GAEB-FCLS算法对2000个像元的FM，GBM和PPNM三类数据的收敛曲线。可以明显看到丰度的RMSE都随着迭代单调下降，而且算法基本在40次迭代左右收敛。其中，PPNM和FM两种模型下的结果精度几乎收敛于0，另外，虽然由于GBM非线性参数γ的微弱投影扰动使收敛精度稍差于其它两种模型，但依然较快地收敛到了较小值。为了比较各算法的运行时间，在端元数为5和50dB高斯白噪声的条件下，分别取像元数为1000，2000，4000，6000和8000然后比较平均结果，如表7所示。此外，表8中还比较了在2000像元，50dB高斯白噪声时，端元数取3到10的算法运行时间。从结果来看，随着像元数的增多所有算法的运行时间也随之增加，而对比其它两种解混精度较好的非线性解混算法GBM-GDA和PPNM-GDA，GAEB-FCLS的计算时间减少较多。随着端元数目的增加，两种GDA算法尤其是PPNM-GDA的运行时间大量增加，而GAEB-FCLS的时间不仅少于前面两者而且受端元数目影响不大。由此可见，本文提出的GAEB-FCLS算法能避免传统BMM解混算法受端元数目增多导致计算复杂度很大的问题。

表7不同像元数目下算法运行时间比较(s)

表8不同端元数目下算法运行时间比较(s)

2、真实遥感图像实验

本节首先采用了机载可见光/红外成像光谱仪(Airborne Visible/InfraredImaging Spectrometer,AVIRIS)获取于1997年6月19日的美国内华达州Cuprite矿区高光谱图像来分析和评价各个算法的丰度估计结果。该图像大小为512×614(如图6(a)所示)，在0.4-2.5μm的波长区间内有224个波段，光谱分辨率为10nm，空间分辨率为20m。解混前剔除噪声影响较大及水汽吸收波段(1-2,104-113,148-167,221-224)后剩余188个波段。另外，还将HYDICE华盛顿国家广场Washington DC Mall图像(图6(b))作为第二个实验数据，它具有0.4-2.4μm间的210个波段，同样去除部分信噪比较低及水汽吸收波段(1-4,76,87,101-111,136-153,198-210)后剩余191个波段。实验中以VCA算法[11]提取端元后用于丰度估计。通过将解混结果的端元物质丰度图与地物真实分布进行对比，并计算图像的重构误差RE来评价各算法解混性能。

从第一幅AVIRIS图像中截取大小为50×80的子区域作为测试数据，如图6(a)所示。利用HySime算法[13]并根据历史研究和实地调查结果[14][15]可知，AVIRIS的该子区域中主要包括蒙脱石，沙漠地表和明矾石3种端元物质，并且以VCA算法提取端元。图7是AVIRIS图像各算法解混结果对应的丰度图，与其他算法相比，PPNM-GDA和GAEB-FCLS两种算法都能更好地反演得到明矾石接近于其地物真实的分布。表9种对这些算法的图像重构误差进行了比较，GAEB-FCLS在不同BMM假设下解混结果的RE都要小于对应的非线性解混算法。而对于第二个HYDICE数据，则在实验中将其大小为100×100的子区域图像用于算法测试(如图6(b)所示)。利用与AVIRIS数据相同的方法，可知该场景中主要存在水体，屋顶，树木，道路和草地5种地物。从图8可知，算法GAEB-FCLS分别在FM，GBM和PPNM三种模型假设下的丰度反演图，都能正确反映出5种地物的真实分布。从表9来看，具有过拟合问题的KFCLS算法对于HYDICE图像的RE最小，但算法GAEB-FCLS与对应BMM下的非线性解混算法的结果RE相比都相对较小。

表9真实数据的算法重构误差对比(×10^-2)

综上可知，对于模拟和实际高光谱数据来说，本发明提出的算法相对于其他类似算法而言，都具有较好的估计精度和较快的收敛速度，可以快速而精确地实现高光谱图像的解混。

参考文献：

[1]Bioucas-Dias J M,Plaza A,Dobigeon N,et al.Hyperspectral unmixingoverview:Geometrical,statistical,and sparse regression-based approaches[J].IEEE J.Sel.Topics Appl.Earth Observ.Remote Sens.,2012,5(2):354–379.

[2]Keshava N,Mustard J F.Spectral unmixing[J].IEEE SignalProcess.Mag.,2002,19(1):44–57.

[3]Heylen R,Parente M,Gader P.Areview of nonlinear hyperspectralunmixing methods[J].IEEE J.Sel.Topics Appl.Earth Observ.Remote Sens.,2014,7(6):1844–1868.

[4]Dobigeon N,Tourneret J,Richard C,et al.Nonlinear unmixing ofhyperspectral images:models and algorithms[J].IEEE Signal Process.Mag.,2014,31(1):82–94.

[5]Hapke B W.Bidirectional reflectance spectroscopy.I.Theory[J].J.Geophys.Res.,1981,86:3039–3054.

[6]Fan W,Hu B,Miller J,et al.Comparative study between a newnonlinear model and common linear model for analysing laboratory simulatedforest hyperspectral data[J].Int.J.Remote Sens.,2009,30(11):2951–2962.

[7]Halimi A,Altmann Y,Dobigeon N,et al.Nonlinear unmixing ofhyperspectral images using a generalized bilinear model[J].IEEETrans.Geosci.Remote Sens.,2011,49(11):4153–4162.

[8]Altmann Y,Halimi A,Dobigeon N,et al.Supervised nonlinear spectralunmixing using a post-nonlinear mixing model for hyperspectral imagery[J].IEEE Trans.Image Process.,2012,21(6):3017–3025.

[9]Ma L,Chen J,Zhou Y,et al.Two-step constrained nonlinear spectralmixture analysis method for mitigating the collinearity effect[J].IEEETrans.Geosci.Remote Sens.,2016,54(5):2873–2886.

[10]Heinz D C,Chang C.Fully constrained least squares linear spectralmixture analysis method for material quantification in hyperspectral imagery[J].IEEE Trans.Geosci.Remote Sens.,2001,39(3):529–545.

[11]Nascimento J M P,Dias J M B.Vertex component analysis:Afastalgorithm to unmix hyperspectral data[J].IEEE Trans.Geosci.Remote Sens.,2005,43(4):898–910.

[12]Broadwater J,Chellappa R,Banerjee A,et al.Kernel fullyconstrained least squares abundance estimates[C].Barcelona:2007.

[13]Bioucas-Dias J M,Nascimento J M P.Hyperspectral subspaceidentification[J].IEEE Trans.Geosci.Remote Sens.,2008,46(8):2435–2445.

[14]Swayze G,Clark R,Sutley S,et al.Ground-truthing AVIRIS mineralmapping at Cuprite,Nevada[J].Summaries of the Third Annual JPL AirborneGeosciences Workshop,1992,AVIRIS Workshop JPL Publication:47–49.

[15]Bethel J,Lee C,Landgrebe D A.Geometric registration andclassification of hyperspectral airborne pushbroom data[J].InternationalArchives of Photogrammetry and Remote Sensing,2000,33:183–190.。

Claims (5)

1.一种基于双线性混合模型的高光谱图像非线性解混方法，其特征在于：在双线性混合模型的基础上，引入一个非线性端点使高光谱遥感图像的非线性解混问题投影转换为简单的线性解混，然后采用线性解混算法在修正投影位置的同时迭代地求解该问题；具体步骤如下：

(1)确定双线性混合模型的非线性端点

基于不同的光谱混合模型，高光谱遥感图像X＝[x₁,x₂,…,x_m]∈R^n×m中的每个像元x_j∈R^n×1(j＝1,2,...m)可以表示为：

$x_j^{LMM}=\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{a}_i{s}_{ij}+{\mathrm{\ϵ}}_j---\left(1\right)$

$x_j^{FM}=\underset{i=1}{\overset{r}{\Σ}}{a}_i{s}_{i,j}+\underset{i=1}{\overset{r-1}{\Σ}}\underset{k=i+1}{\overset{r}{\Σ}}(a_i.*a_k)s_{i,j}{s}_{k,j}+{\mathrm{\ϵ}}_j---\left(2\right)$

$x_j^{GBM}=\underset{i=1}{\overset{r}{\Σ}}{a}_i{s}_{i,j}+\underset{i=1}{\overset{r-1}{\Σ}}\underset{k=i+1}{\overset{r}{\Σ}}(a_i.*a_k){\mathrm{\γ}}_{i,k,j}{s}_{i,j}{s}_{k,j}+{\mathrm{\ϵ}}_j---\left(3\right)$

$x_j^{PPNM}=\underset{i=1}{\overset{r}{\Σ}}{a}_i{s}_{i,j}+b_j\underset{i=1}{\overset{r}{\Σ}}\underset{k=1}{\overset{r}{\Σ}}(a_i.*a_k)s_{i,j}{s}_{k,j}+{\mathrm{\ϵ}}_j---\left(4\right)$

式(1)到式(4)分别对应于线性混合模型LMM和三种双线性混合模型(BMM)：FM，GBM及PPNM；其中，A＝[a₁,a₂,...,a_r]∈R^n×r表示端元矩阵，r是端元数，s_ij是端元a_i在像元x_j中的丰度，ε_j是噪声误差；丰度皆满足非负及“和为1”约束：

$s_{i,j}\≥0,\underset{i=1}{\overset{r}{\Σ}}{s}_{i,j}=1---\left(5\right)$

非线性参数γ_i,k,j(0≤γ_i,k,j≤1)与b_j用于调控像元中的非线性影响程度；

首先计算双线性混合模型中的r个非线性混合中点，然后对数据降维以确定r个超平面，最后计算这些超平面的交点，即得到非线性端点，具体流程为：

输入高光谱数据X∈R^n×m和端元矩阵A∈R^n×r；

步骤1.1：FM和GBM按式(6)而PPNM按式(7)分别计算各自对应的r个非线性混合中点ω₁,...,ω_r-1,ω_r(即各端元丰度值相等)：

${\mathrm{\ω}}_q=\frac{1}{r-1}\underset{i=1}{\overset{r}{\Σ}}{a}_i+\frac{1}{{(r-1)}^2}\underset{i=1}{\overset{r-1}{\Σ}}\underset{k=i+1}{\overset{r}{\Σ}}(a_i.*a_k),(i\≠q,k\≠q,q=1,...,r)---\left(6\right)$

${\mathrm{\ω}}_q=\frac{1}{r-1}\underset{i=1}{\overset{r}{\Σ}}{a}_i+\frac{1}{{(r-1)}^2}\underset{i=1}{\overset{r}{\Σ}}\underset{k=1}{\overset{r}{\Σ}}(a_i.*a_k),(i\≠q,k\≠q,q=1,...,r)---\left(7\right)$

步骤1.2：利用PCA算法将数据降到r维，确定顶点集{ω₁,a₂,...,a_r}，{a₁,ω₂,...,a_r}，…，{a₁,a₂,...,ω_r}在该空间中所张成的r个r-1维超平面H₁…H_r；

步骤1.3：求解线性方程组(8)，非线性端点p即为其唯一解：

$\det \left(\left[\begin{array}{ccccc}1& 1& ...& 1& 1\\ U^T{\mathrm{\ω}}_1& U^T{a}_2& ...& U^T{a}_r& p\end{array}\right]\right)=0,...,\det \left(\left[\begin{array}{ccccc}1& 1& ...& 1& 1\\ U^T{a}_1& U^T{a}_2& ...& U^T{\mathrm{\ω}}_r& p\end{array}\right]\right)=0---\left(8\right)$

其中，U^T为数据前r个特征向量构成的矩阵，det表示行列式；

输出非线性端点p；

(2)修正投影位置并迭代估计丰度

先利用p点计算投影坐标s′_j，然后代入模型计算所对应的非线性混合部分，再确定使像元到其所在直线距离最短的系数，从而得到修正后的新投影位置，最后用FCLS算法求解新投影位置的丰度，并以此迭代进行直到收敛，具体流程为：

输入：高光谱数据X∈R^n×m和端元矩阵A∈R^n×r，非线性端点p；

输出：丰度矩阵S∈R^r×m；

步骤2.1：以最小二乘法求解优化问题(9)得投影重心坐标s′_j

$\begin{array}{cccc}\min & \frac{1}{2}\left|\right|U^T{x}_j-\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{U}^T{a}_i{h}_{i,j}-{\mathrm{ph}}_{r+1,j}|{|}^2,& s.t.& \underset{i=1}{\overset{r+1}{\Σ}}{h}_{i,j}=1\end{array}---\left(9\right)$

这里，

步骤2.2：对于每一个观测像元x_j,j＝1,...,m执行如下操作：

2.2a)初始化迭代次数t＝0和丰度

2.2b)计算对应投影点的BMM非线性混合向量

对于GBM和FM，定义为

对于PPNM，定义为：

2.2c)求解最优化问题(10)得系数并以式(11)计算新投影位置

$\begin{array}{cc}min& \frac{1}{2}\left|\right|x_j-\underset{i=1}{\overset{r}{\mathrm{\Σ}}}{a}_i{s}_{i,j}^{\left(t\right)}-{\mathrm{\λ}}_j^{\left(t\right)}{\hat{x}}_j^{\left(t\right)}|{|}^2\end{array}---\left(10\right)$

$y_j^{\left(t\right)}=x_j-{\mathrm{\λ}}_j^{\left(t\right)}{\hat{x}}_j^{\left(t\right)}---\left(11\right)$

2.2d)用线性丰度估计算法求解的丰度s^(t+1)；

2.2e)如达到最大迭代次数或满足收敛精度则输出丰度矩阵S＝S^(t+1)，否则回到2.2b)；

步骤2.3：输出结果S＝(s₁,s₂,…,s_m)，并用线性规划求得非线性参数γ与b。

2.根据权利要求1所述的解混方法，其特征在于：步骤2.2c)中的系数为最小二乘解：

${\mathrm{\λ}}_j^{\left(t\right)}={(x_j-{\mathrm{\Σ}}_{i=1}^r{a}_i{s}_{i,j}^{\left(t\right)})}^T{\hat{x}}_j^{\left(t\right)}/\left({\hat{x}}_j^{\left(t\right)T}{\hat{x}}_j^{\left(t\right)}\right).$

3.根据权利要求1所述的解混方法，其特征在于：步骤2.2d)中采用FCLS算法求解s^(t+1)。

4.根据权利要求1所述的解混方法，其特征在于：步骤2.2e)中的最大迭代次数为200。

5.根据权利要求4所述的解混方法，其特征在于：步骤2.2e)中，通过判断每次迭代与前一次迭代之间的相对变化是否小于10^-8来判断是否收敛。