小样本条件下基于矩阵乘法和秩分析的LDPC参数估计方法

刘 倩，张 昊，宋莹炯，王 刚

（1.信息工程大学信息系统工程学院,河南郑州 450001；2.信息工程大学密码工程学院,河南郑州 450001）

1 引言

信道编码类型众多，其中的LDPC（Low-Density Parity-Check）码［1，2］由于其逼近香农限的传输性能和强纠错能力，以及适用于硬件并行运算的置信度传播解码算法被广泛应用于各种通信标准中.LDPC 码不同于其他线性码，（如Hamming 码、循环码、卷积码和Turbo 码等），这类编码由其稀疏校验矩阵来确定，一般具有较长的码长.由于计算复杂度或空间复杂度的原因，原本适用于短码的编码参数识别算法［3］很难再对LDPC码起作用.因此，对LDPC 码的识别分析方法自成一派.

目前对LDPC 码的参数识别算法主要集中于闭集识别［4～8］，即信号接收方手里已有若干种LDPC 编码方式构成的集合，且发送方选取的编码方式为集合中的一种，只需要利用接收码字序列将编码方式挑出即可.开集识别时，由于信号截获方没有任何先验知识，识别难度更大导致研究成果较少［9～16］.已有工作中有基于对偶码的方法，通过寻找具有一定的汉明重量（或小汉明重量）的对偶码获得校验向量，与此同时恢复码长、码率和起点［9］.也有假设编码长度、码率及起点均已知，在此基础上恢复稀疏校验矩阵［10～13］.或者假设截获信号序列长度足够长，在此基础上利用秩分析法（分析构造码字矩阵的秩率［14］或近似秩率［15，16］）估计码长和起点，这两种方法均是基于高斯列消元方法.近两年还出现了利用秩的概率分布获取正确的编码长度和起点的算法［17，18］.算法中需要多次从编码矩阵中随机获取方阵，计算它们的秩并得到其经验分布规律.如果其与已有的同阶随机方阵的秩的分布规律相同，则被认为是随机方阵，此时码长和起点不正确.否则便被认为具有某种代数结构，此时码长和起点正确.关于已有工作［14～18］的详细介绍和对其原理及缺点的具体分析将在节2.4节中给出.

以上方法均要求截获比特流的长度是足够供信息截获者进行分析的，但是，在信息对抗过程中，信息对抗方是被动接收信息，并不能保证截获所得数据量充足，此时已有的方法全部失效.本文在截获比特数量远少于现有方法所需的比特数量的情况下（不含或含少量错误比特），基于码字空间与其对偶空间的正交性、完整码字比特间的线性相关性、矩阵乘积的秩不大于每个因子矩阵的秩等性质，提出了一种估计LDPC 码的码长、起点等参数的矩乘秩减算法，并对算法的计算复杂度和正确识别概率进行了分析.

2 数学模型和理论基础

2.1 符号说明

接下来在二元域F2上展开问题的研究.用C表示所采用的纠错码，n，k分别表示纠错码的码长和维数.信息向量和码字向量等矢量分别用粗体u和c表示，变量用小写斜体字母表示.大写斜粗体字母A和AT分别表示码字矩阵及其转置，矩阵A的零空间记为Kel(A)，A(i，：)和A(：，i)分别表示矩阵A的第i行和第i列.W和W⊥分别代表纠错码码字空间和其对偶空间.d表示正确的起点位置，l0和t0分别表示利用我们提出的算法得出的码字长度和起点位置，l和t是在一定区间内变化的变量，将截获比特流在不同的起点t下按列数l排列成的码字矩阵记为Dl，t.

2.2 码字空间的封闭性及与其对偶空间的正交性

一般的线性分组码的编码方式是由生成矩阵定义：给定生成矩阵Gk×n，输入由k个比特组成的信息块u=(u1，u2，…，uk)，得到一个完整的码字c=(c1，c2，…，cn)=u·Gk×n.已知Gk×n的行向量线性无关，分别记之为g1，g2，…，gk，则

可知所有码字c均表示为生成矩阵行向量的线性组合，因此生成矩阵行向量张成一个k维线性空间，称为码字空间，记为W.当Gk×n为系统形式(Ik×k，Pk×(n-k))时，则生成系统码.

