一种基于加权窗求直射与反射信号SNR干涉峰值的方法与流程

一种基于加权窗求直射与反射信号snr干涉峰值的方法
技术领域
1.本发明涉及遥感领域，特别是涉及一种基于加权窗求直射与反射信号snr干涉峰值的方法。


背景技术：

2.利用gnss(global navigation satellites system)反射信号进行地球陆地物理参数的反演是遥感领域的新型技术之一，具有信号源广、设备便携，受大气影响小等优势。利用gnss反射信号对土壤湿度进行测量能够有效降低成本，提高测量覆盖面，将传统的单点测量转换到面的测量，减少了所需要的测量设备，使得测量数据更有效利用在农林、地下水循环等相关领域。
3.土壤湿度，又称为土壤含水量，是用来表示土壤干湿程度大小的物理量。土壤湿度在研究陆地能量循环和水循环过程中是一项非常重要的变量，除了可以影响净辐射能量向潜热和感热转换的比率,还可以影响大气降水如雨、雪等转变为蒸发、下渗以及地表径流的比例。尽管土壤含水量在全球水资源中只占到0.0005％，但是对农作物生长的重要性是不可忽略的。
4.在传统的卫星信号接收中，接收机收到的除了视距传播的直射信号外，还有经过多个路径反射后的信号，又称为多径反射信号。在导航定位应用中通常把多径反射信号当作干扰进行抑制，因为其会对定位准确性产生影响，并且难以完全消除。但实际上反射信号与反射面的物理特征如土壤湿度具有相关性。单一天线接收到的gnss信号包括了直射分量和反射的多径分量，由于地基情况下直射跟反射信号频率近似相等，只是传播路径长度不同，这两个分量会在接收机天线处发生较为稳定的干涉现象，形成干涉信号。在低高度角时，干涉振荡现象十分明显。测绘级接收机会将这种干涉信号以信噪比(snr,signal
‑
to
‑
noise ratio)的形式记录存储在设备记录中。干涉现象的变化决定于土壤表面的反射特性，因此，可以利用这种干涉现象反演反射面土壤湿度等一些重要物理属性。snr的数学模型可用下式表示:
[0005][0006]
只有反射信号的特性与反射面的物理特性有关，所以需要去掉干涉信号中的直射信号，得到反射信号的snr
m
。传统的方法没有考虑周期信号中噪声所产生的毛刺影响，往往会输出多个伪峰值。本发明根据归一化干涉幅度进行加权，设置一个比较窗，该窗长度单位为卫星仰角，对窗内的数据进行比较，得到该窗内的局部峰值，作为该窗的特征值，然后将所有同样单位窗的特征值作为新的序列，改变窗的大小重复上述步骤进行比较。


技术实现要素：

