1、随机矩阵  stochastic matrix
定义：

　　A square matrix P is said to be stochastic if its elements are all nonnegative and all row sums are one.

　　如果一个方阵P的每一项均为非负的并且每行的和均为1，则称其为随机矩阵stochastic ；
2、双随机矩阵  doubly stochastic matrix
定义：

　　If in addition to being stochastic, all column sums are one, the matrix is said to be doubly stochastic. An obvious example of a doubly stochastic matrix is the n × n matrix in which each entry is 1/n.

　　如果在是stochastic矩阵的基础上，每列的和也均为1，那么此矩阵为 doubly stochastic。一个最简单的双随机矩阵例子是一个每项均为  1/n  的  n x n矩阵。
参考：website
3、bistochastic matrix
有些地方定义bistochastic matrix：the rows and columns of the matrix A sum to 1.（若矩阵A的行和与列和均为1，则A为bistochastic）

采用双门限法的思路：



首先通过噪声段（随机性大，混乱度高）的谱熵大于语音段，区分出语音的浊音与噪声，保留浊音，另一部分主要是发音起始和终止时的清音阶段与噪声，再由清音的短时过零率低于噪声段，将清音与噪声分离，保留清音，由此就得到了完整的语音段，实现了语音段的端点检测，可以自动去除噪声段。
谱熵定义：
按照矩阵运算的思路编辑程序：先计算分帧矩阵的fft结果，转为dB单位（实验测试了使用dB作为单位效果明显）平方后得到矩阵Y对Y每列的前N/2行相加得概率计算的分母cigmaY再把cigY拓展成N行，便于与Y做./的矩阵运算得到概率矩阵P代入最终的谱熵计算公式即可。代码如下：freq=fft(frame_w,N); %将数据由时域转换到频域spect=real(10log10(freq));%将频域结果单位换dB

Y=spectconj(spect);%计算能量sumY=sum(Y(1:end/2,:));sumY=sumY(ones(1,N),：);

P=Y./sumY;%计算样本点概率

H=-sum(P(1:end/2,:).*(log2(P(1:end/2,:))));%代入谱熵定义的公式
检测过程根据双门限法的思想，对分帧矩阵做一次帧数的遍历，假如某帧的谱熵小于了高谱熵门限，则可能进入了语音段，继续遍历，假如某帧的谱熵小于了低谱熵门限，确定已经进入了语音段，记录下该帧的位置，从该位置向前遍历，假如某帧的短时过零率大于了过零率门限，确定该帧为噪声首部结束的位置，记录改帧位置为noiseEnd,跳出向前遍历。  接着对噪声尾部处理的思路类似，假如某帧的谱熵小于低谱熵门限，而它的下一帧谱熵高于门限，确定该帧在语音片段，开始向后遍历，加入过零率大于了过零率门限，记录该帧位置为noiseBegin，跳出所有遍历。
思路还是非常清晰的，代码写起来不难。

先给一遍实验流程中得出的图：












下边上完整代码（仅供参考，实验中用到的语音数据为：oh.mat，你自己做可能需要修改相应的门限参数）

此外实验中需要的语音信号分帧函数，我有篇索引分帧法的博文里边有讲解和代码

传送门

%Ocross：短时过零率

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%输入参数：

% frame_w              分帧矩阵

%输出参数：

% zerocross       过零率序列

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

function [zerocross]=Ocross(frame_w)

[frame_length,frame_number]=size(frame_w);

zerocross=zeros(1,frame_number);

for i=1 : frame_number

u=frame_w(:,i);        % 取出一帧

for j=1 : frame_length-1

if u(j)* u(j+1)< 0       % 判断是否为过零点

zerocross(i)=zerocross(i)+1;   % 是过零点，记录1次

end % end过零判断

end     % end单帧循环

end         % end帧数循环
% Experiment 4 endpointDetection

% Apr.18 2020

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

clear

close all

clc
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%

load(‘oh.mat’,‘data’);

fs=10000;

N=1024;
width = 3;     % Width in inches

height = 3;    % Height in inches

alw = 0.75;    % AxesLineWidth

fsz = 13;      % Fontsize

lw = 1.2;      % LineWidth 1.5

msz = 7;       % MarkerSize
frame_time=20e-3;

[frame_m,frame_w,frame_length,frame_shift,frame_number]=enframe(data,fs,20e-3,10e-3,‘hamming’);

