正定矩阵因子分解法（PMF）

3.2.4.1 方法建立

就全国范围而言，我国地下水质量总体较好，根据国家《地下水质量标准》（GB/T 14848—93），我国63%的地区地下水可直接饮用，17%经适当处理后可供饮用，12%不宜饮用，剩余8%为天然的咸水和盐水，由此可见，不宜饮用的地下水和天然咸水、盐水占到了20%，对于这些地下水型水源地饮用水指标并不一定受到污染而存在超标现象，其水质可能受到地下水形成演化影响更为明显，因此，考虑选择反映地下水形成、演化的地下水水化学类型常规指标，进行影响因素解析。地下水水质指标在取样与分析过程中，由于取样和样品处理、试剂和水纯度、仪器或让量度和仪器洁净、采用的分析方法、测定过程以及数据处理等过程均会产生测量误差（系统误差，随机误差，过失误差）。从取样到分析结果计算误差都绝对存在，虽然在各个过程中进行质量控制，但无法完全消除不确定性的影响，为确保分析结果的可靠性，采用PMF法对地下水水质指标考虑一定的不确定性误差，使分析数据能够准确地反映实际情况。

PMF（Positive Matrix Factorization）与主成分分析（PCA）、因子分析（FA）都是利用矩阵分解来解决实际问题的分析方法，在这些和锋方法中，原始的大矩阵被近似分解为低秩的V＝WH形式。但PMF与PCA和FA不同，PCA、FA方法中因子W和H中的元素可为正或负，即使输入的初始矩阵元素全是正的，传统的秩削减算法也不能保证原始数据的非负性。在数学上，从计算的观点看，分解结果中存在负值是正确的，但负值元素在实际问题中往往是没有意义的。PMF是在矩阵中所有元素均为非负数约束条件之下的矩阵分解方法，在求解过程中对因子载荷和因子得分均做非负约束，避免矩阵分解的结果中出现负值，使得因子载荷和因子得分具有可解释性和明确的物理意义。PMF使用最小二乘方法进行迭代运算，能够同时确定污染源谱和贡献，不需要转换就可以直接与原始数据矩阵作比较，衫棚局分解矩阵中元素非负，使得分析的结果明确而易于解释，可以利用不确定性对数据质量进行优化，是美国国家环保局（EPA）推荐的源解析工具。

3.2.4.2 技术原理

PMF：模型是一种基于因子分析的方法，具有不需要测量源指纹谱、分解矩阵中元素非负、可以利用数据标准偏差来进行优化等优点。目前PMF模型此方法成功用于大气气溶胶、土壤和沉积物中持久性有毒物质的源解析，已有成熟的应用模型 PMF1.1，PMF2.0，PMF3.0等。PMF模型基本方程为：

Xnm＝GnpFpm+E （3.7）

式中：n——取样点数；

m——各取样点测试的成分数量；

p——污染源个数；

Xnm——取样点各成分含量；

Gnp——主要源的贡献率；

Fpm——源指纹图谱。

基本计算过程如下：

1）样品数据无量纲化，无量纲化后的样品数据矩阵用D表示。

2）协方差矩阵求解，为计算特征值和特征向量，可先求得样品数据的协方差矩阵，用D′为D的转置，算法为：

Z＝DD′ （3.8）

3）特征值及特征向量求解，用雅各布方法可求得协方差矩阵Z的特征值矩阵E和特征向量矩阵Q，Q′表示Q的转置。这时，协方差矩阵可表示为：

Z＝QEQ′ （3.9）

4）主要污染源数求解，为使高维变量空间降维后能尽可能保留原来指标信息，利用累计方差贡献率提取显著性因子，判断条件为：

地下水型饮用水水源地保护与管理：以吴忠市金积水源地为例

式中：n——显著性因子个数；

m——污染物个数；

λ——特征值。

5）因子载荷矩阵求解，提取显著性因子后，利用求解得到的特征值矩阵E和特征向量矩阵Q进一步求得因子载荷矩阵S和因子得分矩阵C，这时，因子载荷矩阵可表示为：

S＝QE1/2 （3.11）

因子得分矩阵可表示为：

C＝（S′S）-1S′D （3.12）