[0007]
本发明的目的是提供一种基于加权窗求直射与反射信号snr干涉峰值的方法，降低噪声对干涉信号所带来的的毛刺影响。
[0008]
为实现上述目的，本发明提供了如下方案：
[0009]
一种基于加权窗求直射与反射信号snr干涉峰值的方法，所述方法包括：
[0010]
获取接收机接收到的直射信号与反射信号的干涉信号；
[0011]
对所述干涉信号进行单周期提取；
[0012]
对单周期提取后的干涉信号进行归一化处理，得到归一化序列；
[0013]
基于所述归一化序列计算窗的长度；
[0014]
基于所述窗的长度计算每次迭代的峰值，得到峰值序列；
[0015]
检测所述峰值序列的个数是否为1，若为1，则表示全部搜索完成，输出此时的峰值序列，若不为1，则重复上述步骤直到峰值序列的个数为1。
[0016]
可选的，所述干涉信号的表达式为：
[0017][0018]
其中，a
d
表示直射信号辐射值，a
m
表示反射信号辐射值，表示直射信号与反射信号的相位差。
[0019]
可选的，对所述干涉信号进行单周期提取具体包括：
[0020]
计算所述干涉信号的均值；
[0021]
判断当前序列值是否大于均值；
[0022]
若否，序列左移，返回上一步骤；
[0023]
若是，开始信号截取序列左移，序列计数器加1；
[0024]
判断序列计数器是否大于周期t，若是，则输出截取序列；
[0025]
若否，则返回步骤“开始信号截取序列左移，序列计数器加1”。
[0026]
可选的，计算所述干涉信号的均值具体采用以下公式：
[0027][0028]
其中，θ表示gnss卫星高度角，e[]表示求均值运算，snr为接收机接收到的干涉信号的信噪比。
[0029]
可选的，对单周期提取后的干涉信号进行归一化处理，得到归一化序列具体采用以下公式：
[0030][0031]
其中，x(θ)表示干涉信号的半周期信号，x
p
(θ)表示归一化序列，表示干涉信号的均值。
[0032]
可选的，基于所述归一化序列计算窗的长度具体包括：
[0033]
给定窗的大小其中，[]为向上取整符号，n为归一化序列的长度；
[0034]
计算归一化的功率和γ，γ＝sum(x
p
(θ))；
[0035]
当前序列和l置零；
[0036]
l＝l 当前序列x；
[0037]
判断[l]＝[γ/k]是否成立；
[0038]
若不成立，则序列左移，序列个数λ 1，返回步骤“l＝l 当前序列x”；
[0039]
若成立，则判断窗个数n是否等于窗的大小k；
[0040]
若是，则输出序列λ；
[0041]
若否，则窗的个数n 1，返回步骤“当前序列和l置零”。
[0042]
可选的，基于所述窗的长度计算每次迭代的峰值，得到峰值序列具体包括：
[0043]
建立矩阵m；所述矩阵m为λ
j
×
k矩阵；
[0044]
对矩阵m的每一行求和，得到一个l维列向量；
[0045]
计算所述l维列向量最小值对应的下标；
[0046]
取所述矩阵m的前r列，构成新矩阵m
js
；
[0047]
计算所述新矩阵m
js
的每一列的标准差σ
i
；
[0048]
所有σ
i
＝0的元素对应的下标为原始序列的峰值，记为λ＝[λ1,λ2,
…
λ
k
]。
[0049]
可选的，所述矩阵m的表达式如下：
[0050][0051]
可选的，计算所述新矩阵m
js
的每一列的标准差σ
i
具体采用以下公式：
[0052][0053]
其中，r表示前r列，表示矩阵元素。
[0054]
根据本发明提供的具体实施例，本发明公开了以下技术效果：
[0055]
本发明中的上述方法可以有效去除噪声所带来的毛刺影响；
[0056]
可以提升每次迭代的计算效率，减少所消耗的时间；
[0057]
可以通过改变设计窗的个数来决定计算精度和计算效率，更加灵活方便。
附图说明
[0058]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
[0059]
图1为本发明实施例一种基于加权窗求直射与反射信号snr干涉峰值的方法流程图；
[0060]
图2为本发明实施例序列截取流程图；
[0061]
图3为本发明实施例窗大小计算流程图；
[0062]
图4为本发明实施例局部峰值计算流程图；
[0063]
图5为本发明实施例一种基于加权窗求直射与反射信号snr干涉峰值的系统结构示意图。
具体实施方式
[0064]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0065]
本发明的目的是提供一种基于加权窗求直射与反射信号snr干涉峰值的方法，降低噪声对干涉信号所带来的的毛刺影响。
[0066]
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
[0067]
图1为本发明实施例一种基于加权窗求直射与反射信号snr干涉峰值的方法流程图，如图1所示，所述方法包括：
[0068]
步骤101：获取接收机接收到的直射信号与反射信号的干涉信号。
[0069]
所述干涉信号的表达式如下：
[0070][0071]
其中，a
d
表示直射信号辐射值，a
m
表示反射信号辐射值，表示直射信号与反射信号的相位差。
[0072]
步骤102：对所述干涉信号进行单周期提取。
[0073]
具体流程图如图2所示，步骤如下：
[0074]
计算所述干涉信号的均值；
[0075]
判断当前序列值是否大于均值；其中，当前序列值指的t1时刻的值。
[0076]
若否，序列左移，返回上一步骤；其中，序列左移相当于时间右移，即处理下一个时刻的值。
[0077]
若是，开始信号截取序列左移，序列计数器加1；从第一个均值点开始，截取一个周期t的长度，即为一个单周期。
[0078]
判断序列计数器是否大于周期t，若是，则输出截取序列；
[0079]
若否，则返回步骤“开始信号截取序列左移，序列计数器加1”。
[0080]
步骤103：对单周期提取后的干涉信号进行归一化处理，得到归一化序列。
[0081]
具体的，是用均值除以干涉信号的半周期信号，得到归一化序列，公式如下：
[0082][0083]
其中，x(θ)表示干涉信号的半周期信号，x
p
(θ)表示归一化序列，表示干涉信号的均值。
[0084]
步骤104：基于所述归一化序列计算窗的长度，具体流程如图3所示。
[0085]
首先给定窗的大小其中，[]为向上取整符号，n为归一化序列的长度；在实际应用中可以对k进行调整，窗的个数越多，计算越精确，相反的计算时间越长。
[0086]
记上一步骤中的归一化的功率和为γ＝sum(x
p
(θ))，则第一个窗的大小表示为：l1＝x1 x2
…
x
i
，其包含的序列个数为：λ1＝i
‑
1 1，使得[l1]＝[γ/k]，同理可求得第二个窗
的大小：l2＝x
i
x
i 1