若有两个码字ci=(ci1，ci2，…，cin)=ui·Gk×n和cj=(cj1，cj2，…，cjn)=uj·Gk×n，则

ci⊕cj=(ci1⊕cj1，ci2⊕cj2，…，cin⊕cjn)=(ui⊕uj) ·Gk×n仍然是一个合法码字，因此码字空间关于F2上的加法运算⊕封闭.

LDPC 码是一种特殊的线性分组码，仅由一个包含很少个1 的校验矩阵H(n-k)×n的零空间所定义.因此当给定k个信息比特(u1，u2，…，uk)时，需要利用校验矩阵求解n-k个校验比特(uk+1，uk+2，…，un)，即系统码字c=(c1，c2，…，cn)=(u1，u2，…，uk，uk+1，uk+2，…，un) 满足c·HT=0，因此LDPC码的码字空间不仅关于⊕封闭，并且与其对偶空间具有正交性.而LDPC 码的识别过程也正是基于这种关系.

2.3 高斯若当列消元法

高斯若当列消元法GJETP（Gauss-Jordan Elimination Through Pivoting），通过行置换和列变换将一个矩阵变换为下三角矩阵，文献［14～16］均是在此基础上给出各自的算法.先给出F2上GJETP 法实施的具体过程.

对于一个s×l阶的矩阵A，i的变化范围是从1 到min｛s，l｝.

（1）如果A的第i列的第i个元素为0，并且存在最小的j(j＞i)使得第j列的第i个元素为1，则交换A的第i列和第j列；

（2）如果A的第i列的第i个元素为0，并且不存在第j(j＞i)列其第i个元素为1，但有最小的j(j＞i)使得A的第j行中第i个元素为1，则交换A的第i行和第j行；

（3）如果A的第i列的第i个元素为1，则将其加到第i行中所有列数大于i的位置上为1 的列上.

2.4 已有工作和数学模型

在一个完整码字中，校验比特是某些信息比特的线性组合.如果接收比特足够多，且由无误码比特序列构建的矩阵A的每一行恰好为一个完整的码字，则其列与列之间的比特服从相同的线性约束关系，因此至多只有k列线性无关，A的秩至多为k.在通过GJETP 将A变换为有n-k列全零列的矩阵的过程中，设列变换矩阵为Q.若编码是系统码，则存在Q=(Q1Q2)使得AQ1=VM×k，AQ2=0M×(n-k)，记作

图1 当l ≠βn和l=βn时，矩阵A在形式上的差别

为了更好的做出对比，将矩阵的秩除以l得到归一化秩率［14，15］，便可利用这两种情况下归一化秩率的不同表现来区分预设的编码长度和起点是否正确.为了进一步提高算法的容错能力，降低错误比特对码字矩阵秩的影响，文献［16］改为寻找矩阵的“近似秩”.其做法是利用GJETP 将码字矩阵化成下三角阶梯型矩阵后，寻找矩阵中列重大于某一阈值的列，其数目便为矩阵的近似秩.同样地，近似秩归一化后得到近似秩的归一化秩率，将使得其取得最小值的矩阵的列数和起点视为正确的参数.这里把文献［14，15］中的方法叫做秩准则，把文献［16］中的方法叫做近似秩准则，这两种方法都为基于GJETP的秩分析法.

但在利用上述秩分析法获取码字长度和起点等参数的过程中，一般要求矩阵A的行数s大于其列数l.如果截获码字个数较少使得构建的码字矩阵的行数远远小于列数，则秩分析法失效.以秩准则为例，设码字为(n，k)线性分组码，截获的码字序列为Z，其中含有M个比特.假定码字长度和起点分别为l和t，构造码字矩阵Dl，t，其中Dl，t的第i个行向量为在≤k的情况下分析问题.此时，如果l≠βn，相应的矩阵Dl，t为随机矩阵，由于其行数远小于列数，因此有rank(Dl，t)≤k.而当l=βn时，矩阵Dl，t是结构矩阵，其行向量是至多k个线性无关的行向量（生成矩阵G的行向量）的线性组合，因此也有rank(Dl，t)≤k.此时无法由秩准则区分假定的码字长度和起点是否正确.