6）非负约束旋转，由步骤5求得的因子载荷矩阵S和因子得分矩阵C分别对应主要污染源指纹图谱和主要污染源贡献，为解决其值可能为负的现象，需要做非负约束的旋转。

7）首先利用转换矩阵T1对步骤5求得的因子载荷矩阵S和因子得分矩阵C按下式进行旋转：

地下水型饮用水水源地保护与管理：以吴忠市金积水源地为例

C1＝T1C （3.14）

式中：S1——旋转后的因子载荷矩阵；

C1——旋转后的因子得分矩阵；

T1——转换矩阵，且T1＝（C∗C′）（C∗C′）-1（其中：C∗为把C中的负值替换为零后的因子得分矩阵）。

8）利用步骤7中旋转得到的因子载荷矩阵S1构建转换矩阵T2对步骤5中旋转得到的因子载荷矩阵S1和因子得分矩阵C1继续旋转：

S2＝S1T2 （3.15）

地下水型饮用水水源地保护与管理：以吴忠市金积水源地为例

式中：S2——二次旋转后的因子载荷矩阵；

C2——二次旋转后的因子得分矩阵；

T2——二次转换矩阵，且T2＝（S′1+S1）-1（S′1+ ）（其中： 为S1中的负值换为零后的因子载荷矩阵）。

9）：重复步骤7、8，直到因子载荷中负值的平方和小于某一设定的误差精度e而终止，最终得到符合要求的因子载荷矩阵S，即主要污染源指纹图谱。

3.2.4.3 方法流程

针对受体采样数据直接进行矩阵分解，得到各污染源组分及其贡献率的统计方法（图3.5）。

图3.5 方法流程图

（1）缺失值处理

正定矩阵因子分析是基于多元统计的分析方法，对数据有效性具有一定的要求，因此在进行分析之前首先对数据进行预处理。根据已有数据的特征结合实际情况主要有以下5种处理方法。

1）采样数据量充足的情况下直接丢弃含缺失数据的记录。

2）存在部分缺失值情况下用全局变量或属性的平均值来代替所有缺失数据。把全局变量或是平均值看作属性的一个新值。

3）先根据欧式距离或相关分析来确定距离具有缺失数据样本最近的K个样本，将这K个值加权平均来估计该样本的缺失数据。

4）采用预测模型来预测每一个缺失数据。用已有数据作为训练样本来建立预测模型，如神经网络模型预测缺失数据。该方法最大限度地利用已知的相关数据，是比较流行的缺失数据处理技术。

5）对低于数据检测限的数据可用数据检测限值或1/2检测限以及更小比例检测限值代替。

（2）不确定性处理

计算数据不确定性。

地下水型饮用水水源地保护与管理：以吴忠市金积水源地为例

式中：s——误差百分数；

c——指标浓度值；

l——因子数据检出限。

（3）数据合理性分析

本研究所用数据在放入模型前以信噪比S/N（Signal to Noise）作为标准进行筛选，信噪比S/N为：

地下水型饮用水水源地保护与管理：以吴忠市金积水源地为例

式中：xij——第i采样点第j个样品的浓度；

sij——第i采样点第j个样品的标准偏差。

信噪比小，说明样品的噪声大，信噪比越大则表示样品检出的可能性越大，越适合模型。

（4）数据输入及因子分析

与其他因子分析方法一样，PMF不能直接确定因子数目。确定因子数目的一般方法是尝试多次运行软件，根据分析结果和误差，Q值以及改变因子数目时Q值的相对变化等来确定合理的因子数目。

3.2.4.4 适用范围

PMF对污染源和贡献施加了非负限制，并考虑了原始数据的不确定性，对数据偏差进行了校正，使结果更具有科学的解释。PMF使用最小二乘方法，得到的污染源不需要转换就可以直接与原始数据矩阵作比较，PMF方法能够同时确定污染源和贡献，而不需要事先知道源成分谱。适用于水文地质条件简单，观测数据量较大，污染源和污染种类相对较少的地区，运用简便，可应用分析软件进行计算。

3.2.4.5 NMF 源解析