% 调用分帧加窗函数，获取分帧矩阵、分帧加窗矩阵、帧长、帧数

timeAxis=(1:frame_number)*20e-3;

zerocross=Ocross(frame_w);

plot(data)

legend(‘原始语音时域图’,‘Location’,‘best’)

xlabel(‘时间（采用点数）’);

ylabel(‘幅度’);
mag=sum(abs(frame_w));%幅度

f = (0:N/2-1)/Nfs;  %计算频率序列

freq=fft(frame_w,N); %将数据由时域转换到频域

spect=real(10log10(freq));%将频域结果单位换算为dB

Y=spect.conj(spect);%计算能量

sumY=sum(Y(1:end/2,：));

sumY=sumY(ones(1,N),：);

P=Y./sumY;%计算样本点概率

H=-sum(P(1:end/2,:).(log2(P(1:end/2,:))));%代入谱熵定义的公式

figure();

plot(timeAxis,H)

legend(‘谱熵时域图’,‘Location’,‘best’)

xlabel(‘时间/s’);

ylabel(‘幅度’);

figure()

subplot(211)

plot(timeAxis,mag)

subplot(212)

plot(timeAxis,H)

EntropyHigh=max(H)*0.995;%谱熵高门限

EntropyLow=min(H)*1.06;%谱熵低门限

hold on

plot([timeAxis(1), timeAxis(end)], [EntropyHigh, EntropyHigh], ‘r’, ‘LineWidth’,lw, ‘MarkerSize’, msz);

plot([timeAxis(1), timeAxis(end)], [EntropyLow, EntropyLow], ‘g’, ‘LineWidth’,lw, ‘MarkerSize’, msz);

legend(‘谱熵与赋值对比时域图’,‘高谱熵门限’,‘低谱熵门限’,‘Location’,‘best’)

xlabel(‘时间/s’);

ylabel(‘幅度’);

figure()

T=24;%设置过零率门限

plot(timeAxis,zerocross)

hold on

plot([timeAxis(1), timeAxis(end)], [T, T], ‘r’, ‘LineWidth’,lw, ‘MarkerSize’, msz);

legend(‘过零率’,‘门限’,‘Location’,‘best’)

xlabel(‘时间/s’);

ylabel(‘幅度’);
figure()

plot(data)

hold on

for i=1:frame_number

if(H(i)<EntropyLow && H(i-1)>EntropyLow )

noiseEnd=i-1;

for j=noiseEnd：-1:1

if(zerocross(j)>T)

noiseEnd=j;

break

end

end

end

if(H(i)<EntropyLow && H(i+1)>EntropyLow )

noiseBegin=i+1;

for j=noiseBegin:frame_number

if(zerocross(j)>T)

noiseBegin=j;

break

end

end

break

end

end

endIndex=noiseEnd*(frame_length+1)/2;

beginIndex=noiseBegin*(frame_length+1)/2;

noise1=data(1:endIndex);

noise2=data(beginIndex:length(data));

plot(noise1,‘r’)

hold on

plot(beginIndex:length(data),noise2,‘g’)

legend(‘原始语音’,‘首噪声段’,‘尾噪声段’,‘Location’,‘best’)

xlabel(‘时间（采样点数）’);

ylabel(‘幅度’);

figure()

oh_clean=data(endIndex:beginIndex);

plot(oh_clean)

legend(‘cleanData’,‘Location’,‘best’)

xlabel(‘时间（采样点数）’);

ylabel(‘幅度’);

save(‘oh_clean’)

一、观测值随机模型

用于描述观测值与未知参数之间关系的模型称为函数模型或数学模型；

而描述观测值本身统计特性的模型称为随机模型，主要通过一个适当的协方差矩阵来定义。

随机模型分类：详见参考3

等权随机模型
	
卫星高度角随机模型
	
信噪比随机模型
	
基于验后残差的随机模型
1、等权随机模型 ：

相对定位，认为观测值精度相等，双差中测量噪声，方差-协方差矩阵为：







另外一种表示方式：其中α^2为单差观测值方差，双差观测值的协方差矩阵如下：



2、高度角随机模型：

利用卫星高度角为变量的函数模型对观测量的方差进行估计：

其中Elev为高度角，函数/f通常有指数函数模型、正切函数模型、余弦函数模型、正弦函数模型；一般采用正弦函数模型。



