CN114329333A - 航空γ能谱数据NASVD降噪有效奇异值数目确定方法 - Google Patents

Abstract

本发明涉及一种航空γ能谱数据NASVD降噪有效奇异值数目确定方法，构建输入谱线数据矩阵后，利用NASVD算法得到主成分矩阵和幅度矩阵，接着利用“滑窗法”对获取的谱线主成分进行分类，对分类后的谱线主成分进行谱线重构，得到系列重构谱线。然后对系列重构谱线中相邻谱线之间的相关性进行分析计算，得到重构谱线的系列相关系数。以相关系数的绝对值作为判断依据，确定有效奇异值的数目。本发明通过一种简单、高效的方法能够快速确定航空γ能谱数据NASVD降噪中有效奇异值的数目，提高了降噪质量、缩短了处理时间周期，避免人工靠经验或推荐值确定有效奇异值数目产生过降噪或产生结果系统性误差，避免返工和重复计算，降低操作者的工作强度。

Description

航空γ能谱数据NASVD降噪有效奇异值数目确定方法

技术领域

本发明涉及一种航空γ能谱降噪方法，具体地说是一种航空γ能谱数据NASVD降噪有效奇异值数目确定方法。

背景技术

航空γ能谱测量是根据天然或人工放射性核素的γ射线能量差异，利用航空γ能谱测量系统确定地表岩石、土壤和大气中的钾、铀、钍和其它放射性核素含量及其分布的一种方法，是开展战略性矿产勘查、环境辐射监测与评价、核设施监测、核应急航空监测的一种重要技术手段，在众多应用领域中发挥着重要作用。近年来随着资源与环境对航空测量的急速增长，尤其对航空γ能谱测量技术提出了新的更高的要求。

γ能谱降噪是航空γ能谱测量中的关键技术之一，是提高放射性核素含量估算精度和减小谱定性定量分析误差的关键步骤，多年里一直是国内外研究者的研究热点。目前，航空γ能谱采用的最主要、最有效的谱线降噪方法是国际原子能机构IAEA-1363报告中推荐的NASVD（Noise Adjusted Singular Value Decomposition，噪声调整奇异值分解）。NASVD是一种基于多元统计分析的多道能谱降噪方法，即通过主成分分析技术从原始γ能谱数据集中提取出主成分谱，通过选择低阶主成分谱重构谱线来实现降噪。自1997年NASVD方法提出以来，国内外学者对其降噪效果以及与传统方法比对进行了深入研究，大量研究表明该方法不仅能够显著地减少γ能谱中的噪声，而且在克服传统降噪方法不足方面具有显著的优势，因而在国内外大量的航测项目中得到广泛应用。

NASVD方法降噪的关键是如何选取最佳的低阶主成分个数（即有效奇异值数目）来重构谱线。有效奇异值如果选择个数过少，则有可能将能谱信号分量当作噪声去掉，造成重要能谱异常的丢失；选择个数过多，可能导致能谱信号中噪声信息成分过多，造成重要微弱异常淹没在噪声中而不能被有效识别。因此，有效奇异值数目选取至关重要，选取不当将可能造成整个测区或区域航测结果分析和解释的偏差，由此可能造成不可挽回的损失，尤其是拥有海量γ能谱数据的重特大航测项目。

目前，有效奇异值数目确定方法主要是通过观察主成分谱线形状进行人工判断。在“国际原子能机构IAEA-1363报告”中，关于有效奇异值个数选取问题仅提出不能仅依靠特征值的大小进行判断，同时还要考虑参考空间信息，且也仅给出了推荐值。在“EJ/T1032-2005航空伽马能谱测量规范”中，关于自行设计的滤波器，要求滤波类型和参数的选择要经过反复试验进行确定。查阅的国内外相关文献中，涉及内容大多数是降噪效果比对及改进方法，而关于如何选取有效奇异值个数内容，只简单描述为观察确定，而未给出确定的理论依据。在国际专业处理软件Praga4中，开发了NASVD方法处理模块，但有效奇异值数目还是需要依靠观察主成分谱线形状由人工确定。