…
x
s
，其包含的序列个数为：λ2＝s
‑
i 1，使得[l2]＝[γ/k]，记所求的所有窗包含的序列个数为：λ＝[λ1,λ2,
…
,λ
k
]，具体步骤如下：
[0087]
首先给定窗的大小其中，[]为向上取整符号，n为归一化序列的长度；
[0088]
计算归一化的功率和γ，γ＝sum(x
p
(θ))；
[0089]
当前序列和l置零；
[0090]
l＝l 当前序列x；
[0091]
判断[l]＝[γ/k]是否成立；
[0092]
若不成立，则序列左移，序列个数λ 1，返回步骤“l＝l 当前序列x”；
[0093]
若成立，则判断窗个数n是否等于窗的大小k；
[0094]
若是，则输出序列λ；
[0095]
若否，则窗的个数n 1，返回步骤“当前序列和l置零”；
[0096]
步骤105：基于所述窗的长度计算每次迭代的峰值，得到峰值序列，具体流程如图4所示。
[0097]
得到所有窗包含的序列个数后，接着建立一个λ
j
×
k矩阵m；
[0098][0099]
其中，矩阵中每一个元素的定义如下：
[0100][0101]
对矩阵的每一行求和，得到一个l维列向量：其中，γ
kj
表示行求和值，λ
j
表示矩阵行数，m
i,k
表示矩阵元素；
[0102]
计算所述l维列向量最小值对应的下标：μ
j
＝argmin(γ
kj
)；
[0103]
取所述矩阵m的前r列，构成新矩阵m
js
；
[0104]
计算所述新矩阵m
js
的每一列的标准差σ
i
：其中，r表示前r列，表示矩阵元素；
[0105]
所有σ
i
＝0的元素对应的下标为原始序列的峰值，记为λ＝[λ1,λ2,
…
λ
k
]。
[0106]
步骤106：检测所述峰值序列的个数是否为1，若为1，则表示全部搜索完成，输出此时的峰值序列，若不为1，则重复上述步骤直到峰值序列的个数为1。
[0107]
图5为本发明实施例一种基于加权窗求直射与反射信号snr干涉峰值的系统示意
图，如图5所示，该系统包括：
[0108]
单周期提取模块、归一化功率处理模块、窗设计模块、局部峰值求解模块、迭代控制模块。首先需要给定一个窗的个数k，其个数决定了计算精度和计算效率。干涉信号首先通过单周期提取模块，产生一个周期内的干涉信号；然后通过归一化功率处理模块，得到基于单周期均值的归一化功率；接着窗设计模块根据前置模块所得到的的归一化功率来设计窗的大小，归一化功率越小的位置，所设计的窗大小越大，窗大小与归一化功率成反比；根据前置模块所设计的窗，局部峰值求解模块对窗内的序列利用二分法进行峰值的寻找，将寻找到的峰值存储起来，记为该窗的1阶特征值；迭代控制模块在收到所有窗内序列一次迭代全部计算完成的信号后，将该次迭代的全部特征值合并成一个新的序列，重复进行以上的工作，直到窗的加权值为归一化的1，即窗已经覆盖了全部的序列，此时该窗的特征值即为峰值。
[0109]
所述的归一化处理模块首先计算单个周期信号的均值，然后对单个周期内的全部周期信号进行归一化，并存储。
[0110]
所述的窗设计模块根据窗的个数k与归一化的功率值设计窗的大小，窗的个数为k个，窗的大小与归一化的功率值成反比。
[0111]
所述的局部峰值求解模块在单个窗内进行峰值寻找，在单个窗内利用二分法不断逼近峰值，由于所涉及的信号为正弦函数，所以设计二分法的加权系数比值为正切函数的比值，即为1：tan(1/2)，对单个窗函数所求解的峰值进行存储，记为n阶特征值，n为迭代的次数。
[0112]
所述的迭代控制模块控制前述三个模块，当局部峰值求解模块的所有窗全部计算完成后，向迭代控制模块发送迭代完成的信号，收到此信号后，迭代控制模块将局部分值求解模块的全部n阶特征值记为新的序列，然后控制局部峰值求解模块进行n 1次迭代。若第n次的局部峰值求解模块窗的加权系数为1，即已经覆盖到了全部序列信号后，迭代结束，迭代控制模块输出此时的n阶特征值。
[0113]
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
[0114]
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。