第34卷第6期 系统工程与电子技术 Vol.34 No.6 2012年6月 Systems Engineering and Electronics June 2012 文章编号：1001506X(2012)06124105 大气低频噪声混合模型的MCMC参数估计 应文威，蒋宇中，刘月亮 (海军工程大学电子工程学院，湖北武汉430033) 摘要：大气噪声是低频通信中的主要干扰，且其有严重非高斯分有特性，对非高斯噪声模型的参数估计对 于提高低频接枚机的性能具有重要意义。设计了估计非高斯混合模型参数的马尔可夫链蒙特卡罗(Markov chain Monte Carlo,MCMC)算法，该算法通过构建见叶斯层大模型，利用Gibbs.袖样和M-H抽样更新选代参数。利用 乘积特性，将稳定分布作为等价的高斯分布来处理，并在层次模型中设置多个额外参数，以增强其灵活性。仿真 实验与实测数摄表明，这算法选代收敛快、精度高，有很高的实用价值。 关健词：混合模型；马尔可夫链蒙特卡罗：非高斯噪声；口稳定分布 中图分类号：TN911.23 文献标志码：AD01:10.3969/1.issn.1001-506X.2012.06.29 Parameter estimation for mixture model of atmospheric noise through MCMC method YING Wen-wei,JIANG Yu-zhong,LIU Yue-liang (College of Electronic Engineering,Naval University of Engineering,Wuhan 430033,China) Abstract:Atmospheric noise is the main interference in a low-frequency communication system,which is highly impulsive.So the work for estimating the parameters of model of non-Gaussian noises is of great signifi- cance to improve the performance of the low-frequency receiver.This paper proposes a Markov chain Monte Carlo (MCMC)method to estimate the parameters of a mixture model.The method updates the parameters through a Gibbs sampler and M-H algorithm,which are based on the Bayesian hierarchical model.The a stable distribution in the mixture model is equivalent to the normal distribution by using the product properties.An ex- tra layer is added to the hierarchy for full flexibility.The result shows that the new method has a good perform- ance,high precision and can be excellently applied in practice. Keywords:mixture model;Markov chain Monte Carlo (MCMC);non-Gaussian noise;a stable distribution 0 引言 注。对该混合模型的参数估计，传统上采用最大期望(x pectation-maximization,EM)算法或最大似然估计)。但 低频通信系统中的最优接收机一直以来都是基于高斯 由于该混合棋型参数较多，且对称α稳定分布的概率密度 噪声的假设，然而这类接收机在实际应用，特别是在低频段 设有闭合表达式，计算难度大，精度也不够理想。 通信系统中，往往达不到预期的性能。原因在于低领段受 马尔可夫链蒙特卡罗(Markov chain Monte Carlo, 雷电等形响，其噪声特性往往具有非高斯特征)，概率密 MCMC)算法，是一类在统计计算中有优异性能的方法。本 度函数在尾部较高斯分布收敛更慢，具有重尾的特点。为 文通过设计MCMC多层次算法来估计混合模型参数，迭代 了使接收机能够进行最优接收，则必须对噪声模型进行参 收敛快、精度高、具有明显优势。 数估计，非高斯分布和高斯分布的混合模型对自然界的非 高斯噪声有很强的代表性，有极为广泛的应用。本文研 1a稳定分布 究的是对称a稳定(symmetric a stable.,SaS)分布和高斯分 1.1定义 布的混合模型】。该模型涵盖了在非高斯模型应用广 如果一个随机变量服从稳定分布，则其特征函数满足： 泛的稳定分布和混合高斯分布)，因而受到学者的广泛关 (t)exp (j-1 nt '[1+isgn (t)w(t.a)]}(1) 收稿日期：2011-09-264修回日期：2012-03-03. 基金项日：国防预研基金资助课题 作者葡介：应文成(1987-)，男，博土研究生，主要研究方向为大气低频噪声研究及相关通信技术。E-mail:yingwenwei@sina.cm 万方数据
第34卷第6期 系统T程与电子技术 V01．34 No．6 2012年6月 Systems Engineering and Electronics June 2012 文章编号：100l 506X(2u12)06 1241 05 大气低频噪声混合 模型的MCMC参数估计 应文威，蒋宇中，刘月亮 (海军工程大学电子工程学院，湖北武汉430033) 摘要：大气噪声是低频通信中的主要干扰，且具有严重非高斯分布特性，对非高斯噪声模型的参数估计对 于提高低频接收机的性能具有重要意义。设计了估计非高斯混合模型参数的马尔可夫链蒙特卡罗(Markov chain MonteCarlo，MCMC)算法，该算法通过构建贝叶斯层次模型。利用Gibbs抽样和M—H抽样更新迭代参数。利用 乘积特性，将稳定分布作为等价的高斯分布来处理。并在层次模型中设置多个额外参数，以增强其灵活性。仿真 实验与实测数据表明，该算法迭代收敛快、精度高。有很高的实用价值。 关键词：混合模型；马尔可夫链蒙特卡罗；非高斯噪声；n稳定分布 中图分类号：TN 911．23 文献标志码：A DOI：10．3969／j．issn．1001-506x．2012．06．29 Parameter estimation for mixture model of atmospheric noise through MCMC method YING Wen—wei，JIANG Yu—zhong，LIU Yue--liang (College of Electronic Engineering，Naval University of Engineering，Wuhan 430033，China) Abstract：Atmospheric noise is the main interference in a low-frequency communication system，which is highly impulsive．S0 the work for estimating the parameters of model of non-Gaussian noises is of great signifi· cance tO improve the performance of the low—frequency receiver．This paper proposes a Markov chain Monte Carlo(MCMC)method tO estimate the parameters of a mixture model．The method uNates the parameters through a Gibbs sampler and M—H algorithm，which are based on the Bayesian hierarchical model．TheⅡstable distribution in the mixture model is equivalent tO the normal distribution by using the product properties．An ex— tra layer is added tO the hierarchy for full flexibility．The result shows that the newmethod has a good perform— ante，high precision and can be excellently applied in practice． Keywords：mixture rnodd；Markov chain Monte Carlo(MCMC)；non-Gausman noise；a stable distribution 0引 言 低频通信系统中的最优接收机一直以来都是基于高斯 噪声的假设，然而这类接收机在实际应用，特别是在低频段 通信系统中，往往达不到预期的性能。原因在于低频段受 雷电等影响，其噪声特性往往具有非高斯特征[1…，概率密 度函数在尾部较高斯分布收敛更慢。具有重尾的特点。为 了使接收机能够进行最优接收，则必须对噪声模型进行参 数估计。非高斯分布和高斯分布的混合模型对自然界的非 高斯噪声有很强的代表性，有极为广泛的应用[I．7]。本文研 究的是对称口稳定(symmetric口stable，泌)分布和高斯分 布的混合模型哺】。该模型涵盖了在非高斯模型应用广 泛的n稳定分布和混合高斯分布¨]，因而受到学者的广泛关 注。对该混合模型的参数估计，传统上采用最大期望(ex— pectation-maximization，EM)算法或最大似然估计Ⅲ。但 由于该混合模型参数较多，且对称a稳定分布的概率密度 没有闭合表达式。计算难度大，精度也不够理想。 马尔可夫链蒙特卡罗(Markov chain Monte Carlo， MCMC)算法，是一类在统计计算中有优异性能的方法．本 文通过设计MCMC多层次算法来估计混合模型参数，迭代 收敛快、精度高、具有明显优势。 1口稳定分布 1．1定义 如果一个随机变量服从稳定分布，则其特征函数满足： 9(￡)=exp{jot—l弦l。[1+j肛gn(￡)∞(￡．口)]}(1) 收稿日期：201l—09—26；修回日期：2012一03一03。 基金项目：国防预研基金资勘课题 作者简介：应文威(1987一)，男，博士研究生，主要研究方向为大气低频噪声研究及相关通信技术。E-mail：yingwenwei@sina．err 万方数据
·1242· 系统工程与电子技术 第34卷 因此，4可以看成是独立变量，满足： -tan(g),a≠1 P(名=j》=地，j=1,2 (7) (t,a）= (2) 21nt,。=1 2.2贝叶斯层次模型和先验 贝叶斯推断通过先验分布来推嘶后验分布 易看出，稳定分布由4个参量决定：a是特征指数， P(A|B)oc P(B A)P(A) (8) a的值越小，远离中心值的出现摄率越大；9是对称参数，当 根据贝叶斯层次理论，参变址集{w,a,Y,d}的联合概 一0变成对称a稳定分布：Y是尺度参数，代表了杂散程 率密度为 度，同能敏密切相关：d是位移参数。对于SS来说，当10,若X~f(x: 对于指数参数a,选取先验使得。的值在(0,2]上均匀 a1,0,0,Y,Y~f(xa,1,0,(cos(aa/2)). 分布，为 则X六仍服从a稳定分布，且满足X片~f(xa1a, pala)=是-tae0,2】 (12) 0,0,y)。 由于a稳定分布具有SMiN(scale mixture of normal) 参量w的共轭先验选取为Dirichlet分布： 的性质2)，则a稳定分布可以表示为X=Y6,Y为正态分 D(8,6》 (13) 布，g为正的独立随机变量。又因为当a=2时，a稳定分布 使超参数B,K具有更强的灵活性，多增加一个额外层。它 即为高斯分布。因此综合上述有：如果x是对称α稳定分 们的共轭先验为 (14) 布的独立同分布抽样，满足x4一f(a,0,0,Y）,则可得下列的 B~I(g.h)~T(r,e) 混合模型的直接非循环图(directed acyclic graph, 等价式子： DAG)如图1所示。为表述简单，设9={a,Y,}和6的超 ~N(0,Ay) (4) 参数n={axe,a,3,r,e,g,h},则式(8)可扩蔑为 ~f（xra/2,1,0,2(cos(a/4)片) .(5) P(w,9,7z,8,x)=P(P(8)P(x|8,z)X 2混合模型 P(z|w)P(w18》P(8| (15) 2.1混合模型 对称α稳定分布和高斯分布的混合模型如下： p(x）=fxa,0,0,y》）+g(x0,d) (6) 式中，f(xa,0,0,y)为非高斯SaS部分g(x0,d)为高斯 部分：w,和为权重因子。显然地，满足两十=1。因 此，需要估计的参数集为{w,a,Y,o2}。 假设观测到的数据集为X={工：|1≤k≤N),引人标签 变量Z={1≤≤N,其中名的取值为： (1)当x属于对称a稳定分布部分的时候，x4的取值 为1； 0:参：□：定值。 (2)当x属于高斯分布部分的时候，4的取值为2。 图】混合模型的DAG 万方数据
·1242· 系统下程与电了技术 第34卷 w(t．口)== 因此，缸可以看成是独立变量，满足： ，口≠1 ⅢP(zI=J)=m，j=1，2 (7) (2) ‘“ 一”。 … … ． 2．2贝叶斯层次模型和先验 ， 口2 l 易看出。a稳定分布由4个参量决定【B]：a是特征指数。 a的值越小，远离中心值的出现概率越大；口是对称参数，当 口=0变成对称a稳定分布；，，是尺度参数，代表了杂散程 度，同能量密切相关{8是位移参数。对于&S来说，当10，若x～，(z； a1，0，0，y)，y～，(z；毗，1，0，(COS(m2／2))i)。 则x讨仍服从n稳定分布．且满足x晴～，(z㈣却。 0·0，y)。 由于口稳定分布具有SMiN($ca|e mixture of normal) 的性膨“““，则n稳定分布可以表示为X=110，Y为正态分 布，口为正的独立随机变量。又因为当a=2时，a稳定分布 即为高斯分布。因此综合上述有：如果也是对称a稳定分 布的独立同分布抽样，满足缸～f(a，0，0，y)，则可得下列的 等价式子： 孔～N(0，丸y2) (4) k～f(x；也／Z。1，O。2(cos(腑2／4))专) ．(5) 2 混合模型 2．1混合模型 对称a稳定分布和高斯分布的混合模型如下： p(z)="Wlf(x；a，0．0t】，)+t地g(z；0，，) (6) 式中，f(x}a，0，0。’，)为非高斯S口S部分；g(x}O，矿)为高斯 部分；w．和毗为权重因子。显然地，满足硼。+妣=l。因 此，需要估计的参数集为{w，a，y，引。 假设观测到的数据集为X一{z。l 1≤^≤N}，引入标签 变量z={≈Il≤女≤N}，其中≈的取值为： (1)当以属于对称a稳定分布部分的时候，≈的取值 为1； (2)当而属于高斯分布部分的时候，≈的取值为2。 贝叶斯推断通过先验分布来推断后验分布 P(A B)。c P(B A)P(A) (8) 根据贝叶斯层次理论，参变量集{_lI，，a，y，矿}的联合概 率密度为 p(a，y，，，w，z，x)=p(x{口，y，，，w，z)× 户(口·y-，，w，z) (9) 选取合理的共轭先验能够充分利用先验信息，提高贝 叶斯模型的性能。选取以下先验： 由于高斯分布的参数具有共轭先验，因此对于式(4)中 高斯部分的参数一有 r2～F(ao．p (10) 式中，r(·)为Gamma分布；％，口分别为形状参数和反尺 度参数。 利用第1．2节中的结论，则对称a稳定分布参数y2有 厂2～r(￡，Ⅳ) (11) 对于指数参数a。选取先验使得。的值在(o，2]上均匀 分布，为 P(口I n)：_1：百1，口∈(0，zl (12) 参量w的共轭先验选取为Dirichlet分布： w～D(d，占) (13) 使超参数口，*具有更强的灵活性，多增加一个额外层。它 们的共轭先验为 口～r(g，^)，Ⅳ～r(r，e) (14) 混合模型的直接非循环图(directed aeyclic graph， DAG)如图l所示。为表述简单，设日：{口，y，矿)和0的超 参数_={a，x，E，a，卢，r，P，g，h}，则式(8)可扩展为 P(w，0·n·z·占，工)=P(7／)P(d)P(x 1日，z)× ．P(z l’I，)P(_I．．{8)P(O l") (15) O：参量；口：定值。 图1混合模型的DAG ) 坚0 o m( h 1 ，●●●●●●●，、●●●●●【 呈。 万方数据
第6期 应文威等：大气低频噪声混合模型的MCMC参数估计 ·1243· 3MCMC算法 P(r:I w.Y,,a) (20) 在本节将设计具体的MCMC算法来估计混合模型的 P(I.I wYo,d) 参数。在MCMC算法中，Gibbs抽样和Metropolis-Has tings(MH)算法是两种常用的计算方法。一般来说，若是 若不接收新值，则。=。。 3.5通过Gibbs抽样更新参数K,B 后验概率较易计算则常采用Gibbs.,不然采用M-H较多。 分别对：，日的全条件分布进行Gibbs采样来更新参数 对于本文，采用Gibbs和M-H的混合算法，算法流程如下， 的值，它们的全条件分布分别为 步骤1通过Gibbs抽样更新权重系数w: e|…~(r+ae+y2) (21) 步骤2通过Gibbs抽样更新参数Y; 月|…~r(g+a,h+g) (22) 步骤3通过Gibbs抽样更新参数子， 3.6更新标签变量z 步骤4通过M-H算法更新a: 变量指示了观测数据x,所属的概率分布，对于其全 步骤5通过Gibbs抽样更新参数(e,), 条件分布满足 步骤6更新标签变量z: wif(I.l a,y) 步骤7通过M-H算法更新变量A: P(x,=1|）=fxa,)+gx1可 3.1通过Gibbs抽样更新权重系数 (23) 直接计算w后验概率，易得”的全条件后验分布仍服 wf(z.a) Pz=2引)“fx,1a,刀+ga1o 从Dirichlet分布，即 (24) w|…~D(6十m1,8+mz) (16) 通过生成U(0,1)的均匀随机数，确定其所在区间来判 式中，聊…表示在除焊外的所有变量的条件下：n:代表属 定的取值。 于对称α稳定分布部分的观测数据数目：代表属于高斯 3.7M-H算法更新变量1 部分的观测数据数目。通过生成Dirichlet分布随机数更新 由于变量入的引入，使得对称a稳定分布可以写成入： 的值。 条件下的高斯分布。由于观测数据只有一部分属于对称α 3.2通过Gibbs抽样更新参数Y 分布，因此仅更新那一部分对应的值。易得变址A,的全条 通过参数Y的全条件分布更新参数Y的值，易算得 件分布为 Y2的全条件分布为 P(A;a,z=1)oc N(z.10.A.y)X y1~r+0.5m0.52三+ f(2,a/2,1,0,2(cos(元a/4)) (25) (17) 采用MH算法更新λ，的值，进用建议分布为 3,3通过Gibbs抽样更新参数d A~faa/2,1,0,2(cos(xa/4)) (26) 通过参数。2的全条件分布更新参数。的值，的全 则易得接受概率为 条件分布为 R=min N(z.10,) (27) g1…ra+0.5%0.52+p (18) N(,10,y) ,=1=g 生成U(0,1)的均匀随机数，若小于接受微率则接受新 3.4通过M-H算法更新a 值，否则不变， 对于参数。，由于其全条件据率分布不是常用的概率分 布，因此采用M-H算法更新a的值。 4仿真及实测结果 (I)从第t次选代时，从建议分布q(。+|d)产生一个 在仿真实验中，按照预设的权重生成混合模型。其中α 新值+。在本文中，选用正态分布。即a+~N(a,0.3), 稳定分布随机值生成采用常用的Chambers-Mallows-Stuck (2)以接受概率min(1,A)接收新值a+1,其中A为 算法。对于预设参数，设置为r=0.2,c=1,g=0.2,h=1, e=2,a。=2。设置预烧期为N。=200,总的送代次数为N= ⅡPgl,yo,a)XP(%y,e,e“)g(d|d) A= 1000。在实验中，将混合模型参数设为设置为a=1.5,w= II P(wYa,d)x P(w Y)q(a"Id) (0.7,0.3),y=7,d=4。仿真数据数目设为2000。仿真 结果如图2一图5所示。其中，心的送代收敛情况由于限 (19) 于篇幅并没有给出。但因为，=1一心，因此不彰响读者 由于α是独立的且先验服从均匀分布，建议分布是对 理解。将各个参量在预烧期之后的选代值作平均，得到结 称的，于是式(19)可化为 果为a=1.52,=0.71,y=7.1d2=3.9。 万方数据
第6期 廊文威等：大气低频噪声混合模型的MCMC参数估计 ·1243。 3 MCMC算法 在本节将设计具体的MCMC算法来估计混合模型的 参数。在MCMC算法中，Gibbs抽样和Metropolis-Hastings(M-H)算法是两种常用的计算方法。一般来说，若是 后验概率较易计算则常采用Gibbs，不然采用M—H较多。 对于本文，采用Gibbs和M-H的混合算法，算法流程如下： 步骤1通过Gibbs抽样更新权重系数wI 步骤2通过Gibbs抽样更新参数y； 步骤3通过Gibbs抽样更新参数矿； 步骤4通过M—H算法更新af 步骤5通过Gibbs抽样更新参数(K，p I 步骤6更新标签变量z， 步骤7通过M-H算法更新变量A； 3．1通过Gibbs抽样更新权重系数w 直接计算w后验概率，易得Ⅳ的全条件后验分布仍服 从Diriehlet分布。即 w I…～D(占+n1，占+n2) (16) 式中，wl¨·表示在除W外的所有变量的条件下；”，代表属 于对称a稳定分布部分的观测数据数目；”：代表属于高斯 部分的观测数据数目。通过生成Diriehlet分布随机数更新 Ⅳ的值。 3．2通过Gibbs抽样更新参数r 通过参数y1的全条件分布更新参数y的值，易算得 y_2的全条件分布为 川…～I’(e十o．5吣0．5。蚤．詈h)(17) 3．3通过Gibbs抽样更新参数一 通过参数一-2的全条件分布更新参数，的值，口-2的全 条件分布为 ，2|．…rf‰+0．5n：，0．5∑z；+∞ (18) 3．4通过M-H算法更新口 对于参数a．由于其全条件概率分布不是常用的概率分 布，因此采用M-H算法更新a的值。 (1)从第t次迭代时，从建议分布q(∥“10)产生一个 新值一“。在本文中，选用正态分布。即口“1～N(口‘，0．3)。 (2)以接受概率min(1，A)接收新值口f+1，其中A为 ⅡP(xl l岷，nd，∥1)XP(岷，y洒矿1)口(一I，1) A一三尘≠L———————————————————————一 ⅡP(xl I址．，y，d，d)x P(岷，y，口，矿1)q(∥1 0) I。1”1 (】9) 由于a是独立的且先验服从均匀分布，建议分布是对 称的，于是式(19)可化为 A= ⅡP(x (20) ⅡP(x。I K，ym口‘) t。l—i‘1 若不接收新值。则n“1=a‘。 3．5通过Gibbs抽样更新参数K，一 分别对￡，口的全条件分布进行Gibbs采样来更新参数 的值，它们的全条件分布分别为 Ⅳ|．…r(r+口，e+r2) (21) 口|．…F(g+口，h+a-2) (22) 3．6更新标签变量z 变量≈指示了观测数据z。所属的概率分布。对于其全 条件分布满足 盹，=，|¨归可醋筹酵辆 (23) P(≈一z I…)=石了ii彳等焉宇芊‰ (24) 通过生成U(O，1)的均匀随机数，确定其所在区间来判 定z。的取值。 3．7 M-H算法更新变量A 由于变量A。的引入，使得对称a稳定分布可以写成A； 条件下的高斯分布。由于观测数据只有一部分属于对称a 分布。因此仅更新那一部分对应的值。易得变量A，的全条 件分布为 P(^，I口，墨，，，，≈=1)oc N(x。l 0，A；y2)X ，(A。；口／2，1，0，2(COS(m2／4))“4) (25) 采用M-H算法更新A。的值，选用建议分布为 时1～fCh；a／2，1，0，2(cos(m／4))”。) (26) 则易得接受概率为 r R—min 1 【 N(z。I O，A¨矿) ’雨瓦T丽 (27) 生成U(O，1)的均匀随机数，若小于接受概率则接受新 值，否则不变。 4仿真及实测结果 在仿真实验中。按照预设的权重生成混合模型。其中a 稳定分布随机值生成采用常用的Chambers-Mallows—Stuck 算法。对于预设参数，设置为r=0．2。e=1，g=0．2。h=1· t=2’口0=2。设置预烧期为N。=200，总的迭代次数为N= l 000。在实验中，将混合模型参数设为设置为口=1．5，w= (o．7，0．3)，)，=7，矿=4。仿真数据数目设为2 000。仿真 结果如图2～图5所示。其中，锄的迭代收敛情况由于限 于篇幅并没有给出。但因为劬=1一"，，因此不影响读者 理解。将各个参量在预烧期之后的迭代值作平均，得到结 果为口；1．52，∞l=O．71。y=7．1，一一3．g。 万方数据
·1244· 系统下程与电子技术 第34卷 5 在图2一图5中，容易发现在达代200次后，各个参量 0 都能趋于收敛，总体上收敛速度快，结果精度较高。在实验 50 中，a的更新过程需要计算。稳定分布的概率帝度值，以往 一般采用Nolan提出的数值积分进行计算，但耗时长，计算 过大。在本文中，通过对特征函数做快速傅里叶变换，求特 20 定序列的概率密度值，再插值求各个点的概率密度值，这样 可以大幅降低时间消耗。 10 困6为实测大气噪声数据，地点为河南某地，采样援率 0 为40960Hz,带宽80Hz,测量仪器为GSM06,该仪器为地 1002003004005006007008009001000 选代次数 球物理电磁探测领域事实上的标准，选择该地的主要原因 一：速代值一：实际值。 为该地区属前寒武纪玄武岩地区，电导率低，绝缘性好，对大 地电磁脉冲的彩响较弱，易获得质量较好的信号。特别地， 图2Y的选代收效情况 选择这样的场地，对保留脉冲的细节部分特别有利。从图6 2.0 为可以看到明显的脉冲干扰。图7为实际值及估计的幅度 1.8 概率分布(amplitude probability distribution,APD)图(X轴为 16 -lg(-lgP(x>x),Y轴单位为B),其中，圆周代表实 1.4 际数据，曲线代表通过本文方法估计出的结果。预测的结果 是 为a=1.5,w=(0.36,0.64),y=8.49,d2=30.95。从图7中 可以看出，MCMC方法确实能够有效估计模型参数，具有较 0.8 高精度，同时也再次验证了对称α稳定分布和高斯分布混合 0.6 模型对大气低频噪声的适用性。 01002003004005006007008009001000 600 选代次数 500 一：透代值：：实际值。 400 图3：的运代收效情况 w 100 06 -100 05 -200 -300 .400 03 500 1000150020002500300035004000 采样点 图6大气低频噪声 0.1 0 1002003004005006007008009001000 60 选代次戴 —:透代值；一…：实际值。 50 40 图4助的速代收敛情况 30 50 20 45 3 0 ooogo'o 2 -10 0 0 1 -30 5 01002003004005006007008009001000 0.010.1151020304050607080 90 95 98 选代次数 P中x/% :选代值：：实际值。 。:实际值： :估计曲统。 图5d的送代收效情况 图7实际值及估计的APD图 万方数据
·1244· 系统T程与电子技术 第34卷 -●—_●_——-—_—-—_—●___-●—■———_—__-——●————●—●—●_●—●___—————●—●_————■■——_—____●-__I一一III●————_-__●●●—_—-_______ 迭代次数 一：迭代值；…-：实际值。 图2 y的迭代收敛情况 选代次数 ——：迭代值；～．：实际值。 图3 a的迭代收敛情况 4505 篓 聪 翟10 ： 造代次数 迭代值；…：实际值。 ”l的迭代收敛情况 迭代次数 ——：迭代值；…一：实际值。 图5矿的迭代收敛情况 在图2～图5中，容易发现在迭代200次后，各个参量 都能趋于收敛。总体上收敛速度快，结果精度较高。在实验 中，a的更新过程需要计算a稳定分布的概率密度值，以往 一般采用Nolan提出的数值积分进行计算，但耗时长，计算 过大。在本文中，通过对特征函数做快速傅里叶变换，求特 定序列的概率密度值，再插值求各个点的概率密度值，这样 可以大幅降低时间消耗[”]。 图6为实测大气噪声数据。地点为河南某地，采样频率 为40 960 Hz，带宽80 Hz，测量仪器为GSM-06，该仪器为地 球物理电磁探测领域事实上的标准。选择该地的主要原因 为该地区属前寒武纪玄武岩地区，电导率低，绝缘性好，对大 地电磁脉冲的影响较弱，易获得质景较好的信号。特别地， 选择这样的场地，对保留脉冲的细节部分特别有利。从图6 为可以看到明显的脉冲干扰。图7为实际值及估计的幅度 概率分布(amplitude probability distribution，APD)图(X轴为 一lg(一lgP(}zI>凰))，Y轴单位为dB)，其中，圆罔代表实 际数据。曲线代表通过本文方法估计出的结果。预测的结果 为a=1．5，_Il’=(0．36。0．64)，y一8．49，矿一30．95。从图7中 可以看出，MCMC方法确实能够有效估计模型参数．具有较 高精度，同时也再次验证了对称a稳定分布和高斯分布混合 模型对大气低频噪声的适用性。 l～．1‘…．IJ ．．f一．“LL ’…。r”7”’’ ’I…”7『 0$00 1 000l 500 2 000 2 500 3 000 3 500 4 000 采样点 图6大气低频噪声 奄-一。。 -——‘。 O _。h _●●●_ L-_．n C -_ ●●●● 。刍0．冁 … ．j． 黪玉j ●一● _…1一： jO o O 只u|IⅪ0)／％ 0：实蕲值； ：估计曲线。 图7实际值及估计的APD图 鲫湖枷瑚瑚瑚。啪瑚姗枷 掣馨 一 一 一 一 万方数据
第6期 应文威等：大气低频噪声混合模彩的MCMC参数估计 ·1245· IEEE Trans.on Communication 2011,49(3):506-517. 5结论 [5 Schulte S.A fuzzy impulse noise detection and reduction method[J]. 本文设计了基于MCMC算法的大气低频噪声混合模 IEEE Trans.on Image Processing.2006,15(5):1153-1162. 型的参数估计方法。通过利用。稳定分布的乘积性质，将 [6]Lind L.Mufti N.Efficient method for modeling impulse noise in 混合模型中的a稳定分布等价为高斯分布，便于MCMC算 a communication system [J].Electronies Letters,1996,32(3): 法的更新选代。设计了贝叶斯层次模型，充分利用各个参 1440-1441. 数的先验，以Gibbs和M-H混合抽样更新参数。同时在层 [7]Ekrem E.Koca M.Robust ultra-wideband signal acquisition[J]. 次模型中引人了多个超参数，增强了灵活性。在算法中，利 IEEE Trans.on Wireless Communications.2008,7(11):4656 用FFT算法计算α稳定分布的概率密度，加速了算法运行 -4669. 的速度。仿真实验和实测数据表明，该算法迭代收效快，辖 [8]Nikias C L.Shao M.Signal processing with alpha-stable dis- 度高，有很高的实用价值，对日后设计高性能的非高斯低频 tribution and applications[M].New York:Wiley,1995. [9]Wang X D.Poor H V.Robust multi-user detection in nop- 接收机起了重要的铺垫作用。 Gaussian channels [J].IEEE Trans.on Signal Processing 参考文献： 1999,47(2):289-305. [10]Samorodnitsky G,Taqqu M.Stable non-Gaussian random [1]Abraham D A.Detection-threshold approximation for non- process:stochastic models with in finite wariance [M].New Gaussian backgrounds [J].IEEE Journal of Oceanic Engineer- York:Chapman Hall,1994. ig,2010,35(8):335-341. [2]Goken C.Gezici S,Arikan O.Optimal signaling and detector [11]Godsill S,Kuruoglu EE.Bayesian inference for time series design for power-constrained binary communications systems o- with heavy-tailed symmetric alpha stable noise process[C]/ ver non-Gaussian channels [J].IEEE Trans.on Communica- Proc.of the Applications of Heavy Tailed Distributions in tions Letters,2010,14(2):100-111. Economics.Engineering and Statistics,1999:101-128. [3]Middleton D.Non-Gaussian noise models in signal processing for [12]Salas-Gonzales D.Kuruoglu EE.Modeling with mixture of telecommunications:new methods and results for class A and symmetric stable distributions using Gibbs sampling].Sig- elass B noise models[J.IEEE Trans.on Information Theory, nal Proce3ing,2010,17(2):774-783. 1999,45(4):1129-1149. [13]Menn C,Rachev S T.Calibrated FFT-based density approxi- [4]Stopler D.Zamir R.Capacity and error probability in single-tone mations for stable distribution [Com putational Statistics and multitone multiple access over an impulsive channel [J]. and Dat4 Analysis,2006,50(8):1891-1904. 万方数据
第6期 应文威等：大气低频噪声混合模型的MCMC参数估计 ·1245· 5结论 本文设计了基于MCMC算法的大气低频噪声混合模 型的参数估计方法。通过利用a稳定分布的乘积性质，将 混合模型中的a稳定分布等价为高斯分布，便于MCMC算 法的更新迭代。设计了贝叶斯层次模型，充分利用各个参 数的先验。以Gibbs和M—H混合抽样更新参数。同时在层 次模型中引入了多个超参数，增强了灵活性。在算法中，利 用FFT算法计算。稳定分布的概率密度．加速了算法运行 的速度。仿真实验和实测数据表明，该算法迭代收敛快，精 度高。有很高的实用价值，对日后设计高性能的非高斯低频 接收机起了重要的铺垫作用。 参考文献： [13 Abraham D A．Detection-threshold approximation for 11011- Gaussian backgrounds[J]．IEEEJoumal ofOceanic Engineer— ing，2010，35(8)：335—341． [z]Goken C，Gezici S，Arikan O．Optimal signaling and detector design for power-constrained binary communications systems o— ver non—Gaussian channels[J]．IEEE Tram．on Commanica— tions Letters，2010，14(2)：100—111． [33 Middleton D．Non-Gaussian noise models in signal processing for telecommunications：new methods and results for class A and class B noise models[-J]．IEEE Tram．on In formation Theory， 1999，45(4)11129一1149． [4]Stopler D，Zamir R．Capacity and error probability in single-tone and multitone multiple access over all impulsive channel[J]． IEEETrans．onCommunication，2011t 49(3)：50fi一517． [s]Schulte&A fuzzy impulse noise detection and reduction method[J]． IEEETram．onImage Processing，2006，15(5)：1153—1162． [6]Lind L，Mufti N．Efficient method for modeling impulse noise in a communication system口]．ElectronicsLetters，1996，32(3)： 1440一1441． [7]Ekrern E。Koca M．Robust ultra—wideband signal acquisition[J]． 1EEE Trans．on Wireless Communications。2008，7(11)：4656 —4669． [83 Nikias CL，Shao M．Signal processing with alpha—stable dis— tribution and applications[M]．New York：Wiley．1 995． [9]Wang X Dt Poor H V．Robust multi—user detection in non— Oaussian channels[J]．IEEE Trans．on Signal Processing。 1999，47(2)：289—305． [10]Samorodnitsky G，Taqqu M．Stable non-Gaussian random process：stochastic models with抽finite variance[M]．New Yorki Chapman＆Hall，1994． [11]Godsill S，Kuruoglu E E．Bayesian inference for time series with heavy-tailed symmetric alpha stable noise process[c]∥ Proc．of the Applications of Heavy Tailed Distributions in Economics·Engineering and Statistics，1999 l 101—128． [12]Salas-Oonzales D，Kuruoglu E E，Modeling with mixture of symmetric stable distributions using Gibbs sampling[J]．Sig— hal Processing。2010·17(2)：774—783． [13]Menn C-Rachev S T．Calibrated FFT_ba∞d density approximations for stable distribution[J]．Computational Statistics and Data Analysis，2006，50(8)l 1891—1904． 万方数据