近似秩准则［16］在估计编码参数问题上表现突出.但是近似秩准则有一个要求，那就是截获比特数量必须非常大，使得构造的码字矩阵为“高瘦型”.矩阵越高，近似秩才越准确，从而得出正确参数的概率也越大，并可以利用近似秩进一步确定校验关系.而当码字矩阵为方阵或者“矮胖型”矩阵时，则无法计算近似秩.因此近似秩准则在截获比特数量较少时完全失效（高瘦型矩阵指的是其行数远远大于列数，矮胖型矩阵指的是其行数远远小于其列数）.

在利用随机方阵秩的分布获得编码参数的做法［17，18］中，需要在构造的码字矩阵Dl，t中随机选取行向量构造出许多个方阵.计算方阵的秩得出其经验分布函数，然后将其与同阶随机方阵的秩的分布规律进行比较.如果分布规律相同，则判定预设码字长度和起点不正确，如果不同则认为码字长度和起点正确.但在构造码字矩阵为“矮胖型”的情况下，其行向量不可能构造出方阵，因此算法［17，18］也失效.另一方面，算法需要做大量实验统计不同设定码长和起点下方阵的秩的分布规律，因此其计算量也是一个天文数字，在码长较大时并不实用.

3 本文算法、原理及复杂度分析

在给出本文具体算法之前，先以定理的形式给出一个乘积矩阵的秩的性质［19］.

定理1设矩阵A和B都是l阶方阵，则有rank(AB)≤min{rank(A)，rank(B)}成立.容易将此结论推广到n个方阵相乘的情形.即对n个方阵A1，A2，…，An，有

3.1 本文算法

为了解决在截获比特数量较少的情况下编码参数的识别问题，本文改变研究矩阵的秩的方式.注意到不管码字个数有多么少，一个完整码字中的比特间的线性约束关系总是存在的.因此，我们不去关心全部的校验关系，而将注意力放在寻找满足同一个校验关系的比特上.在l=βn的情形下，设h=(h1，h2，…，hn) ∈W⊥，令集合

从矩阵Dl，t中随机选取列构建一个阶的方阵Sl，设选取列的列标集合为θ.若有τ(h) ⊂θ，则Sl不满秩.由于LDPC 码的校验向量的稀疏性，可知满足同一个校验的比特数目较少，只要这几个比特所在列被选入Sl，Sl就不满秩.若是几个校验关系中的比特所在列均被选入Sl，则Sl秩亏的更多，因此使得Sl秩亏的条件较容易达到.当从Dl中随机选取p个方阵，分别记为Sl1，Sl2，…，Slp，作乘积得Sl1Sl2…Slp，由式（3）可得

这使得乘积矩阵Sl1Sl2…Slp的秩较小的概率大大提升.若l≠βn，则Dl，t为随机矩阵，其列之间并不存在线性结构.随机从中选取列组成的方阵Sl仍然是随机矩阵，因此乘积矩阵Sl1Sl2…Slp也是随机矩阵，其秩在一般情况下不会较小.在此分析基础上，为了更好地区分两种情形，定义归一化秩率函数

在具体实施过程中，先估计l0，再将t0确定，提出了下面的矩乘秩减算法，如算法1所示.

注1关于随机选取方阵的个数p的选择问题.由式（4）可知，随着p的增大，乘积矩阵的秩是单调不增的.只要出现一个乘因子矩阵的秩比较小，那么乘积矩阵的秩就会很小.但若Sl1，Sl2，…，Slp都是随机矩阵，它们的乘积也是随机矩阵，出现秩较小的概率较低.因此，随着p的增大，在预设码长正确和不正确两种情况下，乘积矩阵的秩的差距进一步拉大，码长和起点的正确识别概率也会提升（在理论分析部分详细介绍）.但如果p的取值过大，则会使得计算矩阵乘法的次数较多，计算量显著增加.所以，需要在识别概率和计算复杂度之间做一个权衡，选择合适的p值.