而仅靠人工判断方法在项目实践中却存在很大不足，主要表现在：

1、易发生过降噪现象。降噪效果取决于操作者的工作经验，人工判断缺乏理论依据，难以避免过降噪现象的发生，往往造成重要的航测成果可能丢失或难以识别而被忽略。

2、易发生结果系统性偏差，造成项目降噪处理工作全部返工。项目实践中，NASVD方法通常采用一条/多条测线、单/多架次或单/多/全部测区数据集进行降噪处理，涉及成百万甚至几千万条谱线数据。由于测区环境的复杂性，导致噪声信号大量存在，原始信号的主成分谱线形状往往难以区分，致使有效奇异值数目难以确定，仅依靠人工判断达到最佳降噪效果，对操作人员专业素质要求高且工作效率十分低下，易产生结果系统性偏差且难以溯源，甚至会造成项目降噪处理工作全部返工。

3、重复计算工作量大，操作人员工作强度高，降噪质量难以保证。

4、自动化程度低，处理时间周期长。难以满足堪称海量航空γ能谱数据快速处理和分析的要求。目前，运用新技术对历史航测资料的二次开发涉及近500万平方公里的海量能谱数据，工作量巨大且不能满足数据快速处理和分析的要求。

所以，如果能够通过一种方法准确确定航空γ能谱数据NASVD降噪有效奇异值的数目，将有效避免由于人为因素造成的过降噪现象的发生，提高NASVD方法降噪的自动化水平和工作效率，进一步保障放射性核素含量估计精度和减小航空γ能谱定性定量分析的误差。

发明内容

本发明的目的就是提供一种航空γ能谱数据NASVD降噪有效奇异值数目确定方法，以解决现在航空γ能谱数据NASVD降噪中人工判断有效奇异值数目容易发生过降噪现象、产生结果系统性偏差、重复计算工作量大以及降噪质量难以保证的问题。

本发明是这样实现的：一种航空γ能谱数据NASVD降噪有效奇异值数目确定方法，包括以下步骤。

a.读入原始航空γ能谱数据，构建输入谱线数据矩阵A。

b.利用NASVD算法对输入谱线数据矩阵A进行计算，得到主成分矩阵S和幅度矩阵B，主成分矩阵S的各个列向量依次为各主成分，且各主成分按照对应特征值的大小降序排列。

c.对主成分按照前1个、前2个、前3个…前n个进行分类。

d.分别对前1个、前2个、前3个…前n个主成分进行重构，分别得到重构谱矩阵A₁，A₂，A₃，…A_n。

e.依次计算相邻两个重构谱矩阵的相关系数|r_i|。

f.判断|r_i|的大小，当|r_i|第一次等于1时，i值即被确定为有效奇异值的数目。

进一步地，在步骤b中，主成分矩阵S和幅度矩阵B的求解步骤如下。

1）对输入谱线数据矩阵A各行元素求和，公式如下：

R=(r₁₁，r₂₁，… r_m1)^T

式中，R表示输入谱线数据矩阵A各行元素之和组成的列向量，r_i1表示输入谱线数据矩阵A中第i行元素之和，其中i=1，2，…，m。

2）对输入谱线数据矩阵A各列元素求和，公式如下：

C=(c₁₁，c₁₂，…，c_1n)

式中，C表示输入谱线数据矩阵A中各列元素之和组成的行向量，c_1j表示输入谱线数据矩阵A中第j列元素之和，其中j=1，2，…，n。

3）归一化单位谱，公式如下：

式中，N表示对C进行归一化处理后的行向量。

4）求取单位方差矩阵，公式如下：

式中，A_new表示m×n单位方差矩阵。

5）对A_new进行奇异值分解，公式如下：

A_new=UΣV^T

式中，U表示m×n正交矩阵，V表示n×n正交矩阵，Σ表示由奇异值λ组成的n×n对角矩阵，Σ=diag(λ₁，λ₂，…，λ_n)，且λ₁≥λ₂≥…λ_n≥0。