NMF在实现上较PMF算法简单易行，非负矩阵分解根据目的的不同大致可以分为两种：一是在保证数据某些性质的基础上，将高维空间的样本点映射到某个低维空间上，除去一些不重要的细节，获得原数据的本质信息；二是在从复杂混乱的系统中得到混合前的独立信息的种类和强度。因此，基于非负矩阵分解过程应用领域的不同，分解过程所受的约束和需要保留的性质都不相同。本书尝试性地将NMF算法应用于水质影响因素的分离计算中（表3.2）。

表3.2 RMF矩阵分解权值表

依照非负矩阵分解理论的数学模型，寻找到一个分解过程V≈WH，使WH和V无限逼近，即尽可能缩小二者的误差。在确保逼近的效果，定义一个相应的衡量标准，这个衡量标准就叫作目标函数。目标函数一般采用欧氏距离和散度偏差来表示。在迭代过程中，采用不同的方法对矩阵W和H进行初始化，得到的结果也会不同，算法的性能主要取决于如何对矩阵W和H进行初始化。传统的非负矩阵算法在对矩阵W和H赋初值时采用随机方法，这样做虽然简单并且容易实现，但实验的可重复性以及算法的收敛速度是无法用随机初始化的方法来控制的，所以这种方法并不理想。许多学者提出改进W和H的初始化方法，并发展出专用性比较强的形式众多的矩阵分解算法，主要有以下几种：局部非负矩阵分解（Local Non-negative Matrix Factorization，LNMF）、加权非负矩阵分解（Weighted Non-negative Matrix Factorization，WNMF）、Fisher非负矩阵分解（Fisher Non-negative Matrix Factorization，FNMF）、稀疏非负矩阵分解（Sparse Non-negative Matrix Factorization，SNMF）、受限非负矩阵分解（Constrained Non-negative Matrix Factorization，CNMF）、非平滑非负矩阵分解（Non-smooth Non-negative Matrix Factorization，NSNMF）、稀疏受限非负矩阵分解（Nonnegative Matrix Factorization with Sparseness Constraints，NMF-SC）等理论方法，这些方法针对某一具体应用领域对NMF算法进行了改进。

本书尝试应用MATLAB工具箱中NNMF程序与改进的稀疏非负矩阵分解（SNMF）对研究区11项指标（同PMF数据）进行分解，得到各元素在综合成分中的得分H，初始W0，H0采用随机法取初值。r为分解的基向量个数，合适的r取值主要根据试算法确定，改变r值观察误差值变化情况，本书利用SMNF算法计算时，r分别取2，3，4，采用均方误差对迭代结果效果进行评价，结果显示当r取2，4时误差值为0.034，取3时误差值为0.016，因此r＝3是较合理的基向量个数。采用NNMF算法进行计算时，利用MATLAB工具箱提供的两种计算法分别进行计算，乘性法则（Multiplicative Update Algorithm）计算结果误差项比最小二乘法（Alternating Least-squares Algorithm）计算误差值小且稳定，但总体NNMF计算误差较大，改变初始W0，H0取值和增加迭代次数误差均未明显减小，调整r取值，随着r值的增大误差逐渐减小。

对比SNMF和NNMF算法所得权值结果，两种方法所得权值趋势一致，但得分值有所不同，由于SNMF算法对矩阵进行了稀疏性约束，计算结果中较小的权值更趋近于0，两次结果中在三个基向量上总体权值较大的元素项为T-Hard、 、Mg2+、Ca2+、 ，从盲源分离的角度来看该几种元素对地下水具有较大的影响，但从地下水水质影响因素来看，该方法对数据的分析偏重于突出局部数据的特征，在各因素相关性较大但含量不高的情况下，容易忽略了关键的影响因素。从权值得分来看，SNMF法解析的第一个基向量上的元素包括EC、T-Hard、NH4—N、 、 、TDS；第二基向量主要有Na+、Mg2+、Cl-；第三个基向量 、Ca2+，从结果可以看出该方法进行矩阵分解并未得到可合理解释的源项结果，方法有待进一步研究及验证。

得到非负矩阵分解（NMF）的基图像矩阵和系数矩唯誉衡阵后，Hoyer [1] 提出可以利用L1范数和L2范数之间的差异度来衡量分解后矩阵的稀疏度。

L1范数 ：所有元素的绝对值之和。

L2范数 ：所有元素指做的平方之和的平方根。

计算虚陪x的稀疏度，分三步：

稀疏度计算（matlab）