为了更好地定义低高度角观测值地中误差，又常采用改进后地正弦函数模型：



式中，a、b都是经验值，一般取a=3/4mm，b=3mm。

3、载噪比随机模型：

（1）SIGMA-ε 随机模型



Ci的计算为：



式中，Bi为相位跟踪环带宽（Hz）；λ为载波相位波长（ｍ）。通常取Cl1=0.00224m^2Hz，Cl2=0.00077m^2Hz。

 

（2）SIGMA-Δ 随机模型



式中，Δ为模板值与实际观测值之间的差值；α为经验系数，一般取2。

4、信号强度随机模型：





二、GNSS系统观测值质量分析

      数据质量是GNSS卫星导航定位的精度和可靠性的重要保障，它受到外业观测环境、接收机自身以及卫星健康状况等诸多因素的影响。如何对GNSS原始观测值的质量作出较为准确的评价，从实测数据中剔除某些质量较差的数据成为高精度GNSS定位首先考虑且首要解决的问题，也是GNSS数据预处理的关键。

数据质量分析指标：详见参考5

多路径效应
	
电离层延迟
	
周跳比
	
信噪比
	
数据完整率
	
钟跳
1、伪距多路径效应分析：参考5、6、9

双频伪距、载波相位观测值：

以GPS的L1、L2频点伪距、载波相位为例：



其中，mp1、mp2分别表示L1、L2频点伪距和载波相位观测值的多路径效应组合；P是伪距观测值、λ是频点波长、Φ是载波相位观测值/cycle，α=(f1^2/f2^2)；

由于，载波相位中含有未知的模糊度参数/N，因此，上式求的多路径mp1、mp2，包含周跳信息；故在不发生周跳的情况下，可利用多个历元的数据取平均值，再将包含模糊度参数的序列减去该平均值，变得出2个伪距的多路径效应值，可用来衡量多路径效应影响的程度。

其中IGS数据质量检测分析显示，对于多路径效应而言，2/3的IGS站的mp1平均值小于0.5，而2/3的mp2平均值小于0.75。mp1、mp2越小，说明抗多路径效应能立越强。

2、电离层延迟及延迟变化率分析：参考7



 

3、周跳比分析：参考5

O/slips 是观测值和周跳的比，反映的是数据周跳的情况。

如果某颗卫星再单位历元的两个频率上，检测到周跳情况，则认为该历元产生周跳现象，TEQC软件采用电离层残差法进行周跳的探测与评定，运行的结果文件中以O/slips值来表示观测值与周跳比，O/slips值的参考值为200，O/slips值越小，说明出现周跳越严重。

还有另外一种形式，CSR：



IGS 的数据质量检测分析显示，超过半数的IGS站的CSR平均值小于5，2/3以上的CSR平均值小于10；

4、信噪比/SNR分析：

信噪比，定义为信号功率S与噪声功率N之比，可用来衡量测距信号质量的优劣，信噪比值越大，说明观测信号的质量越好。一般情况下，RINEX文件会直接给出信号的信噪比。

5、数据完整率分析：参考6

      数据完整性不仅仅指解析得到的观测历元数占整体应观测到的历元数的比例，还应该反映观测历元内出现的所有卫星的完整观测值，包括P1或者C/A码数据、P2或者C2码数据、L1和L2载波相位数据，W及L1和L2的SNR大于或等于指定的阐值等情况。因此不但要统计缺失历元数，还要在可解析的历元内对观测量的完整性进行统计和标记，并刪除不满足条件的观测数据。对于满足条件的观测数据计数器进行自加，则数据的时间完整性为满足条件的观测数据个数与当前历元和初始历元的差值的比值。



Have(i)为第i颗卫星的完整观测值数目，Expert(i)为第i颗卫星理论观测值数目。

6、钟跳分析：

钟跳与周跳十分类似，都会引起观测数据的跳变，但两者在原理和实质上是不同的，钟跳是接收机钟差的突然变化，可以引起所有卫星伪距或载波相位观测值的同时跳变。

不同的接收机有两种典型的钟跳：

频繁跳跃：在每个历元对接收机时钟进行修正，改正量一般较小，一般情况下可以不考虑；

毫秒跳跃：在接收机钟差达到一定数值时才进行修正。

三、观测值的线性组合

GNSS中，主要利用载波相位、伪距观测值；为了达到某种目的，而进行一系列的变形、组合：