6）计算主成分矩阵和幅度矩阵，公式如下：

式中，S表示n×n主成分矩阵，其各个的列向量依次为各主成分，按照对应特征值的大小降序排列，diag(N)表示由向量N组成的对角矩阵。

并计算幅度矩阵，公式如下：

式中，B表示m×n幅度矩阵，diag(R)表示由向量R组成的对角矩阵。

进一步地，在步骤d中，首先重构输入谱线数据矩阵A，公式如下：

式中，b_i表示B中的第i个列向量，s_i ^T表示S^T中的第i个行向量。

按照步骤c中的分类，选取主成分矩阵S中代表信号的前k个主成分，其中k=1，2，3，...，n，得到重构的航空γ能谱数据集A_k：

当k=1时，得到重构谱矩阵A₁；当k=2时，得到重构谱矩阵A₂，依次类推，得到重构谱矩阵A₁，A₂ ，A₃，…A_n。

进一步地，在步骤e中，利用相关系数计算公式，分别计算A_i和A_i+1的相关系数r_i，其中，i=1，2，…，n-1。

进一步地，在步骤e中，如果相关系数第一次等于1之后，后续|r_i|的值在0.99990~1.000000值之间变化，则航空γ能谱数据符合测量误差要求，可继续进行后续操作；否则，说明测量设备可能发生的能谱峰漂超出了峰漂允许的误差范围，需要对能谱数据进行峰漂修正或对航测设备进行检查。

本发明首先对实测的原始航空γ能谱数据进行NASVD方法计算，获取实测谱线中的相互正交的谱线主成分，并按各自对谱线形状贡献的大小进行降序排列，接着利用“滑窗法”对获取的谱线主成分进行分类，对分类后的谱线主成分进行谱线重构，得到系列重构谱线。然后对系列重构谱线中相邻谱线之间的相关性进行分析计算，得到重构谱线的系列相关系数。以相关系数的绝对值作为判断依据，当相关系数的绝对值第一次等于1时，即可明确确定有效奇异值的数目i。最后采用前i个低阶主成分进行谱线重谱，最终得到降噪后的谱线。

主成分采用“滑窗法”进行分类，其基本思想是依据信号相关而噪声不相关的原则，将代表信号的低阶主成分依次叠加后进行谱重构，求取其相邻重构谱的相关系数，随着不断叠加，相关系数不断增大，直至等于1，由此可以得到叠加低阶主成分时的序号，即确定有效奇异值的数目。

本发明通过一种简单、高效的方法能够快速确定航空γ能谱数据NASVD降噪中有效奇异值的数目，提高了降噪质量、缩短了处理时间周期，避免人工靠经验或推荐值确定有效奇异值数目产生过降噪或产生结果系统性误差，从而避免返工和重复计算，降低操作者的工作强度。

附图说明

图1是本发明的流程图。

图2是本发明采用滑窗法对主成分进行分类的示意图。

图3是本发明静态实测数据集相关系数计算结果。

图4是本发明静态实测数据集的前6个主成分。

图5是本发明静态实测数据集的前50个奇异值。

图6是本发明飞行实测数据集相关系数计算结果。

图7是本发明飞行实测数据集的前11个主成分。

图8是本发明飞行实测数据集的前50个奇异值。

图9是本发明静态实测数据集降噪前后的谱线比较。

具体实施方式

如图1所示，本发明包括以下步骤。

a.读入原始航空γ能谱数据，构建输入谱线数据矩阵A。

b.利用NASVD算法对输入谱线数据矩阵A进行计算，得到主成分矩阵S和幅度矩阵B，主成分矩阵S的各个列向量依次为各主成分，各主成分按照对应特征值的大小降序排列。

c.对主成分按照前1个、前2个、前3个…前n个进行分类。

d.分别对前1个、前2个、前3个…前n个主成分进行重构，分别得到重构谱矩阵A₁，A₂，A₃，…A_n。