在ORL人脸数据集上的实验对比：

Reference:

[1] Non-negative Matrix Factorization with Sparseness Constraints

21．MATLAB支持的几种图像文件格式：

⑴JPEG（Joint Photogyaphic Expeyts Group）：一种称为联合图像专家组的图像压缩格式。

⑵BMP（Windows Bitmap）：有1位、4位、8位、24位非压缩图像，8位RLE（Run length Encoded）的图像。文件内容包括文件头（一个BITMAP FILEHEADER数据结构）、位图信息数据块（位图信息头BITMAP INFOHEADER和一个颜色表）和图像数据。

⑶PCX（Windows Paintbrush）：可处理1位、4位、8位、16位、24位等图像数据。文件内容包括文件头、图像数据和扩展色图数据。

⑷TIFF（Tagged Iamge File Format）：处理1位、4位、8位、24位非压缩图像，1位、4位、8位、24位packbit压缩图像，1位CCITT压缩图像等。文件内容包括文件头、参数指针表与参数域、参数数据表和图像数据四部分。

⑸PNG（Portable Network Graphics）：包括1位、2位、4位、8位和16位灰度图像，8位和16位索引图像，24位和48位真彩色图像。

⑹GIF（Graphics Interchange Format）：任何1位到8位的可交换的图像。

⑺HDF（Hierarchial Data Format）：有8位、24位光栅图像数据集。

⑻ICO（Windows Icon resource）：有1位、4位、8位非压缩图像。

⑼CUR（Windows Cursor resource）：有1位、4位、8位非压缩图像。

⑽XWD（X Windows Dump）：包括1位、8位Zpixmaps，XYBitmaps，XYPixmmmaps。

⑾RAS（Sun Raster image）：有1位bitmap、8位索引、24位真彩色和带有透明度的32位真彩色。

⑿PBM（Portable Bitmap）。

⒀PGM（丛尘唯Portable Graymap）。

⒁PPM（Portable Pixmap）。

2．MATLAB支持五种图像类型，即二值图像、索引图像、灰度图像、RGB图像和多帧图像阵列。

图像类型

Double数据

uint8和uint16数据

二值图像

图像为m×n的整数矩阵，元素值范围[0,1]

图像为m×n的整数矩阵，元素值范围[0,1]

索引图像

图像为m×n的整数矩阵，元素值范围[0,p]

图像为m×n的整数矩阵，元素值范围[0,p-1]

灰度图像

图像为m×n的浮点数矩阵，元素值范围[0,1]

图像为m×n的整数矩阵，元素值范围[0,255]或[0,65535]

RGB图像

图像为m×n×3的浮点数矩阵，元素值范围[0,1]

图像为m×n×3的整数矩阵，元素值范围[0,255]或[0,65535]

其中，多帧图像阵列是由多帧图像组成的，每一帧图像可以为前四种图像中的一种，但组成一个多帧图像阵列的图像必须为同一种。cat函数可以将具有相同尺寸的几个独立图像存成多帧文件。对于多帧图像也可以从中提取单帧。

2.2.2 图像类型判断及转换

1．在MATLAB中如果要判断一个图像文件的类型，可使用如下指令：

l isbw：若图像为二值图像，则返回真。

l isgray：若图像为灰度图像，则返回真。

l isind：若图像为索引图像，则返回真。

l isrgb：若图像为RGB图像，则返回真。

2．在MATLAB系统中，要将一个类型的图像文件转换成另一个类型的图像文件，只需将前一个文件的图像数据用imread读出，再用imwrite以适当的格式写到后一个图像文件中去即可。另外，要将灰度图像I转换成RGB图像，可运用cat指令。cat函数可以把一些单一的图像合并成图像兄空序列。在图像序列中每个图像需有相同的大小，如果渗培是索引图像，色图也要一致。

此外，MATLAB还提供了若干函数，用于图像类型的转换。这些函数有：

l dither：用抖动法（dithering）转换图像。该函数通过颜色抖动（颜色抖动即改变边沿像素的颜色，使像素周围的颜色近似于原始图像的颜色，从而以空间分辨率来换取颜色分辨率）来增强输出图像的颜色分辨率。该函数可以把RGB图像转换成索引图像或把灰度图像转换成二值图像。

l gray2ind：灰度图像或二值图像向索引图像转换。