同类型同频率观测值的线性组合：
比如：单差、双差等，目的是为了消除一些参数；

同类型不同频率观测值的线性组合：
比如：消除电离层延迟，电离层延迟分析

不同类型观测值的线性组合：
比如：伪距多路径分析、

1、同类型同频率观测值的线性组合

广为大家熟知的就是 单差、双差；

缺陷/热点：差分观测值间具有了相关性，如何设置随机模型/stochastic model ；对应此博客第一部分；

2、同类型不同频率观测值的线性组合

以GPS L1、L2载波相位为例：

给出线性组合观测值的频率、波长、整周模糊度、电离层延迟、组合观测值测量噪声关系：






或：






一般而言，对于用户有用的组合观测值应该满足以下标准：

线性组合后构成的新“观测值”应能保持模糊度的整数特性，以利于正确确定整周模糊度；
	
线性组合后构成的新“观测值”应具有适当的波长；
	
线性组合后构成的新“观测值”应不受或基本不受电离层折射的影响；
	
线性组合后构成的新“观测值”应具有较小的测量噪声；
（1）宽巷组合/Wide Lane

宽巷观测值Φ1-Φ2之差：n=1,m=-1;

利于求解模糊度，噪声为：0.01*sqrt(2)*(0.8619)=1.22厘米

波长为0.86m，模糊度为整数，适于中长基线的模糊度分解。首先用宽巷组合确定双频模糊度之差，然后引入LI组合或L1或L2观测方程，分解出L1模糊度。  



 

（2）无电离层折射延迟组合/ iono-free

具体见：GNSS中多频观测值的组合形式

消电离层组合消除了一阶电离层影响，但模糊度不再为整数了，且观测噪声比L1放大三倍，对中长基线解算有利，可显著改善中长基线解的精度。



（3）窄巷组合

波长为0.107m，模糊度为整数，宽巷组合和窄巷组合中电离层影响的大小相等，符号相反，适用于模糊度分解。

3、不同类型观测值的线性组合

用到了伪距和载波相位；

（1）不同类型的双频观测值间的线性组合之Melboune-Wubbena组合

Melboune-Wubbena组合 （双频伪距和相位组合），推导如下：

消除了电离层延迟、接收机钟差、卫星钟差、卫星至接收机的几何距离；仅受多路径和观测噪声的影响；



其中，ρ为卫星至接收机的距离与所有与频率无关的偏差改正项之和。


推导过程如下：



联合<2><3>代入<1>中，得到：




或：









（2）不同类型的双频观测值间的线性组合之电离层残差组合

在推导Melboune-Wubbena组合中，载波相位Φ1和Φ2单位为：周/cycle;将其乘以波长，得到米/m为单位，进行推导，得到电离层残差组合：

消除了电离层延迟、接收机钟差、卫星钟差、卫星至接收机的几何距离；仅受多路径和观测噪声的影响；

因此，Melboune-Wubbena 和 电离层残差组合，均可以进行模糊度分解和周跳探测；但是周跳探测时，不能对L1、L2频点同时发生的周跳进行检测。



（3）不同类型的单频观测值间的线性组合

由于伪距和载波观测值中电离层延迟大小相同，符号相反，故利用单频伪距P1与载波相位Φ1也能消除电离层延迟：



伪距和相位这种线性组合消除了一阶电离层影响，这种观测量除包含非色散性误差外，还包含一个模糊度参数，比较适用于单点定位，有利于改善定位的精度。

（4）附加：•Geometry-free又称电离层残差组合

这一组合与接收机至卫星的几何距离无关，消除了诸如轨道误差、接收机钟差、卫星误差和对流层误差，仅包含电离层及双频模糊度实数组合，适于电离层研究、数据编辑及周跳探测。



 

四、周跳探测  需要根据误差传播率计算误差

       每颗卫星的观测值在其可视的连续时间内应该是一条光滑的曲线，可通过观测值的线性组合构造合适的周跳检测量，通过求得周跳检测量的时间序列，采用合适的误差判别方法，判断出周跳发生的位置和大小，因此周跳探测与修复的思路如下：

利用观测值的线性组合构造合适的周跳检测量；
	
选择合适的误差判别方法；
	
探测周跳检测量序列中突变的位置，判断出是粗差还是周跳；
	