e.依次计算相邻两个重构谱矩阵的相关系数|r_i|。

f.判断|r_i|的大小，当|r_i|第一次等于1，i值即被确定为有效奇异值的数目。

具体的，首先设原始航空γ能谱数据集中谱线条数为m，谱线道址数为n，得到m×n的输入谱线数据矩阵A。

得到输入谱线数据矩阵A后，利用NASVD算法求取主成分矩阵S和幅度矩阵B，求解步骤如下。

1）对输入谱线数据矩阵A各行元素求和，公式如下：

R=(r₁₁，r₂₁，… r_m1)^T

式中，R表示输入谱线数据矩阵A各行元素之和组成的列向量，r_i1(i=1，2，…，m)表示输入谱线数据矩阵A中第i行元素之和。

2）对输入谱线数据矩阵A各列元素求和，公式如下：

C=(c₁₁，c₁₂，…，c_1n)

式中，C表示输入谱线数据矩阵A中各列元素之和组成的行向量，c_1j (j=1，2，…，n)表示输入谱线数据矩阵A中第j列元素之和。

3）归一化单位谱，公式如下：

式中，N表示对C进行归一化处理后的行向量。

4）求取单位方差矩阵，公式如下：

式中，A_new表示m×n单位方差矩阵。

5）对A_new进行奇异值分解，公式如下：

A_new=UΣV^T

式中，U表示m×n正交矩阵，V表示n×n正交矩阵，Σ表示由奇异值λ组成的n×n对角矩阵，Σ=diag(λ₁，λ₂，…，λ_n)，且λ₁≥λ₂≥…λ_n≥0。

6）计算主成分矩阵和幅度矩阵，公式如下：

式中，S表示n×n主成分矩阵，其各个的列向量依次为各主成分，按照对应特征值的大小降序排列，diag(N)表示由向量N组成的对角矩阵，

同时计算幅度矩阵，公式如下：

式中，B表示m×n幅度矩阵，diag(R)表示由向量R组成的对角矩阵。

然后如图2所示，对得到的主成分按照“滑窗法”进行分类，即按照前1个、前2个、前3个…前n个进行分类，然后对主成分进行重构。

对主成分进行重构时，首先利用以下公式重构输入谱线数据矩阵A：

式中，b_i表示B中的第i个列向量，s_i ^T表示S^T中的第i个行向量。

按照主成分的分类，选取主成分矩阵S中代表信号的前k个（k=1，2，3，...，n）主成分，得到重构的航空γ能谱数据集A_k：

当k=1时，得到重构谱矩阵A₁；当k=2时，得到重构谱矩阵A₂，依次类推，得到重构谱矩阵A₁，A₂ ，A₃，…A_n。

对于得到的系列重构谱矩阵，利用相关系数计算公式，分别计算A_i和A_i+1的相关系数r_i(i=1，2，…，n-1)。

对计算得到的相关系数|r_i|的值进行判断，当|r_i|第一次等于1时，i的值便为有效奇异值的个数值。

如果相关系数|r_i|第一次等于1之后，后续|r_i|的值在0.99990~1.000000值之间变化，则航空γ能谱数据符合测量误差要求，可继续进行后续操作；否则，说明测量设备可能发生的能谱峰漂超出了峰漂允许的误差范围，需要对能谱数据进行峰漂修正或对航测设备进行检查。

确定了有效奇异值的数目i后，选取前i个主成分重构输入谱线数据矩阵A，得到降噪后的谱线数据，然后对降噪后的谱线数据进行后续处理和分析。

本发明首先对实测的原始航空γ能谱数据进行NASVD方法计算，获取实测谱线中的相互正交的谱线主成分，并按各自对谱线形状贡献的大小进行降序排列，接着利用“滑窗法”对获取的谱线主成分进行分类，对分类后的谱线主成分进行谱线重构，得到系列重构谱线。然后对系列重构谱线中相邻谱线之间的相关性进行分析计算，得到重构谱线的系列相关系数。以相关系数的绝对值作为判断依据，当相关系数的绝对值第一次等于1时，即可明确确定有效奇异值的数目i。最后采用前i个低阶主成分进行谱线重谱，最终得到降噪后的谱线。