注2关于起点所在位置对算法的识别概率的影响问题.在估计码字长度的过程中，若预先设定的码字长度l0正确，t0未知，则码字矩阵Dl0的每一行由上一个码字的l0-t0个比特和下一个码字的t0个比特构成.当t0的值靠近0 或者l0时，由于大部分比特在同一个码字中，随机选取的方阵中出现具有线性相关关系的列的可能性较大，此时更容易得到秩较小的矩阵，识别概率较高.当t0的值靠近时，每一个行向量中约一半的比特分别位于在上一个码字和下一个码字中，而两个码字中的比特之间一般不具有线性相关关系，因此随机选取的方阵秩较大的概率较高.此时，与非正确起点下乘积矩阵的秩区分度不大，识别概率较低.

3.2 算法理论分析

设LDPC 码稀疏校验矩阵为H，其行向量记为h1，h2，…，hr，记w(h)=，w(·)代 表向量的Hamming 重量.我们想要在下面两个假设检验中做判决

下面来考虑正确识别概率，这需要二元域上随机方阵的秩的分布规律.由于二元域上随机方阵秩的具体分布规律目前还未知，已有的仅仅是当方阵阶数趋于无穷大时随机方阵秩的分布规律［20］.因此，利用多次蒙特-卡罗实验将有限阶随机方阵的秩的经验分布函数得出.在10000 次实验中，为了映照后面实验部分选择的编码类型，图2（a）中给出了280阶随机方阵的秩的经验分布函数.图2（b）为阶数从200变化到280的随机方阵的秩取不同值时的频率.从图中可知不管方阵阶数如何变化，l阶随机方阵的秩极大概率都落在区间[l-2，l]上，其中以l-1比例最大.

图2 随机方阵的秩的分布规律

考察正确检测概率

假定随机选取的方阵之间的相关性不大，则可认为其是相互独立的，那么正确检测概率可表示为

由于只要有某两个τ(hi) ⊂θ(i=1，2，…，r)，便有

式（10）中w(h)为稀疏校验向量中重量的最大值.可知随着w(h)变小，Pcd变大，随着M变小，Pcd变小.LDPC 码的校验矩阵为稀疏的且没有4 环，校验矩阵的行重一般来说远远小于码长，因此提出的算法特别适用于LDPC码的参数识别.

3.3 算法计算复杂度分析

当l在[lmin，lmax]中变化时，对于给定的l，计算p个阶方阵的乘积需要的乘法数量为

加法数量为

因此，要估计出编码长度l0，一共需要的乘法数量为

加法数量为

比较运算的数量为(lmax-lmin)次.固定编码长度，让起点t在[1，l0]中变化.要得出起点t0，类似的操作下至多需要的乘法数量为

加法数量为

比较运算的数量为l0-1次.

4 模拟仿真实验

4.1 实验设置

在模拟仿真实验中，均以IEEE802.16e 标准中的（576，288）LDPC 码为例.在截获比特流长度M固定的情况下，在截获数据含少量误码和无误码两种情形下，我们分别考察了选取乘积矩阵的个数p对识别成功概率的影响.令M=132480，分别在误比特率Pe为0 和0.0001时令乘积矩阵的个数从1连续变化到10.

由于有这样的结论：若长向量线性相关，则从中截取的短向量一定线性相关.而短向量线性相关时，其延长得到的长向量未必线性相关，因此，M对识别概率也是有影响的.一般情况下，由上面结论可知当比特流长度越长时，码字长度和起点的正确识别概率也会越高.因此在p固定时，我们考察在有误码和无误码情况下M的变化对识别成功率的影响.在比特流长度M=q·n时，令q在180和320之间变化，将本文方法与文献［14～17］在无误码和Pe=0.0001两种情况下进行对比.

为了考察本文算法在不同的误比特率下的表现，我们在截获比特数量M=132480，p=10的情况下，令误比特率Pe从0.0001按间隔为10-4变化到0.0008，考察Pe的变化对识别成功率的影响.还研究了在码字长度确定的情况下，起点t0的位置对归一化秩率的影响，印证了我们在第3节中的分析.

4.2 实验结果分析