选用合适的方法进行周跳大小的修复；
周跳的本质是探测观测序列是否出现突变！ 

方法主要有：高次差法、多普勒观测值法、多项式拟合法、M-W组合法、电离层残差法等；

1、高次差法： 参考13，P37

对载波相位观测值，进行相邻历元间做差；

一次差：即卫星至接收机距离变化/速度：

二次差：速度变化/加速度；

三次差：加速度变化/加加速度； 

连续跟踪的载波相位值，随着次差的增加，结果趋于0。

其中三次差公式为：



需要注意的是：接收机晶振的短期稳定度、采样率；以此判定用高次差法是否合适。

2、多普勒观测值法

利用历元间载波相位差值 - 多普勒在时间段内积分值；

缺点：易受采样率、接收机钟差、电离层、噪声的影响，只能检测周跳比较大的情况；



3、M-W组合法 / 宽巷相位减窄巷伪距组合

宽巷模糊度



宽巷观测值的整周模糊度主要受伪距多路径和噪声影响，可通过多个历元的平滑进行削弱；

当有周跳发生时，宽巷模糊度为ΔNw=（N1+n1）- (N2+n2)，n1、n2为L1、L2载波上的周跳，其值将会出现较大的变化，在实际应用中，联合利用Nw及其均值对周跳的存在性进行判断。利用递推的方法计算出第i个历元的宽巷模糊度均值和方差为：









缺点：1、无法确定周跳发生在哪一个载波上；2、当两个载波上发生周跳大小相同时，将无法探测；

4、电离层残差法  参考14

同一历元双频载波相位测量差：



消除了接收机至卫星间的几何距离、接收机和卫星的钟差、对流层延迟、及部分电离层延迟；剩余误差为：频率间模糊度、电离层影响、多路径效应和观测噪声；

将上式两端同除以λ1，则有：





电离层残差检测方法特点：

只用到了双频载波相位观测量，无需其他信息；
	
只能判断出是否出现周跳，不能判断是哪一个频率出现周跳；
如果，电离层延迟变化缓慢，历元间电离层残差相减得到：



其中：77/60=1.28；即如果 |ΔΦ|>0.28，则说明在t+1历元 存在周跳，但是至于在哪一个频率发生周跳，则无法检测出来。

 

参考：

1、GPS/GLONASS/BDS/Galileo 系统载波相位观测值质量分析

2、北斗与GPS数据质量对比分析

3、基于观测值质量指标的GPS观测量随机模型分析

4、Stochastic Modeling for Static Gps Baseline Data Processing

5、Beidou/GPS/GLONASS多系统卫星定位数据质量比较分析

6、多系统GNSS实时数据质量分析及软件实现 朱静然

7、GNSS观测数据质量分析软件质量设计

8、Assessment of stochastic models for GPS Measurements with different types of receivers

9、GNSS数据质量分析

10、《GPS测量与数据处理》李征航

11、第九讲-观测值的线性组合.ppt

12、GPS观测值的线性组合解析.ppt

13、GNSS观测数据预处理及质量评估

14、基于电离层残差法的周跳探测和修复方法研究

总览图
线性结构数组矩阵
线性表
栈和队列
串

树
二叉树遍历
图








总览图



线性结构、数组、矩阵



线性表


  
定义：唯一的第一个和最后一个 

  存储分为顺序存储和链式存储：顺序存储的优点是随机存取，缺点是插入和删除需要移动大量的元素，空间需要事先分配。 

  链式存储又分为单链表和双链表，加入连接指针的概念。




栈和队列


  
栈，先进后出（last in first out）的特性，栈也分顺序存储和链式存储 

  常见的栈操作，1。创建，2。判断空，3。入栈，4。出栈，5。读取栈顶元素 

  队列，先进先出（first in first out）的特性，队列分顺序存储和链式存储 

  常见的队列操作：1。创建，2。判断空，3，入队，4。出队，5。读取队头元素




串


  
模式匹配KMP算法，关键设置next函数




树



二叉树遍历


  
先序，中序，后序是意思是访问根字节点的顺序，即先序先访问根节点，中序是先访问左再访问根最后是右，后序是左右中的顺序。




图


  
图的存储由邻接矩阵和邻接链表 

  邻接矩阵，有向图是非对称的， 而无向图则是对称的 

  带权的图 

  广度优先和深度优先


本文最后修改时间：2017年11月5日20:09:23（1）