主成分采用“滑窗法”进行分类，其基本思想主要依据的是NASVD方法的核心思想，即从一系列原始谱线数据中提取出相互正交的谱线主成分，各主成分以各自对谱线的形状贡献大小按降序排列，低序主成分体现谱线中绝大部分信号的谱形状，高序主成分体现谱线中不相关的噪声。依据信号相关而噪声不相关的原则，将代表信号的低阶主成分依次叠加后进行谱重构，求取其相邻重构谱的相关系数，随着不断叠加，相关系数不断增大，直至第一次等于1，由此可以得到叠加低阶主成分时的序号，即确定有效奇异值的数目。

在理想情况下，相关系数|r_i|的值随着i的增大逐渐增大直到等于1，且后续|r_i|的值稳定为1，当相关系数|r_i|的值第一次等于1时，表示前i个主成分对谱线形状的贡献最大，采用前i个低阶主成分进行谱线重谱即可，再引入i之后的主成分，就可能会引入噪声。

在非理想情况下，如出现峰漂或数据错误等，可能会导致相关系数|r_i|的值不会出现等于1的情况，但相关系数|r_i|的值随着i的增大是不断趋近于1的，这种情况下，选择相关系数|r_i|的值最接近于1时的i值。

如图3到图9所示，为了验证方法的有效性和适用性，分别对航空放射性测量模型标准装置上静态实测和野外项目实际飞行测量的航空γ能谱数据集进行了方法验证。静态实测和飞行实测使用的飞机为同一架Y-12型飞机，采用的测量系统为相同配置的3箱14条NAI(Tl)晶体探测器的航空γ能谱测量系统，采样时间为1s。每条谱线的道址数均为256道。用于方法验证的谱线数据道址范围为5~254道。

其中，静态实测数据集共3600条能谱数据。飞行实测数据集共980585条能谱数据。

图3为静态实测数据集相关系数计算结果，图6为飞行实测数据集相关系数计算结果。可以直观判断出静态实测数据集有效奇异值数目为4，飞行实测数据集有效奇异值数目为8。图9为静态实测数据集降噪前后的谱线比较，从谱线降噪效果来看，有效奇异值数目为4时，降噪后的谱线与平均谱线几乎重合，降噪效果良好。

图4为静态实测数据集的前6个主成分，图5为静态实测数据集对应的前50个奇异值曲线，根据图4和图5可以看出直接分析主成分或通过奇异值曲线进行判断，都可以得到静态实测数据集有效奇异值数目为4，与通过相关性系数确定得到的结果一致。

图7为飞行实测数据集的前11个主成分，图8为飞行实测数据集对应的前50个奇异值曲线，根据图7和图8可以看出直接分析主成分或通过奇异值曲线进行判断，都可以得到飞行实测数据集有效奇异值数目为8。与通过相关性系数确定得到的结果一致。

同时，由图6可知，当相关系数的值等于1之后，后续的相关系数|r_i|的值小于1且出现浮动，表示可能出现峰漂，但后续|r_i|的值在0.99990~1.000000值之间变化，表示航空γ能谱数据符合测量误差要求，可继续进行后续操作。

本实施例中用于举例的静态实测数据集和飞行实测数据集比较简单，可以通过人工进行判断有效奇异值数目，且人工判断结果和通过相关性系数确定得到的结果一致。对于实际工作中的航空γ能谱测量，得到的数据集要复杂的多，将无法通过人工判断，或不同的人判断得到的结果不一致，而使用本发明的有效奇异值数目确定方法，能够准确地确定有效奇异值数目。

经过多次实验研究表明，本发明方法可以准确、直观实现航空γ能谱NASVD降噪有效奇异值数目的选取，为实现航空γ能谱NASVD自动化降噪处理提供了理论依据。

Claims (5)