从图3中可以看出，码字长度和起点的正确识别概率确实随着p的增大而增大.原因是列向量的选取具有随机性，如果满足几个校验关系的列都被选中，则方阵的秩较小.但若都未被选中，则方阵的秩可能比随机矩阵的秩大，也可能比随机矩阵的秩小.因此，当p较小时不确定性较大，随着p的增大这种不确定性的影响渐渐被减弱.在图3中还可以看出当p=1时，有误码情况下的识别概率稍大于无误码情况下的识别概率，正是这种原因的直观体现.

图3 选择作乘积的矩阵个数p对识别概率的影响

由于识别概率直接受到乘积矩阵的秩的影响，为了更加深入考察p的取值对乘积矩阵的秩的影响，在l=n和l≠n两种情况下按照p取值从1到10的顺序分别做了十次实验.将实验所得的秩求平均值，称其为平均秩.表1给出了两种情况下平均秩在不同的p值下的取值.可以看出，不管预设码长正确不正确，平均秩都随着p的增大而减小，并且码长正确时的平均秩整体上要比码长不正确时的平均秩小，这与我们提出算法的思想相符.

表1 十次实验中，p对乘积矩阵平均秩的影响

从图4（a）中可以看出，在无误码且q≤300时，秩准则方法［14，15］失效，而本文算法在p=10且q=250时正确识别率达到100%.在误比特率Pe=0.0001且q≤309时，秩准则失效，而本文算法在p=10且q=260时正确识别率达到100%.4（b）中展示了利用近似秩准则［16］以及秩分布方法［18］在Pe=0.0008 时的结果，4（b）可作为4（a）的延展.可以看出虽然它们的容错能力强于本文方法，但它们需要码字最少个数时的q也要远大于本文方法.图5 中显示，在这个变化过程中码长和起点的正确识别概率从89%下降到20%.虽然本文算法容错能力仅在10-4量级，但是在截获比特数量较少而其他现有算法均失效的情况下，这样的性能也是难能可贵的.图6 给出了在Pe=0.0001 且预先设定码长l变化时，从相应的分析矩阵中随机选取p=10个方阵，其乘积矩阵的归一化秩率的值.为了更清楚的显示归一化秩率随l变化的取值情况，我们截取了正确码长l0=576 所在的一段区间［558，594］.可以看出当l=576 时归一化秩率的值最小.

图4 截获比特流长度M=q·n的变化对识别概率的影响

图5 误比特率较低时，误比特率Pe对识别概率的影响

图6 乘积矩阵的归一化秩率随着l的变化的取值情况

在确定码字长度之后，分析矩阵中的列都已经固定，差别仅仅是列的位置的不同.此时，正确起点与不正确起点所对应的平均归一化秩率之间的差距进一步缩小，正确起点t0的寻找也不是一件容易的事.为了进一步增加确定性，确保正确起点被选出，可以让起点在[t+1，t+l0+1]中变化.如果某两个位置相差l0，且对应的平均归一化秩率相较于其他值均较小，则可认定此位置为正确起点.表2 给出了（576，288）LDPC 码在第1位即为正确起始位，无误码且M=132480，p=10的情况下，t变化时的平均归一化秩率的部分取值.可以看到，标黑的第1 位和第577 位的值相对其他值较小，同时也可以看出当假定起始位置接近时的平均归一化秩率明显要大.这也印证了在第3节中关于t0位置对平均归一化秩率的影响的分析.

表2 l0已确定，t变化时，归一化秩率的取值情况

5 总结

本文在截获少量比特的背景下，此时其他算法由于所需比特数量众多而全部失效的情况下，对LDPC 码的起点和码长进行了识别.利用完整码字中比特间具有线性相关性和矩阵乘积的秩不大于因子矩阵的秩等性质，从构造的码字矩阵中随机选取列，得到一系列方阵并作乘积，再选出使得乘积矩阵的归一化秩率最小时的分析矩阵列数作为码长.而后，利用类似的方法对列数固定、起点发生变化时构造的分析矩阵进行处理，并且改变分析矩阵的行数多次实验.将使得平均归一化秩率最小时的位置确定为码字起点.我们称之为矩乘秩减算法.与传统算法相比，达到相同的正确识别率时本文算法所需数据量至少减少25%，但计算量没有明显增加.