l grayslice：设定阈值将灰度图像转换为索引图像。

l im2bw：设定阈值将灰度、索引、RGB图像转换为二值图像。

l im2double：将图像数组转换为double型。

l im2uint8：将图像数组转换为uint8型。

l im2uint16：将图像数组转换为uint16型。该函数不支持二值图像序列的转换。

l ind2gray：索引图像向灰度图像转换。

l ind2rgb：索引图像向RGB图像转换。

l mat2gray：将一个数据矩阵转换为灰度图像。

l rgb2gray：RGB图像向灰度图像转换或将彩色色图转换成灰度色图。

l rgb2ind：RGB图像向索引图像转换。包含三种不同方法：均衡量化、最小值量化、色图映射。

l im2java：一般图像向Java图像转换。

l label2rgb：标志图像向RGB图像转换。

2.2.3 图像的查询及读写

在MATLAB中要查询一个图像文件的信息，只要用imfinfo指令加上文件及其完整路径名即可。函数调用格式为：

info = imfinfo(filename,fmt)

info = imfinfo(filename)

参数fmt对应于所有图像处理工具箱中所有支持的图像文件格式。

MATLAB提供了两个重要的用于图像文件的读写的指令，分别是从图像文件中读取数据的imread，以及将数据写入到图像文件中的imwrite。

1．imread的常见调用格式为：

A = imread(filename,fmt)

其作用是将文件名用字符串filename表示的，扩展名用fmt表示的图像文件中的数据读到矩阵A中。如果filename所指的为灰度级图像，则A为一个二维矩阵；如果filename所指的为RGB图像，则A为一个m×n×3的三维矩阵。Filename表示的文件名必须在MATLAB的搜索路径范围内，否则需指出其完整路径。

imread的其他几种重要的调用格式为：

[X,map] = imread(filename.fmt)

[…] = imread(filename)

[…] = imread(URL,…)

[…] = imread(…,idx) (CUR,ICO and TIFF only)

[…] = imread(…,’frames’,idx)(GIF only)

[…] = imread(…,ref) (HDF only)

[…] = imread(…,’BackgroundColor’,BG)(PNG only)

[A,map,alpha] = imread(…) (ICO,CUR and PNG only)

上面一些参数的含义如下：idx是指读取图标（cur、ico、tiff）文件中第idx个图像，默认值为1。’frame’,idx是指读取gif文件中的图像帧，idx值可以是数量、向量或’all’。ref是指整数值。alpha是指透明度。

2．imwrite的常用调用格式为：

imwrite(A,filename,fmt)

imwrite(X,map,filename,fmt)

imwrite(…,filename)

imwrite(…,Param1,Val1,Param2,Val2…)

其中imwrite(…,Param1,Val1,Param2,Val2…)可以让用户控制HDF、JPEG、TIFF等一些图像文件格式的输出特性。

在MATLAB中，默认的保存类型是uint8。由于PNG和TIFF格式支持16位的图像，所以保存这类图像时，保存类型就是uint16。

2.2.4 图像的显示

显示图像的最基本的手段是使用image函数。该函数还产生了图像对象的句柄，并允许对对象的属性进行设置。此外，imagesc函数也具有image的功能，所不同的是imagesc函数还自动将输入数据比例化，以全色图的方式显示。

imshow函数比image和imagesc更常用，它能自动设置句柄图像的各种属性。imshow可用于显示各类图像。对于每类图像，调用方法如下：

l imshow filename：显示图像文件。

l imshow(BW)：显示二值图像，BW为黑白二值图像矩阵。

l imshow(X,map)：显示索引图像，X为索引图像矩阵，map为色彩图示。

l imshow(I)：显示灰度图像，I为二值图像矩阵。

l imshow(RGB)：显示RGB图像，RGB为RGB图像矩阵。

l imshow(I,[low high])：将非图像数据显示为图像，这需要考虑数据是否超出了所显示类型的最大允许范围，其中[low high]用于定义待显示数据的范围。

有关图像显示的函数或其辅助函数，除了上述的以外，MATLAB还提供了一些用于进行图像的特殊显示的函数。

l colorbar：为图像的显示增加一个颜色条，这一用法对于了解被显示图像的灰度级别特别有用。