1.一种航空γ能谱数据NASVD降噪有效奇异值数目确定方法，其特征在于，包括以下步骤：

a.读入原始航空γ能谱数据，构建输入谱线数据矩阵A；

b.利用NASVD算法对输入谱线数据矩阵A进行计算，得到主成分矩阵S和幅度矩阵B，主成分矩阵S的各个列向量依次为各主成分，且各主成分按照对应特征值的大小降序排列；

c.对主成分按照前1个、前2个、前3个…前n个进行分类；

d.分别对前1个、前2个、前3个…前n个主成分进行重构，分别得到重构谱矩阵A₁，A₂ ，A₃，…A_n；

e.依次计算相邻两个重构谱矩阵的相关系数|r_i|；

f.判断|r_i|的大小，当|r_i|第一次等于1时，i值即被确定为有效奇异值的数目。

2.根据权利要求1所述的航空γ能谱数据NASVD降噪有效奇异值数目确定方法，其特征在于，在步骤b中，主成分矩阵S和幅度矩阵B的求解步骤如下：

1）对输入谱线数据矩阵A各行元素求和，公式如下：

R=(r₁₁，r₂₁，… r_m1)^T

式中，R表示输入谱线数据矩阵A各行元素之和组成的列向量，r_i1表示输入谱线数据矩阵A中第i行元素之和，其中i=1，2，…，m；

2）对输入谱线数据矩阵A各列元素求和，公式如下：

C=(c₁₁，c₁₂，…，c_1n)

式中，C表示输入谱线数据矩阵A中各列元素之和组成的行向量，c_1j表示输入谱线数据矩阵A中第j列元素之和，其中j=1，2，…，n；

3）归一化单位谱，公式如下：

式中，N表示对C进行归一化处理后的行向量；

4）求取单位方差矩阵，公式如下：

式中，A_new表示m×n单位方差矩阵；

5）对A_new进行奇异值分解，公式如下：

A_new=UΣV^T

式中，U表示m×n正交矩阵，V表示n×n正交矩阵，Σ表示由奇异值λ组成的n×n对角矩阵，Σ=diag(λ₁，λ₂，…，λ_n)，且λ₁≥λ₂≥…λ_n≥0；

6）计算主成分矩阵和幅度矩阵，公式如下：

式中，S表示n×n主成分矩阵，其各个的列向量依次为各主成分，按照对应特征值的大小降序排列，diag(N)表示由向量N组成的对角矩阵，

并计算幅度矩阵，公式如下：

式中，B表示m×n幅度矩阵，diag(R)表示由向量R组成的对角矩阵。

3.根据权利要求2所述的航空γ能谱数据NASVD降噪有效奇异值数目确定方法，其特征在于，在步骤d中，首先重构输入谱线数据矩阵A，公式如下：

式中，b_i表示B中的第i个列向量，s_i ^T表示S^T中的第i个行向量；

按照步骤c中的分类，选取主成分矩阵S中代表信号的前k个主成分，其中k=1，2，3，...，n，得到重构的航空γ能谱数据集A_k：

当k=1时，得到重构谱矩阵A₁；当k=2时，得到重构谱矩阵A₂，依次类推，得到重构谱矩阵A₁，A₂ ，A₃，…A_n。

4.根据权利要求1所述的航空γ能谱数据NASVD降噪有效奇异值数目确定方法，其特征在于，在步骤e中，利用相关系数计算公式，分别计算A_i和A_i+1的相关系数r_i，其中i=1，2，…，n-1。

5.根据权利要求1所述的航空γ能谱数据NASVD降噪有效奇异值数目确定方法，其特征在于，在步骤e中，如果相关系数|r_i|第一次等于1之后，后续|r_i|的值在0.99990~1.000000值之间变化，则航空γ能谱数据符合测量误差要求，可继续进行后续操作；否则，说明测量设备可能发生的能谱峰漂超出了峰漂允许的误差范围，需要对能谱数据进行峰漂修正或对航测设备进行检查。

Images