l getimage：获取图像数据。

l immovie：将多帧索引图像制作成连续图像格式。其调用格式为：

mov = immovie(D,map)

不过这种功能只对索引图像有效，其中D为多帧索引图像阵列，map为索引图像的对应色阶。对于其他类型图像，则需要首先将其转换为索引图像。

l montage：多帧图像的一次显示。它能将每一帧分别显示在一幅图像的不同区域，所有子区的图像都用同一个色彩条。

l movie：播放多帧连续图像。

l subimage：在一个图形区域内显示多个图像。

l truesize：调整图像显示的尺寸。

l warp：显示图像的纹理表面图。前面提到的图像显示手段都只能在二维平面上显示，MATLAB6.5的一个强大功能是能将平面图像显示在空间三维曲面上。这是由warp函数的纹理成图功能来实现的，该功能能通过双线性插值将平面图像投影到三维曲面上。

l zoom：将图像或二维图形进行放大或缩小显示。zoom本身是是一个开关键，zoom on用于打开缩放模式，zoom off用于关闭该模式，zoom in用于放大局部图像，zoom out用于缩小图像。

至于多个图像的显示，则可分为两个方面：在不同的图形窗口显示不同的图像，可以用figure指令来实现；在同一个图形窗口显示多图，可以用subplot来实现。

第三章 图像运算

3.1 图像的点运算

点运算将输入图像映射为输出图像，输出图像每个像素点的灰度值仅由对应的输入像素点的值决定。它常用于改变图像的灰度范围及分布，是图像数字化及图像显示的重要工具。点运算因其作用性质，也被称为对比度增强、对比度拉伸或灰度变换。在真正进行图像处理之前，有时可以用点运算来克服图像数字化设备的局限性。

点运算实际上是灰度到灰度的映射过程。点运算不会改变图像内像素点之间的空间关系。设输入图像为A(x, y)，输出图像为B(x, y)，则点运算可表示为：

B(x, y) = f [A(x, y)](3-1)

点运算可完全由灰度变换（gray-scale transformation）函数s=f (r)决定，后者描述了输入灰度级与输出灰度级之间的映射关系。图像的点运算分为线性点运算和非线性点运算两种。

3.1.1 线性点运算

线性点运算是指灰度变换函数f为线性函数时的运算。

如图3.1，当a＞1时，输出图像对比度增大；当a＜1时，输出图像对比度降低；当a＝1,b＝0时，输出图像就是输入图像的简单复制；当a＝1,b≠0时，仅使输出图像的灰度值上移或下移，其效果是使整个图像更亮或更暗。如果a为负值，暗区域将变亮，亮区域将变暗，点运算完成了图像求补。

除了调节对比度以外，还有一种典型的线性点运算的应用就是灰度标准化。设灰度图像为I[W][H]，其中W表示图像宽度，H表示图像的高度，那么灰度图像的平均灰度和方 图3.1 线性函数

差由如下计算公式得到：

平均灰度：

(3-2)

方差：

(3-3)

可以将其变换为具有相同均值和方差的变换函数（线性映射），其形式如下：

(3-4)

其中σ0和u 0为给定的变换参数。灰度标准化可以用来生成一些常用的平均模型。

3.1.2 非线性点运算

非线性点运算对应于非线性映射函数，典型的映射包括平方函数、对数函数、截取函数（窗口函数）、域值函数、多值量化函数等。

阈值化处理是最常用的一种非线性点运算，它的功能是选择一阈值，将图像二值化，然后使用生成的二进制图像进行图像分割及边缘跟踪等处理。

直方图均衡化也是一种非常常用的非线性点运算。它是指将一个已知灰度分布的图像使用某种非线性灰度变换函数进行计算，使运算结果变成一幅具有均匀灰度分布的新图像。经过直方图均衡化的点运算处理后，实际的直方图将呈现参差不齐的外形，这是由于灰度级的可能个数是限造成的。在一些灰度级处可能没有像素，在另外一些灰度级处则像素很拥挤。

点运算的MATLAB实现：

I=imread('rice.png')

rice=double(I)

rice2=rice*0.5+50

J=uint8(rice2)

subplot(1,2,1),imshow(I)

subplot(1,2,2),imshow(J)

转载的地址忘啦。呵呵

---