阅读视图

发现新文章,点击刷新页面。
🔲 ☆

摩托驾驶证科目一易错点

限速相关

  • 无中心线,城市道路30公路40
  • 有中心线,城市道路50公路70
  • 划有两条以上车道的,城市道路70
  • 封闭的机动车专用道,限速80
  • 运载易燃易爆、剧毒、放射性等,在封闭的机动车专用道和公路上,限速60;在城市道路上限速50
  • 限30
    • 遇雾、雨、雪等能见度50m以内,限速30km/h
    • 机动车进出非机动车道时,限速30
    • 通过急弯路时,限速30
    • 三轮、拖拉机、电瓶车等,限速30

距离相关

  • 夜间会车,距离150m意外,将远光灯改为近光灯
  • 车辆故障,车后50-100m设置警告标志
  • 口5站3:交叉路口、隧道口50m内不能停车;公交站、消防栓、消防队门前30m内不能停车
  • 距离掉头点50-100m驶入最左侧车道

扣分相关

骑摩托车不戴头盔,扣1分

3分6分9分12分
一般城市道路上逆行高速上违法占用应急车道高速上停车高速上逆行
不规避校车不按交通信号灯指示通行
驾驶证被暂扣期间驾车准驾不符
未取得校车驾驶资格驾驶校车
伪造、变造机动车驾驶证
不悬挂车牌、遮挡号牌使用伪造、变造机动车号牌
饮酒后驾车
超速50%以上

驾驶证相关

  • 多久不得申请驾驶证:假一吊二撤三醉五逃终身,撤、毒、代三年

    • 虚假材料,一年内不得重新申请驾驶证

    • 以欺骗、贿赂等手段取得驾驶证,被依法撤销的,三年内不得重新申请驾驶证

    • 考试过程中有贿赂、舞弊行为的,一年内不得再次申领驾驶证

    • 注射毒品的,三年内不得申请驾驶证

  • 饮酒后驾车,暂扣6个月驾驶证,罚款1000-2000元

  • 超速50%以上,吊销驾驶证,罚款200-2000元

罚款相关

  • 转让不登记,罚200
  • 隐瞒、欺骗手段补领机动车驾驶证,罚款200-500
  • 提供虚假材料申请驾驶证,处500元罚款
  • 代罚三倍不超五万,代审三倍不超两万
  • 饮酒后驾车,暂扣6个月驾驶证,罚款1000-2000元
  • 超速50%以上,吊销驾驶证,罚款200-2000元
  • 交通事故,机动车已采取必要措施,赔偿金额,高架桥上5%不超过一万,普通道路10%不超过五万
  • 违反倒车规定,未按交警指挥行驶,抛洒,罚100
  • 改装助力,罚款50

拘役相关

  • 交通事故致人重伤、死亡,三年以下
  • 致人重伤、死亡并逃逸,3-7年

时限相关

  • 驾驶证有效期超过一年以上未换证的,驾驶证会被注销
  • 驾驶证有效期满之前,90日内申请换证
  • 申请人信息变化的,要在30日内申请换证
  • 摩托车注册登记之日起,10年以内,每2年检验1次;超过10年的,每1年检验1次

灯光相关

  • 雾天:雾灯、危险报警闪光灯;不能开启远光灯,浓雾会反射灯光

路标相关

  • 弯道:一急二反三连续
  • 路面不平两个包
  • 栅栏有人,火车头无人
🔲 ☆

三维旋转、欧拉角、四元数

关于三维旋转、欧拉角和四元数,以前只是简单的调用库函数做转换,并没有深入理解,最近在网上找资料发现还是良莠不齐,很多文章写得比较杂乱,这里根据自己学习后的理解整理出一篇文章以供参考。

1. 欧拉角

欧拉角描述了三维空间中的旋转,想象有一架飞机悬在空中朝着正东方向水平静止,现在要将其姿态原地变换到任意一个其他方向,从最接近人的直观来看,在机头原始位置\(A\)和旋转后位置\(A'\)连一条线即为最短路径,沿着该最短路径旋转飞机即可。如果用数学语言描述,那么就是假设飞机中心位置为\(O\),那么\(AOA'\)构成一个平面,以该平面法向量为轴旋转飞机,即可将飞机从原始姿态旋转到新的姿态。尽管这种方法很直观,但是一开始欧拉并没有采用这种表示方法,或许是觉得不够规整,旋转轴的确定和表示稍显麻烦。那么欧拉是如何做的呢?如下图所示[1],我们知道笛卡尔坐标系有\(xyz\)三个轴,那么绕着这三个轴分别做旋转不就可以了吗?这时只需要给出每次旋转的角度即可。

更具体清晰的描述如下:首先假设世界坐标系\(xyz\),飞机自身的局部坐标系为\(XYZ\)(刚体坐标系),世界坐标系始终保持固定不变,而飞机局部坐标系在飞机旋转过程中会发生变化。假设飞机初始坐标系为\(X_1Y_1Z_1(=xyz)\),想要旋转到一个新的姿态\(X_2Y_2Z_2\),那么欧拉给出的一个结论是:任意两个姿态之间,只需要绕\(x,y,z(X,Y,Z)\)轴做3次旋转即可完成。不过,具体到旋转过程中,绕哪几个轴旋转的先后顺序和角度都会对最终结果有影响<spanclass="hint--top hint--rounded" aria-label="欧拉角细节/旋转顺序/内旋外旋- 知乎(zhihu.com) ">[2]

  1. 内旋还是外旋:可以注意到上面动图中首先绕z轴旋转,第一次旋转之后局部坐标系的\(XY\)轴与世界坐标系的\(xy\)轴已经不对齐了,因此第二步是以局部坐标系的轴进行旋转还是以世界坐标系为轴进行旋转?前者被称为内旋,后者为外旋。显然其他参数相同的情况下,两种方式得到的结果不同。
  2. 旋转轴的顺序:欧拉角实际上有两大类,一类是Proper Eulerangles,其旋转顺序包括 (z-x-z, x-y-x, y-z-y, z-y-z, x-z-x,y-x-y),第二类是Tait–Bryan angles,包括 (x-y-z, y-z-x, z-x-y, x-z-y,z-y-x, y-x-z)。可以看到Proper Euler angles只涉及两个转轴.而Tait–Bryanangles涉及三个转轴。目前机器人以及SLAM领域中常用的是Tait–Bryanangles。后面的论述也都采用Tait–Bryan angles。
  3. 旋转的角度:指定绕 \(xyz(XYZ)\)轴三个的旋转角度分别为 \((\alpha,\beta,\gamma)\),需要注意的是如果旋转轴的顺序不同,即使每个轴对应的旋转角度不变,最终的姿态角也可能不同。

最后给一个小的总结,也就是说要想用欧拉角的方式描述一个三维旋转,实际上需要5个参数,即:1)内旋or外旋;2)旋转顺序;3)每次旋转的角度。

另外纠正一个重要的误解:欧拉角并不是描述一个静态的姿态角,而是描述一种过程!(相信很多人和我一样一直没意识到这点)姿态本身是唯一确定的,但是变换到这个姿态的过程是不唯一的,因此欧拉角与内旋/外旋、旋转顺序都密切相关!

2. 旋转矩阵

我们知道旋转可以用矩阵表示,如果考虑外旋的话,那么假设某个点在世界坐标系中的初始位置为\(\boldsymbol{p}_0\),如果跟随世界坐标系按照某个Tait–Bryanangles做外旋,不妨设为 \(x(\alpha)\toy(\beta)\toz(\gamma)\),其旋转后的坐标(在世界坐标系中)是可以表示为 \(\boldsymbol{p}'= \boldsymbol{R} *\boldsymbol{p}_0\),其中 \(\boldsymbol{R}\)为旋转矩阵,那么这个旋转矩阵可以写为 \(\boldsymbol{R}_{\rmextrinsic}=\boldsymbol{R}_z(\gamma)\boldsymbol{R}_y(\beta)\boldsymbol{R}_x(\alpha)\),其中\[\begin{align}\boldsymbol{R}_x(\alpha) &= \left[\begin{array}{ccc}1 & 0 & 0 \\0 & \cos (\alpha) & -\sin (\alpha) \\0 & \sin (\alpha) & \cos (\alpha)\end{array}\right] \\\boldsymbol{R}_y(\beta) &= \left[\begin{array}{ccc}\cos (\beta) & 0 & \sin (\beta) \\0 & 1 & 0 \\-\sin (\beta) & 0 & \cos (\beta)\end{array}\right] \\\boldsymbol{R}_z(\gamma) &= \left[\begin{array}{ccc}\cos (\gamma) & -\sin (\gamma) & 0 \\\sin (\gamma) & \cos (\gamma) & 0 \\0 & 0 & 1\end{array}\right]\end{align}\]外旋比较简单,但内旋的旋转矩阵就比较复杂,首先给出一个结论:如果内旋操作为\(Z(\gamma) \to Y'(\beta)\toX''(\alpha)\),那么其对应的旋转矩阵为 \[\boldsymbol{R}_{\rmintrinsic}=\boldsymbol{R}_{X''}(\alpha)\boldsymbol{R}_{Y'}(\beta)\boldsymbol{R}_{Z}(\gamma)=\boldsymbol{R}_z(\gamma)\boldsymbol{R}_y(\beta)\boldsymbol{R}_x(\alpha)\] 也就是和前面外旋 \(x(\alpha)\toy(\beta)\to z(\gamma)\) 的结果是完全一样的。也就是说:

每种特定顺序的外旋等价于其相反顺序的内旋,例如 \(x(\alpha)\to y(\beta)\to z(\gamma)\) 与\(Z(\gamma) \to Y'(\beta)\toX''(\alpha)\) 的旋转是等价的。

上面 \(\boldsymbol{R}_{\rmintrinsic}\) 的表达式中,第1次旋转绕 \(Z(=z)\) 轴旋转,因此 \(\boldsymbol{R}_{Z}(\gamma)=\boldsymbol{R}_z(\gamma)\),然而在做第2次旋转的时候,此时初始的\(Y\) 轴已经被旋转到了 \(Y'\) 处,因此绕 \(Y'\) 轴的旋转需要根据 \(Y'\) 的方向来确定,\(\boldsymbol{R}_{Y'}(\beta)\)表达式推导较为繁琐,而第3次旋转类似,其旋转轴变为了 \(Z''\),表达式则更为复杂。但是经过推导,会发现他们神奇地和\(\boldsymbol{R}_z(\gamma)\boldsymbol{R}_y(\beta)\boldsymbol{R}_x(\alpha)\)完全等价!实际上这并不是巧合,如何理解呢?下面我将参考博客<spanclass="hint--top hint--rounded" aria-label="Extrinsic& intrinsic rotation: Do I multiply from right or left? | by DominicPlein | Medium">[7]给出解释,稍微有点绕,不想深究细节的可以只记住上面的结论,跳过这段。

首先定义一些符号,有两个坐标系,即世界坐标系 \(W\) 和旋转后刚体坐标系 \(B\),我们一开始只考虑单步旋转,也即 \(W\)绕某个轴(外旋,同时也是内旋)旋转1次即为 \(B\),刚体上有一个点 \(A\) 跟着刚体坐标系做旋转,旋转后变为点\(A'\)。记旋转前后点 A在世界坐标系\(W\)中的坐标为 \({}^{\rm W}\boldsymbol{p}_{A}, {}^{\rmW}\boldsymbol{p}_{A'}\),在 \(B\) 中的坐标为 \({}^{\rmB}\boldsymbol{p}_{A'}\)。那么可以有下面两个等式:

  1. 记从 \(W\)\(B\) 的旋转矩阵为 \({}^{\rm W}\boldsymbol{R}_{\rmB}\)(外旋矩阵的形式,例如如果\(W\)到\(B\)是绕着\(x\)轴旋转\(\alpha\),那么就有\({}^{\rm W}\boldsymbol{R}_{\rmB}=\boldsymbol{R}_{x}(\alpha)\)),那么我们有 \({}^{\rm W}\boldsymbol{p}_{A'} = {}^{\rmW}\boldsymbol{R}_{\rm B} {}^{\rm W}\boldsymbol{p}_{A}\);
  2. 另一方面,由于点 \(A\)跟着刚体坐标系做旋转,因此其在刚体坐标系中的坐标始终保持不变,也就是说我们有\({}^{\rm B}\boldsymbol{p}_{A'} = {}^{\rmW}\boldsymbol{p}_{A}\)

基于以上两个等式,可以得到 \({}^{\rmW}\boldsymbol{p}_{A'} = {}^{\rm W}\boldsymbol{R}_{\rm B} {}^{\rmB}\boldsymbol{p}_{A'}\),注意到此时发生了质的变化,我们把在同一个点两个不同坐标系中的坐标,通过旋转矩阵\({}^{\rm W}\boldsymbol{R}_{\rm B}\)联系在了一起!实际上,矩阵\({}^{\rmW}\boldsymbol{R}_{\rm B}\)有两种含义:1)从 \(W\) 到 \(B\)的旋转矩阵;2)空间中的同一个点,将其在\(B\)坐标系下的坐标映射到\(W\)坐标系下坐标的过渡矩阵

那么再思考一下,所谓的外旋实际上同时也是关于坐标系 \(W\) 的内旋!如果把世界坐标系 \(W\) 直接替换成任意一个刚体坐标系 \(C\),那么我们就有 \({}^{\rm C}\boldsymbol{p}_{A'} = {}^{\rmC}\boldsymbol{R}_{\rm B} {}^{\rmB}\boldsymbol{p}_{A'}\)。还剩下一个问题就是,\({}^{\rm C}\boldsymbol{R}_{\rm B}\)的表达式是什么?实际上,假如 \(C\toB\) 是绕 \(x\) 轴内旋 \(\alpha\),那么就有 \({}^{\rm C}\boldsymbol{R}_{\rm B} =\boldsymbol{R}_{x}(\alpha)\)!(思考一下为什么~)

现在考虑多步旋转,内旋操作为 \(Z(\gamma)\to Y'(\beta)\toX''(\alpha)\),旋转过程中的坐标系分别为 \(W\to B\to C\to D\),某一个点 A经过多此旋转变为 \(A\to A' \to A''\to A'''\),旋转前 \(A\) 的坐标为 \({}^{\rm W}\boldsymbol{p}_{A}\),由于 \(A\)始终跟着刚体坐标系旋转,因此其在局部坐标系中的坐标始终不变,也就是 \({}^{\rmD}\boldsymbol{p}_{A'''}={}^{\rmC}\boldsymbol{p}_{A''}={}^{\rm B}\boldsymbol{p}_{A'}={}^{\rmW}\boldsymbol{p}_{A}\),那么可以有如下关系:

  1. 先在坐标系 \(C\)的视角下,考虑最后一次旋转操作,有 \({}^{\rmC}\boldsymbol{p}_{A'''} = {}^{\rm C}\boldsymbol{R}_{\rm D}{}^{\rm D}\boldsymbol{p}_{A'''} = \boldsymbol{R}_{x}(\alpha){}^{\rm D}\boldsymbol{p}_{A'''} = \boldsymbol{R}_{x}(\alpha){}^{\rm W}\boldsymbol{p}_{A}\)
  2. 然后在坐标系 \(B\)的视角下,考虑倒数第二次旋转操作,有 \({}^{\rmB}\boldsymbol{p}_{A'''} = {}^{\rm B}\boldsymbol{R}_{\rm C}{}^{\rm C}\boldsymbol{p}_{A'''} = \boldsymbol{R}_{y}(\beta){}^{\rm C}\boldsymbol{p}_{A'''} =\boldsymbol{R}_{y}(\beta)\boldsymbol{R}_{x}(\alpha) {}^{\rmW}\boldsymbol{p}_{A}\)
  3. 接着在坐标系 \(W\)的视角下,考虑倒数第三次旋转操作,类似的有 \({}^{\rm W}\boldsymbol{p}_{A'''} =\boldsymbol{R}_{z}(\gamma)\boldsymbol{R}_{y}(\beta)\boldsymbol{R}_{x}(\alpha){}^{\rm W}\boldsymbol{p}_{A}\)。

因此我们就得到了内旋操作 \(Z(\gamma) \toY'(\beta)\to X''(\alpha)\) 对应的旋转矩阵恰好就是\(\boldsymbol{R}_z(\gamma)\boldsymbol{R}_y(\beta)\boldsymbol{R}_x(\alpha)\),和外旋\(x(\alpha)\to y(\beta)\to z(\gamma)\)是完全一样的。

还有另一种解释方法,旋转矩阵为何左乘是相对固定坐标系,右乘是相对当前坐标系?- 知乎 (zhihu.com)

3. 万向锁

提到欧拉角旋转,很多文章会谈到万向锁问题,需要注意的是:欧拉角用来表示物体的旋转状态是没有问题的,万向锁是用欧拉角表示旋转过程中的一种状态,因此其影响到的是连续运动的问题<spanclass="hint--top hint--rounded" aria-label="欧拉角表示旋转会出现的问题——万向锁(GimbalLock)-CSDN博客">[3]。另外,万向锁发生在内旋过程中,外旋过程是不会发生的。万向锁现象:假设内旋\(Z\to Y\to X\),初始刚体坐标系为 \(X_0Y_0Z_0\),其与世界坐标系重合,第1次旋转的时候绕着\(Z_0/z\)轴旋转(任意角度)后得到新的刚体坐标系 \(X_1Y_1Z_1\),然后第2次旋转过程中绕 \(Y_1\)轴旋转了90°/-90°,此时第2次旋转过后的刚体坐标系 \(X_2Y_2Z_2\) 中的 \(X_2\) 轴与世界坐标系中的 \(z\)轴是同向/反向的,那么在做第3次旋转操作的时候,理应绕着此时的 \(X_2\)轴做旋转,但从世界坐标系来看,也就等价于绕着 \(z\)轴旋转,也就是说第1次和第3次旋转是完全可互相替代的,换言之失去了一个旋转的自由度。其他的内旋顺序也是类似的。简而言之,内旋过程中,当第二个旋转轴旋转90°/-90°时,第三个轴的转动效果可以被第一个旋转轴替代。

通常网上文章的另一种解释思路是用万向节,也就是陀螺仪的例子,个人认为这个例子不合适。首先陀螺仪中三轴万向锁中的三个轴并不总是保持互相垂直的,外圈的轴旋转并不会带动内圈轴旋转,因此它和我们直接理解刚体坐标系的旋转根本不是一回事。那么网上广为流传的那个粒子和万向锁有什么联系呢?还是结合刚才内旋\(Z\to Y\to X\)的例子,以及网上流传的众多图片,下面将给出一个解释。

首先对于陀螺仪而言,可以看到其结构如下,中间是一根竖轴穿过一个盘子,而盘子处于高速旋转状态,是陀螺的转子,根据陀螺的定轴性,竖轴也就是陀螺的自转轴在惯性空间内的方向保持不变。[4]

初始时刻陀螺仪的三轴对应于世界坐标系的三轴,绿色可以绕x轴旋转,红色可以绕y轴旋转,蓝色可以绕z轴旋转。

绿色绕x轴旋转红色绕y轴旋转蓝色绕z轴旋转

第一次旋转的时候,紫色绕 \(z/Z_0\)轴旋转,此时内圈带动外圈一起旋转,还没出什么问题;第2次旋转的时候,红色绕\(Y_1\) 旋转,带动外圈的 \(X_1\)也旋转,需要注意的是,经过第一次旋转之后陀螺仪的坐标轴变成了 \(X_2Y_1Z_0(=X_2Y_1z)\),\(Z\) 轴根本没动啊!还和世界坐标系中的 \(z\)轴保持一致;然后再第3次旋转,就只有最外层的绕 \(X_2\)轴旋转了。对于经典配图,也就是下面这张图当中的情况而言,要达到这种状况,就是在做第2次旋转的时候使得\(X_2\)\(z\)轴同向了,这和我们前面解释死锁问题是一致的,这两个轴同向最终导致第3次旋转的时候实际上等价于在旋转\(z\) 轴。

再品一下上面陀螺仪的例子,会发现刚体坐标系和世界坐标系混在一起,而且几个坐标轴还不是同等关系,很容易就把人绕进去了!害人不浅啊!因此要想理解欧拉角万向锁还是去理解前面的内旋\(Z\to Y\to X\) 例子吧。

另一种解释思路是从矩阵运算的角度<spanclass="hint--top hint--rounded" aria-label="关于欧拉角万向锁一直搞不明白万向锁为啥会形成?- 东方白的回答 - 知乎 ">[5],考虑前面提到的旋转矩阵\(\boldsymbol{R}=\boldsymbol{R}_x(\alpha)\boldsymbol{R}_y(\beta)\boldsymbol{R}_z(\gamma)\),当\(\beta=\pm\pi/2\) 的时候,会推导出\[\boldsymbol{R} = \left[\begin{array}{ccc}0 & 0 & 1 \\\sin (\alpha+\gamma) & \cos (\alpha+\gamma) & 0 \\-\cos (\alpha+\gamma) & \sin (\alpha+\gamma) & 0\end{array}\right]\]这里有一点不太合适的,因为万向锁发生在内旋过程中,而这样左乘的矩阵实际上是对应的外旋操作,因此这个解释是有一点小问题的。但也说得通,因为前面提到过内旋和外旋有对应关系,所以可以认为是另一种内旋等价转换为外旋。这里可以看到复合的旋转矩阵\(\boldsymbol{R}\)实际上只做了2个旋转操作,也就是所谓一个自由度被“锁”住了。

前提到内旋和外旋操作是有等价关系的,又提到万向锁只发生在内旋过程中,这是否矛盾呢?我的理解是不矛盾,因为我们说内旋和外旋等价的时候,关注的是结果,也就是两种操作不管过程如何,最终旋转后的坐标系是相同的;而在考虑万向锁问题的时候,我们关注的是旋转的过程。

关于万向锁带来的问题,例如3D动画中的插帧问题,以及云台追踪问题,可以参考博客<spanclass="hint--top hint--rounded" aria-label="欧拉角表示旋转会出现的问题——万向锁(GimbalLock)-CSDN博客 ">[3] 欧拉角表示旋转会出现的问题——万向锁(GimbalLock)-CSDN博客

4. 四元数

其实大家可以发现,欧拉角方式描述三维旋转的一个关键问题是不唯一性,也就是多种不同的旋转操作可以得到相同的结果,被称为Ambiguity。

关于四元数的部分参考大佬写的教程<spanclass="hint--top hint--rounded" aria-label="如何形象地理解四元数?- 知乎 ">[8],后面再写吧。。。

Reference

  1. Eulerangles - Wikipedia↩︎
  2. 欧拉角细节/旋转顺序/内旋外旋- 知乎(zhihu.com)↩︎
  3. 欧拉角表示旋转会出现的问题——万向锁(GimbalLock)-CSDN博客↩︎
  4. 万向锁与欧拉角- 知乎(zhihu.com)↩︎
  5. 关于欧拉角万向锁一直搞不明白万向锁为啥会形成?- 东方白的回答 - 知乎↩︎
  6. 如何通俗地解释欧拉角?之后为何要引入四元数?- 大脸怪的回答 - 知乎↩︎
  7. Extrinsic& intrinsic rotation: Do I multiply from right or left? | by DominicPlein | Medium↩︎
  8. 如何形象地理解四元数?- 知乎 ↩︎
🔲 ⭐

高斯分布常用公式

后面推导中反复用到Woodbury恒等式 \[(A+UCV)^{-1} = A^{-1} - A^{-1}U(C^{-1}+VA^{-1}U)^{-1}VA^{-1}\]

1. 边缘分布

假设 \(\boldsymbol{x}\sim {\mathcalN}(\boldsymbol{\mu}, \boldsymbol\Sigma_1),\boldsymbol{y}|\boldsymbol{x}\sim\mathcal{N}(\boldsymbol{A}\boldsymbol{x},\boldsymbol\Sigma_2)\),那么 \(\boldsymbol{y}\) 的边缘分布也是高斯分布\[\boldsymbol{y} \sim \mathcal{N}(\boldsymbol{\mu}_3, \boldsymbol\Sigma_3)\\\boldsymbol\Sigma_3 = (\boldsymbol\Sigma_2^{-1} -\boldsymbol\Sigma_2^{-1}\boldsymbol{A}(\boldsymbol\Sigma_1^{-1}+\boldsymbol{A}^{\rm T}\boldsymbol\Sigma_2^{-1} \boldsymbol{A})^{-1} \boldsymbol{A}^{\rm T}\boldsymbol\Sigma_2^{-1})^{-1} = \boldsymbol\Sigma_2 +\boldsymbol{A}\boldsymbol\Sigma_1 \boldsymbol{A}^{\rm T} \\\boldsymbol{\mu}_3 = \boldsymbol\Sigma_3 \boldsymbol\Sigma_2^{-\rm T}\boldsymbol{A} (\boldsymbol\Sigma_1^{-1}+\boldsymbol{A}^{\rm T}\boldsymbol\Sigma_2^{-1} \boldsymbol{A})^{-1} \boldsymbol\Sigma_1^{-1}\boldsymbol{\mu} = \boldsymbol{A\mu}\\\]

2. 后验分布

假设有先验 \(\boldsymbol{x}\sim {\mathcalN}(\boldsymbol{\mu}, \boldsymbol\Sigma_1)\),似然 \(\boldsymbol{y}|\boldsymbol{x}\sim\mathcal{N}(\boldsymbol{A}\boldsymbol{x},\boldsymbol\Sigma_2)\),那么关于 \(\boldsymbol{x}\) 的后验分布也是高斯分布\[p(\boldsymbol{x}| \hat{\boldsymbol{y}}) \proptop(\boldsymbol{x},\hat{\boldsymbol{y}}) = p(\boldsymbol{x})p(\hat{\boldsymbol{y}} | \boldsymbol{x}) \sim {\mathcalN}(\boldsymbol{\mu}_4, \boldsymbol\Sigma_4) \\\boldsymbol{\Sigma}_4 = (\boldsymbol\Sigma_1^{-1} + \boldsymbol{A}^{\rmT}\boldsymbol{\Sigma}_2^{-1}\boldsymbol{A})^{-1} \\\boldsymbol{\mu}_4 = \boldsymbol{\Sigma}_4(\boldsymbol{\Sigma}_1^{-1}\boldsymbol{\mu} + \boldsymbol{A}^{\rmT}\boldsymbol{\Sigma}_2^{-1}\hat{\boldsymbol{y}})\] 特殊情况下,假如 \(\boldsymbol{A}=\boldsymbol{I},\boldsymbol{\Sigma}_2=\sigma^2\boldsymbol{I},\boldsymbol{\mu}=\boldsymbol{0}\),此时将有\[\boldsymbol{\Sigma}_4 = \sigma^2 \boldsymbol{\Sigma}_1(\boldsymbol{\Sigma}_1 + \sigma^2\boldsymbol{I})^{-1} \\\boldsymbol\mu_4 =\boldsymbol{\Sigma}_1(\boldsymbol{\Sigma}_1+\sigma^2\boldsymbol{I})^{-1}\hat{\boldsymbol{y}}\]

3. 条件分布

假设 \(\boldsymbol{x}_1,\boldsymbol{x}_2\)服从联合高斯分布 \[\begin{bmatrix} \boldsymbol{x}_1 \\ \boldsymbol{x}_2 \end{bmatrix} \sim{\mathcal N} \left( \begin{bmatrix} \boldsymbol{\mu}_1 \\\boldsymbol{\mu}_2 \end{bmatrix}, \begin{bmatrix}\boldsymbol{\Sigma}_{11} & \boldsymbol{\Sigma}_{12} \\\boldsymbol{\Sigma}_{21} & \boldsymbol{\Sigma}_{22}\end{bmatrix} \right)\] 那么 \(\boldsymbol{x}_2 |\boldsymbol{x}_1\) 也服从高斯分布 \[\boldsymbol{x}_2 | \boldsymbol{x}_1 \sim {\mathcalN}(\boldsymbol{\mu}_5, \boldsymbol{\Sigma}_5) \\\boldsymbol{\Sigma}_5 = \boldsymbol{\Sigma}_{22}- \boldsymbol{\Sigma}_{21} \boldsymbol{\Sigma}_{11}^{-1}\boldsymbol{\Sigma}_{12} \\\boldsymbol{\mu}_5 = \boldsymbol{\mu}_2 + \boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1} (\hat{\boldsymbol{x}}_1 -\boldsymbol{\mu}_1)\]

🔲 ☆

科目三

记录下科目三的主要操作流程与关键点。

0. 细节

  1. 从左至右,离合,刹车,油门;
  2. 前方右侧车道有停止车辆时,需要按两下喇叭,防止他突然开车门;
  3. 直行的时候,驾驶座正前方大约在车道1/3处;
  4. 停车时,离合和刹车要一起踩,若只踩刹车不踩离合容易发动机熄火;
  5. 1挡不能连续行驶超过50m,2挡不能连续行驶超过150m,两者加起来不能超过200m;
  6. 换挡不要求快,慢一点,不要换错;

1. 考试流程

  1. 上车调整座椅,不系安全带,检查灯光回正,观察后视镜、雨刮器的点;
  2. 听到考试指令开始,下车逆时针绕车一周,点一下车身的黑色按钮;
  3. 上车坐好,系安全带;
  4. 灯光考试;
  5. 起步,道路考试;

2. 灯光考试

近光灯:

  1. 夜间在照明良好的条件下行驶。
  2. 夜间在同方向近距离跟车行驶。
  3. 夜间直行通过路口。
  4. 夜间通过有交通信号灯控制的路口。
  5. 夜间与机动车、非机动车会车。

远近光灯交替使用两下:

  1. 夜间超越前方车辆
  2. 夜间通过急弯、坡路、拱桥,人行横道。
  3. 夜间通过没有交通信号灯控制的路口。

大灯关一半开双闪:

  1. 夜间路边临时停车(大灯关一半 示廓灯)双闪 危险警报灯
  2. 夜间在道路上发生交通事故,妨碍交通又难以移动 。

远光灯:

❗夜间在没有路灯、照明不良条件下行驶。

3. 起步

  1. 打转向,按喇叭,挂1挡,松手刹;
  2. 松手刹后的10s内要完成起步;
  3. 观察左后视镜,看后方有没有来车,有的话先慢慢直行,没有的话左转向,慢松离合起步,不要闯动;
  4. 方向盘角度不要太大,到车道一半的时候就可以回正方向盘;
  5. 转向灯恢复正常;

4. 换挡

  1. 踩油门一口气加速到目标速度,然后快速松开;
  2. 离合踩到底;
  3. 换挡;
  4. 松离合,可以稍快一点,不用像科目二那么慢;1、2挡还是要稍微慢一些,3、4挡可以直接松开;

5. 加减挡

  1. 需要逐级加挡到4挡,并保持大于40km/h的速度5s以上;
  2. 然后离合和刹车一起踩到底,当车速降到35时稍微松一点刹车,当车速降到28左右时松开刹车,仍然保持踩着离合,换到3挡;

考试路线上会出现斑马线,因此要在到达斑马线前将车速降到30km/h以下,整体时间还是很紧的,稍微慢一点很可能就不能及时完成。

6. 变更车道

  1. 响起指令以后就打转向灯;
  2. 观察后方有没有车辆,扭头幅度大一些;
  3. 等待3秒,没有车辆就可以换车道;
  4. 结束后记得转向灯回正;

变更车道过程中注意扶好转向灯,不要让他中途掉了。

7. 超车

  1. 先向左变更车道,再向右变更车道;

需要注意的是要在一定距离范围内完成超车,因此要卡好时间,不要拖延。如果后方有来车,可以挂2挡把速度降下来,防止不能在规定时间和距离内完成超车。

8. 路口直行

  1. 距离路口30m内需要减速慢行,如果可以直接通过红绿灯,则距离斑马线20m左右轻踩一下刹车,如果需要等待红灯,则停车挂空挡,然后再起步;
  2. 起步先挂一档,然后踩一下油门提上速度,换二挡;

有时候红绿灯在坡道上,再次起步可能需要半坡起步。

9. 路口转弯

  1. 首先记得打转向灯;
  2. 距离路口30m左右减速到2挡;
  3. 距离路口20m左右踩一下刹车;
  4. 转弯前扭头看一下侧后方有没有车或者行人;
  5. 后面操作跟直行过马路基本相同,但是转弯时必须要在2挡;

考试路线中右转弯基本不需要等红路灯,可以直接转。但是如果刚好遇到左边道路车辆直行过路口,可以等他们过完再右转。

10. 调头

  1. 首先记得打左转向灯;
  2. 首先换到2挡;
  3. 掉头前需要踩一下刹车;
  4. 观察对方有无来车,能调头就直接调,不能的话可以先停车等待,然后再1挡起步调头;
  5. 看到肩膀超过调头的牌子以后打死方向盘就能正常调过去;

11. 会车

  1. 轻踩一下刹车;
  2. 向右打一下方向盘;
  3. 回正;

由于需要在道路中向右打一下方向盘,为防止压实线,可以提前靠左一点行驶。

12. 直线行驶

  1. 就是要保持100m直线行驶,横向偏离幅度不超过30cm,目视前方微调方向盘;

13. 靠边停车

  1. 先打右转向灯;
  2. 观察后方有没有来车,扭头幅度大一些;
  3. 调整到2挡,踩下离合和刹车慢慢走;
  4. 进入最右侧车道,先将车后回正,再慢慢调方向盘使车身不断靠近路边实线;切忌向右摆动幅度过大,一方面不容调整车身回正,另一方面容易压实线;
  5. 观察雨刮器,对齐路边沿白线右侧边缘,使车距离车道线控制在30cm以内;或者观察右侧后视镜,后边车门把手距离白线左侧边缘两指左右;
  6. 停车,挂空档,拉手刹,关闭发动机,松开安全带,观察后方有无来车,3s后打开车门,打开车门后10s内人下车并关闭车门;

14. 整体理解

科目三实际上就是要在实际道路上跑一圈,在路途中完成前面提到的各种操作,不同操作之间基本是相互独立,所以到每个地方只需要考虑对应的操作就可以了,不要觉得东西又多又杂,其实并不杂乱。

贯穿全程的一个关键是速度和挡位的选择,因为1挡和2挡连续行驶距离有最大限制,所以基本上都要保持3挡大约22km/h的速度行驶,采用其他挡位的地方只有:起步(包括中途停车后起步)、路口转弯、调头、靠边停车,以及加到4挡行驶。因此速度和挡位的选择也是很清晰的,主要是注意换挡的时候速度匹配、不要换错挡。

另外最容易错的地方还有:1)忘记打转向灯;2)斑马线/路口前忘记踩刹车;3)换到1、2挡松离合过快造成闯动。

靠边停车要找好点,速度要慢。

因为有的操作并不会有提示(例如转弯前要提前变更车道),因此提前熟悉考试路线很有必要。

🔲 ☆

信号调制与解调

写这篇文章的起因是我一直对数字调制和模拟调制这两个概念比较模糊,对于相应的通信系统结构也觉得比较混乱,所以在网上查了很多资料来理清楚这件事情。本文引用了诸多网络上的博客、文章中的文字和图片,引用之处都标注了参考出处,若有侵权请告知我,我会删除侵权部分。

信号的调制简单来说,就是用基带信号去控制(调制)单频载波信号的某个或某几个参数,从而使得调制后的信号中嵌入了希望传递的信息[1]。这个载波往往是一个高频信号,便于传输。解调则是调制的逆过程。之所以要这么做是因为:1)高频信号更容易收发传输,天线尺寸需要是波长的1/4,使用高频信号可以减小天线尺寸;2)通过调制不同频率的载波信号,可以将基带信号搬移到不同频段的载波上,便于同时传输多路不同的基带信号[2]

1. 移动通信发展历史

我们主要关注无线通信。移动通信最开始的目的是打电话,上世纪80年代出现了第一代(1G)移动通信网络,传输的是人的语音信号,1G使用模拟调制,多址技术为FDMA。当时使用的都是大哥大,采用模拟电路,因此手机都比较大,使用的频段在800MHz,所以天线长度大概要10cm,所以也能看到大哥大外面露出来的很长的天线,这个阶段摩托罗拉风光无限。传输速率只有2.4Kbps,要知道人的语音要想保证比较好的质量需要的速率是64Kbps,因此1G的语音通话质量也不太行[22]

上世纪90年代,形成2G移动通信,大致可以分为两类,一类是基于TDMA的GSM标准,主要使用地区是我国和欧洲;另一类是基于CDMA的IS-95标准(也被称为cdmaOne),主要使用地区是美国。欧洲(芬兰)提出的通信标准GSM中采用了GMSK调制,而IS-95则采用QPSK和OQPSK调制,这三种都是数字调制,因此2G也将移动通信带入了数字时代。由于处理和传输的是数字信号,数字电路取代了模拟电路,能够集成在一个小小的芯片上,因此手机体积大大的减小了,诺基亚在这个阶段大杀四方。并且除了语音信号,其他的多媒体数据例如文字、图片也能传输了。GSM,工作在900MHz和1800MHz频点附近;CDMA工作在800MHz频点附近。2G移动通信,带宽有25MHz,峰值数据速率能达到14.4Kbps[20][21]

在2G和3G中间经过了2.5G,是二者的衔接,2.5G相比于2G提供了更高的速率和更多的功能。蓝牙等技术都是2.5G技术[23]

到3G时代,CDMA技术大行其道,作为CDMA技术的提出者,高通握着大量专利,简直是躺着赚钱。为了绕开高通的专利,欧洲日本联合起来成立了3GPP组织,提出了W-CDMA;中国搞了一个TD-SCDMA;高通与韩国组成3GPP2组织,推出了CDMA2000[24]。前两个可以看作是从GSM演化而来;CDMA2000则是从IS-95继承而来。3G的工作频段在2100MHz附近,带宽有25MHz,上行链路采用BPSK调制,下行链路采用QPSK,数据速率能达到3.1Mbps[19][20]。苹果2007年发布的iPhone也带动了3G的推广和发展。

在进入4G时代之前,插入一段关于WiFi和WiMax的历史。1999年,IEEE分别推出了802.11b与802.11a两种WiFi标准,分别使用2.4GHz和5GHz频段,彼此标准不相容。2003 年,IEEE引入正交频分复用技术(OFDM),推出802.11b的改进版802.11g使传输速度从原先的11Mbps提升至54Mbps。现在我们使用的WiFi主要为802.11n,与802.11a、802.11b、802.11g皆兼容,并采用MIMO技术,使传输速度及距离都有所提升,速度甚至可达600Mbps[24]

随着版图不断扩大,IT业巨头们开始觊觎起蜂窝移动通信市场大饼——4G。OFDM说起来也不是新技术,早在1960年代贝尔实验室发明OFDM后,技术框架约在1980年代便已建立完成。然而当时能支持OFDM的硬件不成熟,CDMA又由高通领军一时红火,便淘汰在3G标准之外。简单来说就是CDMA太红,如果Intel和IT大厂没有在WiFi上将OFDM技术发扬光大,电信业没有一家注意到早期不被重视的OFDM。由于WiMax的关系,OFDM才又重新进入电信业和学术界的视野中[24]

2008年时,3GPP提出了长期演进技术 (Long Term Evolution, LTE)作为3.9G技术标准。又在2011年提出了长期演进技术升级版 (LTE-Advanced)作为4G技术标准,准备把W-CDMA汰换掉,转而采用OFDM。全球覆盖率最高的基站正是W-CDMA,因此,各大运营商无不纷纷决定采用LTE-Advanced当作第四代通信技术标准。如同高通败在W-CDMA基站的广覆盖上,LTE可兼容WCDMA,且利用现有基站配套设备,而WiMax基站却要从头建起,Intel也在2010年宣布放弃WiMax,加入LTE阵营[24]

到4G时代,高通摆烂了,CDMA2000也没有后续向4G的演化,所以4G主要有两支,一个是从W-CDMA演化来的LTE-FDD/LTE-FDD-Advanced,另一个是从TD-SCDMA演化来的TD-LTE/TD-LTE-Advanced。4G采用的OFDM技术可以看成是多址技术,也可以看成是调制技术,实际上就是用多个相互正交的子载波同时传递多路数据,大大提高了频谱效率。FDD-LTE的工作频点在1800MHz附近,TD-LTE的工作频点在1890MHz、2350MHz、2600MHz附近,带宽有100MHz,数据速率可达100Mbps[20][26]

2. 调制的分类

调制信号的大致分类如下图所示[3]

按照载波信号(也被称为被调信号)可以分为三类:正弦波调制、脉冲调制与强度调制。调制的载波分别是正弦波,脉冲和光波。

无线通信中一般使用的载波信号都是高频正弦波,而调制过程中改变的就是正弦波的3个参数:幅度、相位、频率[1]。也就是3种基本的调制方式:调幅(AM)、调相(PM)、调频(FM)。除此之外还有一些变异的调制方法,比如正交幅度调制(QAM),单边带调幅(SM)、残留边带调幅(SSB)等[3]\[{\color{red}A_c} \cos(2\pi {\color{orange}f_c} t +{\color{blue}\varphi_c})\]按照基带信号(也被称为调制信号)可以分为两类:模拟调制与数字调制。顾名思义,模拟调制中基带信号是模拟信号,数字调制中基带信号是数字信号。

模拟调制的中所控制的幅度、频率、相位参数是连续变化的,在解调的过程中也需要估计这个连续变化的波形;而数字调制中这些被改变的参数只是一些离散的值。模拟调制与数字调制各自的优缺点为[2]

优点缺点
数字调制抗干扰能力强;
易于加密,保密性强;
便于计算机对数字信息进行处理;
便于集成化。
需要较宽的频带;
进行数/摸转换时会带来量化误差;
要求的技术和设备复杂。
模拟调制直观且容易实现;保密性差,抗干扰能力差。

调制的另一种分类方法是角度调制和幅度调制,其中角度调制包括调频和调相,幅度调制包括调幅AM、双边带调制DSB、单边带调制SSB、残留边带调制VSB和正交幅度调制QAM[3]

根据已调信号的频谱结构是否保留了原来消息信号的频谱模样,可以分为线性调制与非线性调制[5]。幅度调制一般都是线性调制,角度调制都是非线性调制。

3. 解调的分类

相干解调与非相干解调。相干解调(也被称为同步检波)适用于所有线性调制信号的解调[4]

4. 模拟调制与解调

我们都假设基带信号是 \(m(t)\),载波频率是 \(f_c\),调制后的信号是 \(s(t)\)。本部分主要参考了博客[5]

4.1 常规调幅AM

先将基带信号加上一个直流分量,然后乘以载波得到已调信号 \[s(t) = (A_0 + m(t))\cos(2\pi f_c t)\] 对应的波形和频谱如下图所示:

解调阶段只需要通过一个低通滤波器即可检出信号包络:

4.2 抑制载波的双边带调制DSB-SC

上面的调幅方法由于直流的存在,会导致输出的已调信号中含有载波分量\(\cos(2\pi f_ct)\),抑制载波的双边带调制方法就是去除这个直流分量,也简称为双边带调制。\[s(t) = m(t) \cos(2\pi f_c t)\]

由于DSB已调信号的包络不再与基带信号\(m(t)\)成正比,因此不能使用包络检测器进行解调。通常采用的是相干解调。也就是在接收端生成一个与发送端同频同相的相关载波信号,与接收信号相乘,即可将接收信号频谱再搬到基带,通过低通滤波器后即可解调出基带信号。

调制解调
系统框图
频谱

4.3 单边带调制SSB

前面的双边带调制基带信号是对称的,搬移到载波之后实际上有一半的频谱资源是浪费掉的,没有携带有效的信息,因此单边带调制就是在DSB的基础上,要去掉一半没用的频谱,节省频谱资源。

上面通过低通滤波法获得单边带已调信号对滤波器有较高的要求,需要其在载频处具有陡峭的截止特性。这难以实现,因此实际中往往是通过项移法产生SSB信号,如下图所示,推导过程可以参考博客[5]

4.4 残留边带调制VSB

单边带调制需要陡峭的低通滤波器,为了解决这个问题,残留边带调制的想法虽然理想的陡峭滤波器不可实现,但是圆滑滚降的滤波器还是可以的,那就用这样的滤波器。不过这样不能将其中一半的频谱去除干净,那么是否还能够解调出基带信号呢?也是可以的,不过这个滚降滤波器的频谱需要满足一定性质,也就是在载频处具有对称性[5]

4.5 调相PM与调频FM

可以参考博客[6]。调频信号的包络恒定,而噪声往往是作用在信号幅度上,因此FM信号抗噪声能力强;但是FM信号需要占用较大的信号带宽,频谱利用率低。

5. 数字调制

在开始本节之前,需要明晰的一点是,相比于模拟调制,数字调制当中“调制”的概念已经被扩展了。具体而言,在模拟调制当中,基带模拟信号与高频载波通过前面介绍的几种调制方式混合之后,得到的就直接是待发送的射频信号。而在数字调制中通常包含三个阶段:

  1. 第一个阶段将二进制比特流映射为某种效率更高的数字信号(通常被称为码流),这个过程被称为基带调制;
  2. 第二个阶段将码流通过脉冲成型滤波器,将会突变的数字信号变成连续光滑的模拟信号,得到的就是基带信号;
  3. 第三个阶段将基带信号搬移到高频载波上,这个过程被称为载波调制/带通调制/射频调制。
digital-communication-system-diagram

在本节的几个子小节中,5.3IQ调制属于第三阶段载波调制,5.7脉冲成型滤波器属于第二阶段,其余的部分则属于第一阶段基带调制。为了讲述的方便以及前后逻辑的顺畅,这里将他们穿插起来,但是希望读者在脑海里能够明确他们在系统中的位置以及各自的作用。

实际上,笔者在写本文之前对模拟和数字调制感到困惑、不理解也是源于此。

5.1 二进制数字调制

最简单的数字调制系统中,基带信号是0-1二进制数字波形,通过控制开关实现对载波的调制,因此这里面的调幅AM、调相PM、调频FM也分别被称为幅度键控ASK/通-断键控OOK、相位键控PSK、频移键控FSK。他们各自图示如下[8]

2ASK2PSK2FSK

5.2 多进制数字调制

二进制数字调制中,每一个符号只能表示0-1两个数值,为了提高数据传输效率,可以在一个符号内传输更多的比特,从而提高频带利用率。简单来说就是把原来的0-1比特流进行分组,例如两个bit为一组映射到一个符号(symbol、码元)上,那么一个符号就有00,01, 11, 10这4种取值,再去调制:

  • 幅度:0,1,2,3这4个振幅;
  • 相位:0,\(\pi/2\)\(\pi\),\(3\pi/2\)这4个相位;
  • 频率:4个不同的载波频率。

想要传输的信号从bit流映射为了symbol流,symbol同样是只有有限个离散的取值,这个时候引入星座图会更直观更方便。什么是星座图呢?首先我们传输的正弦信号可以表示为一个复信号的实部,也就是\[{A_c} \cos(2\pi {f_c} t + {\varphi_c}) = \operatorname{Re}\{{\color{blue}A_c \exp(\varphi_c)}\exp (2\pi {f_c} t)\}\] 可以看到无论是调制信号幅度还是相位,都是在单频载波\(\exp(2\pi f_c t)\)上乘以了一个复幅度\(A_c\exp(\varphi_c)\)。因此如果是幅度或者相位调制,那么我们只需要看这个复幅度就得到了想要传输的信息(symbol流)。一个复数可以映射为一个x-y二维平面上的点,那么如果我们把所有的symbol取值画出来,那么就是多个离散的点,就像是星座一样,故名星座图。BPSK(也就是2PSK)、QPSK(也就是4PSK)和8PSK示意图如下[9][16]

其实QPSK还有另外一种方式,就是星座图旋转45°[7]

除了MPSK,幅度也是可以控制的,例如下面这个16APSK[10]

不过这里有一个问题就是:我们实际系统中只能产生实信号,而这里转换到复数域虽然看起来很简洁,但是实际系统如何实现呢?答案就是用两路信号,分别表示复信号的实部和虚部,也就是I路(In-phase)和Q路(Quadrature)。这就引出了下面的正交幅度调制QAM。

5.3 IQ调制

说白了IQ调制就是把复的基带信号搬移到高频载波上,仅此而已。前面我们已经提到了symbol可以用复数表示,那么我们要传输的码流就是一个复基带信号,记为\(s_0(t) = a(t) + j b(t)\),其中 \(a(t)\) 和 \(b(t)\)分别是实部和虚部,那么调制信号可以表示为 \[s(t) = \operatorname{Re}\{s_0(t) \exp(2\pi f_c t)\} = a(t)\cos(2\pi f_ct) - b(t) \sin(2\pi f_c t)\] 在信号发射阶段,I路的基带信号就是\(a(t)\),将其乘以\(\cos(2\pi f_ct)\),Q路的基带信号是\(b(t)\),将其乘以\(\sin(2\pi f_ct)\),再将二者求和即可得到要发送的已调信号。可以参考我之前的一篇博客[11],系统框图如下[12][13]

IQ调制IQ解调

IQ调制在数学上是如此的漂亮,实际系统中又是如此的方便好用,因此在现代无线通信系统中几乎是必备的。后面介绍的各种基带调制方式也基本都会配合IQ调制来使用。

5.4 正交幅度调制QAM

有了IQ调制,我们只需要分别控制I、Q两路基带信号的幅度,就能得到星座图上任意一个点,所以各种MPSK、MASK以及MAPSK等都能通过这种方式获得,这种调制方式也被称为正交幅度调制QAM。

其实也不难发现QPSK实际上就是两路(载波)相互正交的BPSK,实际上也是4QAM。

5.5 偏移四相位键控Offset-QPSK

在一般的QPSK基础上有一些变种,比如OQPSK,与QPSK的区别就在于I、Q两路的码流在时间上错开了半个码元周期,这样可以避免出现两个支路同时出现极性翻转的情况,好处是:1)避免180°的相位跳变,这种情况下最多只有90°相位跳变;2)避免I、Q两路信号同时过零点,降低信号峰均比PAR[7]

除此之外,还有其他的变种,例如\(\pi/4\)-QPSK等,不再赘述。

5.6 最小频移键控MSK

前面介绍的不管是二进制还是多进制,码流信号都认为是矩形波,实际上这种矩形波性能并不好。因为矩形波从时域来看会导致相位跳变,波形不光滑、不连续,从频域来看就是矩形波旁瓣幅度较大,会出现频谱扩展,对相邻信道产生干扰。为了解决这个问题,那就把矩形波换掉!换成边缘光滑的波形!

最小频移键控MSK(Minimum-shiftkeying)可以看成是一种特殊的OQPSK,把原来的矩形方波换成半正弦波。如果从另一个角度来看,也可以认为是频率分离为比特率二分之一的连续相位频移键控CP-BFSK(一般通过开关实现的FSK信号在边沿处也会出现波形不连续,而MSK中波形/相位总是连续的)。如下图所示[15]

那么他的数学原理是什么样的呢?首先从时域来看,实际上就是将OQPSK中的矩形方波换成半正弦波,这样相邻码元的相位变化就为0。如下图所示[16][17],一个脉冲的宽度为\(2T\),正弦波的周期为\(4T\)。

MSK对比O-QPSKMSK的波形是连续的

用三角函数表示出来就是[14] \[s(t) = a_I(t)\cos(\frac{\pi t}{2T}) \cos(2\pi f_c t) -a_Q(t)\sin(\frac{\pi t}{2T})\sin(2\pi f_c t)\] 其中\(a_I(t),a_Q(t)\)还是矩形方波,\(\cos(\pi t/2T),\sin(\pi t/ 2T)\)就是周期为\(4T\)的正弦波,一个码元周期内半个正弦波。对这个式子应用三角恒等式,就能得到\[\begin{aligned}s(t) &= \cos\left( 2\pi f_c t + b_k(t)\frac{\pi t}{2T} + \phi_k\right) \\b_k(t) &= \begin{cases} 1, & a_I(t)=a_Q(t) \\ -1, &\text{others} \end{cases} \\\phi_k &= \begin{cases} 0, & a_I(t)=1 \\ \pi, &\text{others} \end{cases}\end{aligned}\]从频域来看,MSK是频率分离为比特率二分之一的连续相位频移键控CPFSK,也就是在一般的BFSK基础上,MSK的两个载波频率间隔恰好是码元波特率的二分之一。假设FSK的两个载波频率为\(f_1,f_2\),码元周期为\(2T\),由于一个码元可以传输2个比特,因此波特率为\(1/T\),那么就有\(|f_1-f_2| = 1/2T\)。

从频谱来看,MSK的旁瓣幅度相比于BPSK和QPSK也更低,因此, MSK情况下的信道间干扰较低[14]

发射接收机如下图所示[16],解调也需要相干载波。

发射机接收机

5.7 脉冲成型滤波器

MSK当中将矩形波换成了半正弦波,在实现的时候直接用一个正弦波发生器加上开关就可以了,但如果我们想换成其他形状的光滑波形呢?就没有这么简单了,这里就引出了数字传输中的脉冲成型滤波器。

数字调制之后得到的码流是矩形波,从时域来看存在跳变,从频域来看频谱宽度无穷大,这在实际系统当中是不可能实现的,原因是我们无法生成跳变的信号,并且信道和收发机的通带也不是无穷带宽的。因此我们需要将码流脉冲转化为模拟信号,这个模拟信号从时域来看是连续光滑的,从频域来看是带限的,如下图所示[29]。思路很简单,将码流脉冲通过一个脉冲成型滤波器,输出信号就是滤波器冲激响应的叠加。不过关键的问题在于这个脉冲成型滤波器要怎么设计?

设计的关键在于不能有符号间串扰(Inter Symbol Interference,ISI)。什么意思呢?频域带限会导致时域扩展,要想当前码元对相邻的码元没有干扰,如下图所示[30],需要其时域波形在相邻码元的采样时刻幅值为0,这就是Nyquist准则。

记经过基带调制的复码流为 \(\{a_k\}\),码元间隔(采样周期)为 \(T_s\),脉冲成型滤波器的冲激响应为 \(g(t)\),那么脉冲成型之后的模拟基带信号为\[s_0(t) = \sum_{k=-\infty}^{\infty} a_k g(t-kT_s)\] 对应的Nyquist无码间串扰条件即为: \[g(kT_s) = \begin{cases}1, & k=0 \\ 0, & k\ne 0 \end{cases}\] 等价的频域条件为(\(G(f)\)为\(g(t)\)的傅里叶变换)[31]\[\frac{1}{T_s} \sum_{k=-\infty}^{\infty} G\left(f-\frac{k}{T_S}\right) =1, \forall f\] 实际系统中常用的就是升余弦滤波器(Raised CosineFilter),不同滚降系数的RC滤波器频谱如下图所示[31],他两侧的边沿是余弦曲线。

除此之外还有Sinc filter(频谱是理想矩形),Gaussianfilter(时域冲激响应是高斯的PDF)。第二代移动通信标准GSM中就采用高斯滤波器。

5.8 高斯最小频移键控GMSK

终于到了GMSK,其被应用于GSM通信标准、蓝牙、卫星通信当中。GMSK是在MSK的基础上,为了改善信号的旁瓣衰减性能,增加了一个高斯滤波器进行整形[14]。其实现方式可以是首先将基带码流通过高斯滤波器得到高斯脉冲,然后再通过MSK调制器。

系统框图如下所示<spanclass="hint--top hint--rounded" aria-label="Turletti, Thierry. "GMSK ina nutshell." Telemedia Networks and Systems Group LCS, MIT-TR(1996).">[18]

下面的图片当中<spanclass="hint--top hint--rounded" aria-label="Turletti, Thierry. "GMSK ina nutshell." Telemedia Networks and Systems Group LCS, MIT-TR(1996).">[18],第一行是已调信号相位轨迹随着非归零码序列(-1,-1,-1,+1,+1,-1,+1,+1,+1,+1,-1,+1,-1,+1,-1,-1)的变化,可以看到GMSK输出的信号相位变化更加平缓。

5.9 正交频分复用OFDM

OFDM本身又是一个非常精妙的调制技术。简单来理解,就是同时传输多路IQ调制的信号,不同路的载波频率不同,在一个符号周期内他们相互正交,互不干扰。

再详细一点的解释。首先对于IQ调制来说,I、Q两路本身是相互正交的,但是他们的载波还是位于同一个频点上,把两路放在一起看就是把基带信号的复频谱搬移到单一载波频点上[11]。而OFDM可以理解为在N个不同的载波频点上,并行传输了N个这样的IQ调制信号。需要注意这些子载波也不是随意选取的,需要满足在一个符号周期内的正交性,子载波之间不能有干扰。示意图如下[27]

但是实际系统中直接这样实现会有很多麻烦,所以人们想出了一种非常巧妙的办法。要传输的串行二进制比特流首先转化为N路并行的比特流,每一路比特流在经过QAM映射为码流,这些码流控制了不同频点载波的复幅度,所以同一时刻不同支路的码元,实际上就是要发送信号的频谱,那么对这些系数进行IFFT,我们就把他们从频域变换到了时域,N路并行的码元映射为了1路串行的N个时刻的(复)信号,然后再将得到时域波形实部和虚部送入一个IQ调制,即可得到要发送的射频信号。示意图如下[27]

OFDM的系统框图如下所示[28]

除此之外,OFDM中还包括了循环前缀以及其他细节的设计,在这里不再详细展开了,可以参考网上相关资料。LTE中OFDM采用的调制方式最高能达到256QAM,NR中能达到1024QAM[25]

6. 数字解调

解调是调制的逆过程,所以数字解调也包括三个阶段:

  1. 第一阶段将信号从载波频段搬回到基带,如果调制的时候采用的是IQ调制,那么解调就采用IQ解调就好了;
  2. 第二阶段从模拟的基带波形中恢复出数字信号,最佳接收方案为匹配滤波并以\(1/T_s\)的频率进行采样;
  3. 第三阶段从数字信号中解调出原本想传输的比特流,需要根据星座图对得到的采样进行判决。

7. 总结

QAM和IQ调制yyds!OFDM yyds!

Reference

  1. https://www.ni.com/zh-cn/innovations/white-papers/06/analog-and-digital-modulation.html↩︎
  2. https://blog.csdn.net/SilenceBurster/article/details/53126468↩︎
  3. Modulation- Wikipedia↩︎
  4. https://blog.csdn.net/m0_51288996/article/details/121340840↩︎
  5. https://blog.csdn.net/weixin_50912862/article/details/114679288↩︎
  6. https://blog.csdn.net/weixin_50912862/article/details/114694510↩︎
  7. https://blog.csdn.net/weixin_50912862/article/details/114735043↩︎
  8. https://www.eecs.yorku.ca/course_archive/2010-11/F/3213/CSE3213_07_ShiftKeying_F2010.pdf↩︎
  9. https://zhuanlan.zhihu.com/p/58119209↩︎
  10. https://en.wikipedia.org/wiki/Amplitude_and_phase-shift_keying↩︎
  11. https://glooow1024.github.io/2021/08/29/communication/IQ-modulation/↩︎
  12. https://www.allaboutcircuits.com/textbook/radio-frequency-analysis-design/radio-frequency-demodulation/understanding-i-q-signals-and-quadrature-modulation/↩︎
  13. https://www.allaboutcircuits.com/textbook/radio-frequency-analysis-design/radio-frequency-demodulation/understanding-quadrature-demodulation/↩︎
  14. https://en.wikipedia.org/wiki/Minimum-shift_keying↩︎
  15. MinimumShift Keying (MSK) - A Tutorial - Qasim Chaudhari↩︎
  16. http://www.dsplog.com/2009/06/16/msk-transmitter-receiver/↩︎
  17. https://ppt-online.org/1039280↩︎
  18. Turletti, Thierry. "GMSK ina nutshell." Telemedia Networks and Systems Group LCS, MIT-TR(1996).↩︎
  19. https://en.wikipedia.org/wiki/Code-division_multiple_access↩︎
  20. https://www.zseries.in/telecom%20lab/telecom%20generations/↩︎
  21. https://baike.baidu.com/item/CDMAOne/7525860↩︎
  22. https://blog.csdn.net/ziv669/article/details/122453973↩︎
  23. 2.5G_百度百科(baidu.com)↩︎
  24. https://blog.csdn.net/zuochao_2013/article/details/78337600↩︎
  25. https://blog.csdn.net/m0_52840978/article/details/123674928↩︎
  26. https://norstc.blog.csdn.net/article/details/80504263↩︎
  27. https://blog.csdn.net/madongchunqiu/article/details/18614233↩︎
  28. https://ecse.monash.edu/staff/eviterbo/OTFS-VTC18/Tutorial_ICC2019___OTFS_modulation.pdf↩︎
  29. https://wirelesspi.com/pulse-shaping-filter/↩︎
  30. https://dsp.stackexchange.com/questions/36340/nyquist-criterion-for-zero-isi↩︎
  31. https://en.wikipedia.org/wiki/Nyquist_ISI_criterion↩︎
🔲 ⭐

IEEE论文爬虫及数据统计

1. IEEE论文爬虫

爬虫代码网上有很多了,这部分是直接用的网上可以跑通的[1]。使用的时候直接调用get_article_info(),其中参数 conferenceID需要手动在 IEEE 上查询会议的 ID 号,参数 saceFileName为希望保存的 csv 文件名。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
# 获取issueNumber
def get_issueNumber(conferenceID):
"""
Get the issueNumber from the website.
"""
conferenceID = str(conferenceID)
gheaders = {
'Referer': 'https://ieeexplore.ieee.org/xpl/conhome/'+conferenceID+'/proceeding',
'User-Agent': 'Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/83.0.4103.97 Safari/537.36'
}
md_url = 'https://ieeexplore.ieee.org/rest/publication/home/metadata?pubid='+conferenceID
md_res = requests.get(md_url, headers = gheaders)
md_dic = json.loads(md_res.text)
issueNumber = str(md_dic['currentIssue']['issueNumber'])
return issueNumber

# 爬取论文及其下载链接
def get_article_info(conferenceID, saveFileName):
"""
Collect the published paper data, and save into the csv file "saveFileName".
"""
# 获取issueNumber
issueNumber = str(get_issueNumber(conferenceID))
conferenceID = str(conferenceID)

# 记录论文数据
dataframe = pd.DataFrame({})
paper_title = []
paper_author = []
paper_year = []
paper_citation = []
paper_abstract = []
paper_ieee_kwd = []

# 从第一页开始下载
pageNumber = 1
count = 0
while(True):
# 获取会议文章目录
toc_url = 'https://ieeexplore.ieee.org/rest/search/pub/'+conferenceID+'/issue/'+issueNumber+'/toc'
payload = '{"pageNumber":'+str(pageNumber)+',"punumber":"'+conferenceID+'","isnumber":'+issueNumber+'}'
headers = {
'Host': 'ieeexplore.ieee.org',
'User-Agent':'Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/88.0.4324.190 Safari/537.36',
'Referer': 'https://ieeexplore.ieee.org/xpl/conhome/'+conferenceID+'/proceeding?pageNumber='+str(pageNumber),
}
toc_res = requests.post(toc_url, headers = headers, data=payload)
toc_dic = json.loads(toc_res.text)
try:
articles = toc_dic['records']
except KeyError:
break
else:
for article in articles:
title = article['highlightedTitle']
paper_link = IEEE_root_url + article['htmlLink']
paper_info = requests.get(url=paper_link, headers=headers, timeout=10)
soup = BeautifulSoup(paper_info.text, 'lxml') # 解析
# 正则表达式 创建模式对象
pattern = re.compile(r'xplGlobal.document.metadata=(.*?)"};', re.MULTILINE | re.DOTALL)
script = soup.find("script", text=pattern) # 根据模式对象进行搜索
try:
res_dic = pattern.search(script.string).group(1)+'"}' # 配合search找到字典,匹配结尾字符串,降低文章摘要中也出现这种字符串的概率
# 解析异常,一般是因为文章 abstract 中出现了字符串 '"};'
json_data = json.loads(res_dic) # 将json格式数据转换为字典
except Exception as e:
print(pattern.search(script.string).group(0))
print(res_dic)
# 保存文章信息
paper_title.append(title)
paper_year.append(json_data['publicationYear'])
print(json_data.keys())
#a = input('input anything...')
if 'author' in json_data.keys():
paper_author.append(json_data['author'])
else:
paper_author.append(None)
if 'abstract' in json_data.keys():
paper_abstract.append(json_data['abstract'])
else:
paper_abstract.append(None)
if 'keywords' in json_data.keys():
paper_ieee_kwd.append(json_data['keywords'][0]['kwd']) # ieee有三种 key words
else:
paper_ieee_kwd.append(None)
count=count+1
#link = 'https://ieeexplore.ieee.org/stampPDF/getPDF.jsp?tp=&arnumber='+article['articleNumber']+'&ref='
#alf.write(title.replace('\n','')+'>_<'+link+'\n')

# 写入csv文件
dataframe = pd.DataFrame({'title':paper_title, 'year':paper_year, 'abstract':paper_abstract, 'key words':paper_ieee_kwd})
dataframe.to_csv(saveFileName, index=True, sep=',')
print('Page ', pageNumber, ', total ', count, 'papers.')
pageNumber = pageNumber+1
# 停一下防禁ip
import time
time.sleep(3)

# 写入csv文件
dataframe = pd.DataFrame({'title':paper_title, 'year':paper_year, 'abstract':paper_abstract, 'key words':paper_ieee_kwd})
dataframe.to_csv(saveFileName, index=True, sep=',')
return

2. IEEE论文数据统计

3. 写一个图形界面

3.1 弹出提示窗口

在写代码过程中有时候需要测试功能是否成功实现,于是想要加一个弹出窗口的函数可以显示调试信息,用以验证想要的功能是否正常实现。主要难点在于根据内容自动调整窗口大小,以获得较好的显示效果。

采用的方法是利用 QLabel.adjust()函数获取文本显示的宽度,并据此调整窗口的大小[2]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
from PyQt6.QtWidgets import (QWidget, QDialog, QLabel, QPushButton)
from PyQt6.QtCore import (QSize, QRect)

class PaperCollector(QWidget):
def __init__(self):
super().__init__()
self.initUI()

def initUI(self):
self.dialog_btn = QPushButton('Click')
self.dialog_btn.clicked.connect(self.click_callback)
self.setGeometry(300, 300, 300, 200)
self.setWindowTitle('IEEE paper collector (by Glooow)')
self.show()

def click_callback(self):
self.show_dialog('You clicked me!')

def show_dialog(self, info):
"""
Pop up dialogs for debug.
"""
hint_dialog = QDialog()
hint_dialog.setWindowTitle('Hint info')
#hint_dialog.setWindowModality(PyQt6.QtCore.Qt.NonModal)

hint_info = QLabel(info, hint_dialog)
hint_info.adjustSize()
padding = 20
max_width = 360
# set the maximum width
if hint_info.size().width() > max_width:
hint_info.setGeometry(QRect(0, 0, max_width, 80))
hint_info.setWordWrap(True)
hint_info.move(padding, padding)

hint_dialog.resize(hint_info.size() + QSize(padding*2, padding*2))
hint_dialog.exec()

3.2 文本框显示爬取日志

我希望在窗口中增加一个文本框,将爬取过程中的日志信息打印出来,便于用户实时监测。

采用的思路是定义一个logging.Logger,将其日志信息同时输出到窗口的文本框和控制台中打印,通过自定义logging.Handler可以实现这一功能[3][5][6]。实现方式为:

  1. 继承 logging.Handler类,并初始化阶段将整个窗口(QWidget类)作为参数传入,便于后续修改窗口的信息;

  2. 自定义实现 emit 函数,在 emit 函数中将log 信息同时输出到窗口文本框、打印到控制台;

  3. 创建 logger 的时候设置Handler[4]

    1
    2
    3
    4
    ex = PaperCollector()
    logger = logging.getLogger("logger")
    handler = LogHandler(ex)
    logger.addHandler(handler)

下面是这部分功能相关的代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
import logging

class LogHandler(logging.Handler):
def __init__(self, parent):
super().__init__()
self.parent = parent

def emit(self, record):
try:
print(self.format(record))
self.parent.print_log(self.format(record))
QApplication.processEvents()
except Exception:
self.handleError(record)

class PaperCollector(QWidget):
def __init__(self):
super().__init__()
self.initUI()

def initUI(self):
"""
Define the UI playout.
"""
# button to start crawing
self.startCrawling_button = QPushButton('Start')
self.startCrawling_button.setToolTip('Click and wait for collecting published paper data.')
self.startCrawling_button.clicked.connect(self.start_collect_paper)
# print log
self.process = QTextEdit(readOnly=True)
self.process.setFont(QFont("Source Code Pro",9))

grid = QGridLayout()
grid.setSpacing(10)
grid.addWidget(self.startCrawling_button, 1, 0)
grid.addWidget(self.process, 2, 0, 3, 3)
self.setLayout(grid)

self.setGeometry(300, 300, 700, 300)
self.setWindowTitle('IEEE paper collector (by Glooow)')
self.show()

def start_collect_paper(self):
global logger
#self.show_dialog('start!')
get_article_info(self.conferenceID_edit.text(), self.saveFile_edit.text(), logger)

def print_log(self, s):
self.process.append(s)

logger = None

def main():
app = QApplication(sys.argv)
ex = PaperCollector()

global logger
logger = logging.getLogger("logger")
logger.setLevel(logging.INFO)
formater = logging.Formatter(fmt="%(asctime)s [%(levelname)s] : %(message)s"
,datefmt="%Y/%m/%d %H:%M:%S")
handler = LogHandler(ex)
handler.setFormatter(formater)
logger.addHandler(handler)

sys.exit(app.exec())


if __name__ == '__main__':
main()

爬取论文的主函数如下,其中一个参数为logger,在函数内部需要打印日志信息的地方添加logger.info(...) 即可。

1
2
def get_article_info(conferenceID, saveFileName, logger):
logger.info('collecting paper......')

3.3 多线程避免卡顿

上述打印日志的方法不能做到实时输出信息到窗口文本框,而是会等到所有论文爬取完毕之后再一股脑的更新,这是因为PyQt的界面线程是主线程,当爬虫开始工作时,也是运行在主线程中,这时主界面就无法更新,看起来就像是卡死了。解决方法就是开一个子线程运行爬虫工作<spanclass="hint--top hint--rounded" aria-label="PyQt- 使用多线程避免界面卡顿 - bailang zhizun的博客 - CSDN博客">[7]

具体实现细节为:

  1. 新建类 SpiderThread 继承 QObject,自定义run 函数,在其中运行爬虫程序;
  2. SpiderThread 类中定义一个_spider_finish = pyqtSignal(),该信号用于告知主线程爬虫子线程已完成工作
  3. PaperCollector 类中定义一个_start_spider = pyqtSignal(str, str, logging.Logger),该信号用于启动爬虫子线程[8][9]
  4. 通过 pyqtSignal.connect分别将各个信号连接到对应的槽(处理函数)上;

下面是这部分功能相关的代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
from PyQt6.QtCore import (QObject, pyqtSignal, QThread)

class SpiderThread(QObject):
_spider_finish = pyqtSignal()

def __init__(self):
super().__init__()
self.flag_running = False

def __del__(self):
print('>>> __del__')

def run(self, conference_ID, save_filename, logger):
get_article_info(conference_ID, save_filename, logger)
self._spider_finish.emit()

class PaperCollector(QWidget):
_start_spider = pyqtSignal(str, str, logging.Logger)

def __init__(self):
super().__init__()
self.initUI()
#sys.stdout = LogStream(newText=self.onUpdateText)

self.spiderT = SpiderThread()
self.thread = QThread(self)
self.spiderT.moveToThread(self.thread)
self._start_spider.connect(self.spiderT.run) # 只能通过信号槽启动线程处理函数
self.spiderT._spider_finish.connect(self.finish_collect_paper)

def start_collect_paper(self):
if self.thread.isRunning():
return

self.startCrawling_button.setEnabled(False)
self.startCrawling_button.setToolTip('I\'m trying very hard to collect papers >_<')
# 先启动QThread子线程
self.thread.start()
# 发送信号,启动线程处理函数
# 不能直接调用,否则会导致线程处理函数和主线程是在同一个线程,同样操作不了主界面
global logger
self._start_spider.emit(self.conferenceID_edit.text(), self.saveFile_edit.text(), logger)

def finish_collect_paper(self):
self.startCrawling_button.setEnabled(True)
self.startCrawling_button.setToolTip('Click and wait for collecting published paper data ^o^')
self.thread.quit()

def stop_collect_paper(self):
if not self.thread.isRunning():
return
self.thread.quit() # 退出
self.thread.wait() # 回收资源
self.show_dialog('stop!')

3.4 流畅中止子线程

有时候我们需要中途停止爬虫工作,比如发现会议ID设置错误、希望先对已经爬取的部分数据进行统计分析等。在上面的实现中,尽管线程正常运行很流畅,但是如果在爬虫运行中途点击停止按钮,程序就会卡死。

在原本的爬虫脚本中,get_article_info()函数内部的爬虫采用了 while(True) 死循环,主线程中直接用self.thread.quit()强制退出,从控制台来看这样确实可以停掉,但是Qt窗口却总是会卡死。原因我也不太清楚,采用的解决方法是:

  1. 定义一个爬虫类 IEEESpider,设置成员变量flag_running,将函数 get_article_info也设置为类成员函数;
  2. get_article_info 中的循环改为while(self.flag_running)
  3. 在主线程中想要停止爬虫子线程的时候,只需要首先设置flag_running=False,那么爬虫子线程在当前一次循环结束后就自动结束,这个时候主线程调用self.thread.quit() 就不会导致界面卡死。需要注意的是设置flag_running=False 一定要 sleep一段时间,以保证爬虫子线程能够结束当前循环,否则还是容易卡死。

下面是这部分功能的代码。

1
2
3
4
5
6
7
class IEEESpider:
def __init__(self):
self.flag_running = False

def get_article_info(self, conferenceID, saveFileName, logger):
while(self.flag_running):
pass
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
class SpiderThread(QObject):
def __init__(self):
super().__init__()
#self.flag_running = False
self.ieee_spider = IEEESpider()

def run(self, conference_ID, save_filename, logger):
self.ieee_spider.flag_running = True
self.ieee_spider.get_article_info(conference_ID, save_filename, logger)
self._spider_finish.emit()

class PaperCollector(QWidget):
def stop_collect_paper(self):
if not self.thread.isRunning():
return
self.spiderT.ieee_spider.flag_running = False
time.sleep(15)
self.thread.quit() # 退出
#self.thread.wait() # 回收资源
self.show_dialog('stop!')

3.5 增加侧边导航栏

前面只有爬取论文的页面,现在我想加上数据分析的页面,那么就需要设置一个侧边导航栏,以切换两种不同的任务。

实现方式为左侧设置多个按钮,右侧添加一个QTabWidget(),将不同的页面设置为子标签页,通过按钮的点击回调函数切换不同的标签页[10]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
class PaperCollector(QWidget):
def sidebarUI(self):
"""
Define the UI playout of sidebar.
"""
self.sidebar_btn_1 = QPushButton('Collector', self)
self.sidebar_btn_1.clicked.connect(self.sidebar_button_1)
self.sidebar_btn_2 = QPushButton('Analyzer', self)
self.sidebar_btn_2.clicked.connect(self.sidebar_button_2)
self.sidebar_btn_3 = QPushButton('Reserved', self)
self.sidebar_btn_3.clicked.connect(self.sidebar_button_3)

sidebar_layout = QVBoxLayout()
sidebar_layout.addWidget(self.sidebar_btn_1)
sidebar_layout.addWidget(self.sidebar_btn_2)
sidebar_layout.addWidget(self.sidebar_btn_3)
sidebar_layout.addStretch(5)
sidebar_layout.setSpacing(20)

self.sidebar_widget = QWidget()
self.sidebar_widget.setLayout(sidebar_layout)

def sidebar_button_1(self):
self.right_widget.setCurrentIndex(0)

def sidebar_button_2(self):
self.right_widget.setCurrentIndex(1)

def sidebar_button_3(self):
self.right_widget.setCurrentIndex(2)

def initUI(self):
"""
Define the overall UI playout.
"""
self.sidebarUI()
self.spiderUI()
self.analyzerUI()
self.reservedUI()

# 多个标签页
self.right_widget = QTabWidget()
self.right_widget.tabBar().setObjectName("mainTab")

self.right_widget.addTab(self.spider_widget, '')
self.right_widget.addTab(self.analyzer_widget, '')
self.right_widget.addTab(self.reserved_widget, '')

# 隐藏标签部件的标签并初始化显示页面
self.right_widget.setCurrentIndex(0)
self.right_widget.setStyleSheet('''QTabBar::tab{width: 0; height: 0; margin: 0; padding: 0; border: none;}''')

# overall layout
main_layout = QHBoxLayout()
main_layout.addWidget(self.sidebar_widget)
main_layout.addWidget(self.right_widget)
main_layout.setStretch(0, 40)
main_layout.setStretch(1, 200)
self.setLayout(main_layout)

self.setGeometry(300, 300, 850, 300)
self.setWindowTitle('IEEE paper collector (by Glooow)')
self.show()

3.6 next ...

接下来考虑:写数据分析页面 ......

Referencce

这里关于参考文献的部分,本来我想按照下面格式来写,希望实现的效果是都像[2][10]一样,每一条引用列出来的是超链接,而不是直接写出来链接地址,但是我发现除了第[2][10]条,其他条这么写话都会像现在的第<spanclass="hint--top hint--rounded" aria-label="PyQt- 使用多线程避免界面卡顿 - bailang zhizun的博客 - CSDN博客">[7]条一样,格式会乱,也不知道为什么。有人知道的话可以告诉我嘛>_<

1
2
3
4
5
6
7
8
9
10
[^1]:[Python爬虫——爬取IEEE论文 - 乐 ShareLe的博客 - CSDN博客](https://blog.csdn.net/wp7xtj98/article/details/112711465)
[^2]:[PyQt 中文教程 (gitbook.io)](https://maicss.gitbook.io/pyqt-chinese-tutoral/)
[^3]:[python日志:logging模块使用 - 知乎](https://zhuanlan.zhihu.com/p/360306588)
[^4]:[python3 自定义logging.Handler, Formatter, Filter模块 - 太阳花的小绿豆的博客 - CSDN博客](https://blog.csdn.net/qq_37541097/article/details/108317762)
[^5]:[python logging output on both GUI and console - stackoverflow](https://stackoverflow.com/questions/41176319/python-logging-output-on-both-gui-and-console)
[^6]:[How to dynamically update QTextEdit - stackoverflow](https://stackoverflow.com/questions/24371274/how-to-dynamically-update-qtextedit)
[^7]:[PyQt - 使用多线程避免界面卡顿 - bailang zhizun的博客 - CSDN博客](https://blog.csdn.net/bailang_zhizun/article/details/109240670)
[^8]:[pyqt 带单个参数/多个参数信号&槽总结 - gong xufei的博客 - CSDN博客](https://blog.csdn.net/gong_xufei/article/details/89786272)
[^9]:[PyQt5 pyqtSignal: 自定义信号传入的参数方法 - Mic28的博客 - CSDN博客](https://blog.csdn.net/qq_39560620/article/details/105711799)
[^10]:[PyQt5 侧边栏布局 • Chang Luo (luochang.ink)](https://www.luochang.ink/posts/pyqt5_layout_sidebar/)
  1. https://blog.csdn.net/wp7xtj98/article/details/112711465↩︎
  2. PyQt中文教程 (gitbook.io)↩︎
  3. https://zhuanlan.gitbook.io/p/360306588↩︎
  4. https://blog.csdn.net/qq_37541097/article/details/108317762↩︎
  5. https://stackoverflow.com/questions/41176319/python-logging-output-on-both-gui-and-console↩︎
  6. https://stackoverflow.com/questions/24371274/how-to-dynamically-update-qtextedit↩︎
  7. PyQt- 使用多线程避免界面卡顿 - bailang zhizun的博客 - CSDN博客↩︎
  8. https://blog.csdn.net/gong_xufei/article/details/89786272↩︎
  9. https://blog.csdn.net/qq_39560620/article/details/105711799↩︎
  10. PyQt5侧边栏布局 • Chang Luo (luochang.ink)↩︎
🔲 ☆

SLAM综述

关于SLAM(simultaneous localization andmapping)问题一个很粗糙的总结。

对于SLAM的发展历程,Leonard和Reid大佬将SLAM到目前为止的发展过程总结为三个阶段:

  1. classicalage(1986-2004):早期阶段,SLAM问题的定义、基于概率框架的建模和求解方法;
  2. algorithm-analysisage(2004-2015):深入研究SLAM问题的一些性质,比如稀疏性、收敛性、一致性等,更多样、更高效的算法也被相继提出;
  3. robust-perceptionage(2015-):开始考虑算法的鲁棒性、可扩展性、资源约束下的高效算法、高层语义认知任务导向等;

SALM系统主要包括前端和后端两部分,前端负责从传感器数据中提取特征、dataassociation(如特征提取、回环检测)等,后端负责最大后验估计、滤波等;

根据感知手段可以分为几类:

  • 测距测角:传感器比如毫米波雷达、声纳等;
  • 视觉相机:各种相机,比如RGB-D相机等;需要从图像中提取特征点,例如SIFT等;纯视觉导航有专门的研究方向,即VO(VisualOdometry),VO+全局地图优化(例如回环检测)=visual SLAM;
  • 激光雷达;
  • 惯性导航:IMU,配合视觉传感器可以实现VIO;

根据地图结构可以分为几类(地图结构与感知手段紧密相关):

  • 基于landmarkd的,传感器通常是测距测角的,或者视觉中提取特征点;
  • 栅格地图(2D)、点云地图(3D),主要是激光雷达;
  • 基于边/表面的几何地图;

根据求解方法可以分为几类:

  • 基于贝叶斯框架递归滤波的:如EKF、PF、Information Filter等;
  • 基于优化的:一种是光束平差法(用于视觉,实际上是最小二乘);另一种是图优化(graphSLAM,实际上就是概率图模型,似乎是现在的主流);
  • Set Membership 类的方法、以及定性方法等;

要考虑的科学问题包括:

  • 地图构建和机器人位姿估计,因此实际上是参数估计/求解问题;
  • 多传感器数据融合,以及与已有地图信息的融合;
  • 回环检测,纠正累积误差;
  • 降低复杂度,比如每次只更新相关的部分地图和状态,而不是所有参数一起更新;
  • 数据关联;

除了基本的SLAM问题,还有一些发展方向:

  • ActiveSLAM:同时考虑机器人的运动控制,与exploration-exploitation问题相关;
  • 多机SLAM,相比于单机需要考虑的问题:
    • 多机协作的架构:集中式(online or offline?)、分布式;
    • 相邻机器之间的数据交互,防止数据伦理问题,即一个数据在多个机器上重复使用多次;
    • 每个机器的地图更新方法;
    • 数据关联问题;
    • 通信负载、数据压缩、地图结构优化等;

Reference

  1. Dissanayake, MWM Gamini, et al. "A solution to the simultaneouslocalization and map building (SLAM) problem." IEEE Transactions onrobotics and automation 17.3 (2001): 229-241.
  2. Durrant-Whyte, Hugh, and Tim Bailey. "Simultaneous localization andmapping: part I." IEEE robotics & automation magazine 13.2(2006): 99-110.
  3. Bailey, Tim, and Hugh Durrant-Whyte. "Simultaneous localization andmapping (SLAM): Part II." IEEE robotics & automationmagazine 13.3 (2006): 108-117.
  4. Bresson, Guillaume, et al. "Simultaneous localization and mapping: Asurvey of current trends in autonomous driving." IEEE Transactionson Intelligent Vehicles 2.3 (2017): 194-220.
  5. Taketomi, Takafumi, Hideaki Uchiyama, and Sei Ikeda. "Visual SLAMalgorithms: A survey from 2010 to 2016." IPSJ Transactions onComputer Vision and Applications 9.1 (2017): 1-11.
  6. Cadena, Cesar, et al. "Past, present, and future of simultaneouslocalization and mapping: Toward the robust-perception age." IEEETransactions on robotics 32.6 (2016): 1309-1332.
🔲 ☆

【随机过程2】连续参数马尔可夫链

7.1 定义与基本概念

定义:设随机过程 \(X=\{X(t),t\ge0\}\),状态空间 \(S\),对任意 \(0\le t_0 < t_1 < \cdots < t_n <t_{n+1}\),\(i_k\in S\),若\(P(X(T_k)=i_k,0\le k\le n) > 0\)\[P(X(t_{n+1})=i_{n+1} | X(t_k)=i_k,0\le k\le n) = P(X(t_{n+1})=i_{n+1} |X(t_n)=i_n)\] 则称 \(X\)连续参数的马氏链。若对任意的 \(s,t\ge0\),\(i,j\in S\) 有 \[P(X(s+t)=j | X(s)=i) = P(X(t)=j|X(0)=i) = P_{ij}(t)\]\(X\)齐次马氏链

对于连续参数马氏链,在原点处有 \(P_{ij}(0)=\delta_{ij}\),因此 \(P(0)=I\)。除此之外假设其在原点处连续,即\(\lim_{t\to 0}P_{ij}(t)=\delta_{ij},\lim_{t\to0}P(t)=I\)。类比离散参数的马氏链,可以得到类似的C-K 方程 \(P(s+t)=P(s)P(t)\)

7.2 转移概率矩阵

根据 C-K 方程可以猜测 \(P(t)=e^{tQ}\),泰勒展开表示为 \(P(t)=I + \sum_{n=1}^{\infty}\frac{t^n}{n!}Q^n\),\(P(t)\)完全由 \(Q\) 确定,并且有 \(P'(0)=\lim_{t\to0}\frac{P(t)-I}{t}=Q\)。

上面只是猜测,那么是否真的存在这样一个 \(Q\) 呢?下面两个定理给出结论。

定理 7.1:对 \(i\inS\),极限 \(-q_{ii}=\lim_{t\to0}\frac{1-P_{ii}(t)}{t}\)存在,但可能是无穷。

定理 7.2:对 \(i\inS\),极限 \(q_{ij}=\lim_{t\to0}\frac{P_{ij}(t)}{t}\)存在且有限。

证明:略。

Remark:实际应用中,\(P(t)\) 很难获得,一般可以求得 \(Q\),此时 \(e^{tQ} = \lim_{n\to\infty}(I +Qt/n)^n\)。如果取 \(n=2^k\)的形式,只需要 \(k\)次矩阵乘法即可。

推论 7.1:对任意 \(i\inS\)\(0\le \sum_{i\ne j}q_{ij}\leq_{ii}\)

证明: \[q_{ii}=\varliminf_{t\to0} \frac{1-P_{ii}(t)}{t} = \varliminf_{t\to0}\sum_{j\ne i}\frac{P_{ij}(t)}{t} \ge \sum_{j\nei}\varliminf_{t\to0}\frac{P_{ij}(t)}{t} = \sum_{j\ne i}q_{ij}\] 推论 7.2:当 \(S\) 为有限状态空间时,\(\sum_{i\ne j}q_{ij}= q_{ii} <\infty\)。

7.3 Kolmogorov前向后向微分方程

定义:如果 \(Q\)满足 \(\forall i\in S\)\(\sum_{j\ne i}q_{ij} = q_{ii} <\infty\),则称 \(Q\)为保守矩阵。

定理 7.3:设马氏链 \(X=\{X(t),t\ge0\}\),\(Q=P'(0)\),当 \(S\) 为有限集时,\(P'(t)=P(t)Q=QP(t)\)。

Remark:当 \(S\)为可数状态时,前向方程与后向方程不一定成立,根据 Fatu 引理有 \(P'(t)\ge P(t)Q,P'(t)\ge QP(t)\)

定理 7.4:当 \(S\)可列个状态\(Q\)为保守矩阵时,后向方程 \(P'(t)=QP(t)\) 成立。

7.4平稳分布与极限分布及其矩阵计算

类似离散马氏链,可以定义状态的连通关系。

定义:若 \(\int_0^\inftyP_{ii}(t)dt = +\infty\),则称状态 \(i\)为常返态,否则称为非常返态。

定义:设 \(i\)为常返态,若 \(\lim_{t\to\infty}P_{ii}(t) >0\),则称为正常返态,若 \(\lim_{t\to\infty}P_{ii}(t) =0\),称为零常返态。

定义:若概率分布 \(\pi=\{\pi_i,i\in S\}\) 满足 \(\pi=\pi P(t)\),称 \(\pi\) 为 \(X\) 的平稳分布。

定理 7.5:设 \(X\)是有限不可约连续马氏链,对任意 \(i,j\inS\),有 \(\lim_{t\to\infty}P_{ij}(t)=p_j\)存在,且与状态 \(i\) 无关。

证明:略。

🔲 ⭐

【随机过程2】马尔可夫过程2 | 状态空间

6.4 状态空间的分解

定义:设 \(A\subsetS\),若对任意 \(i\in A\)\(j\notin A\),都有 \(p_{ij}=0\),称 \(A\) 为闭集。若 \(A\) 的状态是相通的,则 \(A\) 为不可约的。

引理 6.1\(A\)为闭集的充要条件为:任意 \(i\in A\)\(j\notin A\) 都有 \(p_{ij}^{(n)}=0,n\ge1\)。

推论:若 \(A\)是闭集,对任意状态 \(i\in A\) 恒有\(\sum_{j\in A}p_{ij}^{(n)}=1\)

定理 6.5:所有常返态构成一个闭集。

推论:不可约马尔科夫链,所有状态都是常返的,或者所有状态都是非常返的。

定理 6.6:状态空间 \(S\) 可以分解为 \(S = T\cup C = T\cup C_1 \cdots \cup C_h\cdots\),其中 \(T\)表示非常返状态的集合,\(C_i\)为基本的常返闭集,且有

  1. 对任一确定的 \(k,C_k\)中任意两个状态互通;
  2. \(C_k\cap C_l = \varnothing,\forall h\nel\)

定理 6.7:有限状态马尔科夫链具有如下性质:

  1. 状态空间 \(S\) 可分解为 \(S = T\cup C = T\cup C_1 \cdots \cupC_h\),其中 \(T\)表示非常返状态的集合,\(C_i\)为基本的常返闭集;
  2. 非常返状态集合 \(T\)一定不是闭集;
  3. 没有零常返状态;
  4. 必有正常返状态;
  5. 不可约马氏链的的状态都是正常返态;
  6. 任意闭集 \(C_i\) 上的 \(n\) 步转移矩阵为随机矩阵。

6.5 极限特性与平稳分布

6.5.1 极限特性

定理 6.7:若状态 \(j\) 为非常返态或零常返态,则对任意 \(i \in S\),有 \(\lim_{n\to\infty}p_{ij}^{(n)}=0\)。

推论6.7.1:若马尔科夫链有一个零常返态,则必有无穷多个零常返态。

6.5.2 平稳分布

定义:一个定义在 \(S\) 上的概率分布 \(\pi=(\pi_1,\pi_2,...,\pi_i,...)\)称为马尔科夫链的平稳分布,如果有 \(\pi=\piP\)

定理 6.9:若马尔科夫链是不可约的遍历链,则 \(\{\pi_i = 1 / \mu_i\}\) 是 \(\pi=\pi P(\pi_i\ge0,\sum_{i\in S}\pi_i=1)\)的唯一解

推论 6.9.1:不可约遍历链恒有唯一的平稳分布,且 \(\pi_j=\lim_{n\to\infty}p_{ij}^{(n)}\)。

定理 6.10:令 \(C_+\)为马尔科夫链中全体正常返状态构成的集合,则有

  1. 平稳分布不存在的充要条件为 \(C_+=\varnothing\);
  2. 平稳分布唯一存在的充要条件为只有一个基本正常返闭集;
  3. 若马尔科夫链有多于一个基本正常返闭集,则其平稳分布有无穷多个。

定理 6.11:关于有限状态的马尔科夫链

  1. 有限状态马尔科夫链的平稳分布总存在;
  2. 有限不可约非周期的马尔科夫链存在唯一的平稳分布;
  3. 若有多于一个基本正常返闭集,则其平稳分布有无穷多个;

栗子(平衡方程及其应用):对非周期正常返的离散马氏链,平稳分布存在且满足\(\pi P=\pi\),可以写成 \(\pi_j(1-P_{jj}) = \sum_{i\ne j,i\in S} \pi_iP_{ij}\),这被称为离散马氏链的平衡方程。左边表示从状态 \(j\) 流出的量,右边表示从其他状态流入状态\(j\)的量,在讨论一些实际工程问题时,可以借助平衡方程求解系统平稳分布。

6.6 转移矩阵的平均极限

一般情况下,\(n\to\infty\)\(P^n\) 的极限未必存在,主要是因为 \(P^n\)可能有周期性。为了消除周期性,最直接的方法就是取平均。

定理 6.12:设 \(P\)为有限马氏链的转移矩阵,则 \(L:=\lim_{n\to\infty}\frac{1}{n}(I+P+\cdots+P^{n-1})\) 存在,且满足 \(LP = PL = L = L^2\)。

推论:设有限不可约马氏链的转移矩阵为 \(P\),平稳分布为 \(\pi\),\(L=(l_{ij})=\lim_{n\to\infty}\frac{1}{n}(I+P+\cdots+P^{n-1})\),则 \(\pi L=\pi\),且 \(\pi_j=l_{ij} / \sum_{k=1}^N l_{ik}\)。

定理 6.13:设 \(P\)是有 \(m\)个状态的不可约马氏链的转移矩阵,则平稳分布概率为

\(\pi = (1,...,1)(I-P+{\mathbb1})^{-1}\),其中 \({\mathbb 1}\)为全 1 的 \(m\times m\) 矩阵。

证明:关键是要证明矩阵 \((I-P+{\mathbb1})\) 可逆。

Note:马尔科夫链在随机过程里面是比较简单的部分,大部分结论都很直观,凭直观感觉就能得到。

实际上这章马尔科夫链的主要内容就是在讨论几种状态:非常返态、正常返态、零常返态。归结起来,现在想象一个马尔可夫链的状态转移图,有很多节点(状态)和有向边(转移概率)。

首先对于几种状态的区别:

  1. 正常返态就是说从一个状态开始,经过有限步总能再次回到当前状态,并且这个平均回转时间是有限的(中华好男人,常回家看看);
  2. 零常返态就是说从一个状态开始,经过有限步总能再次回到当前状态,但这个平均回转时间是无穷的(渣男海王,开空头支票);
  3. 非常返态就是说从一个状态开始,经过有限步之后就再也不能返回当前状态了(恩断义绝);

如何判断每个节点的状态类型,基本只需要下面两条原则:

  1. 如果两个状态相通,那么他们同为常返态或非常返态。
    1. 那么如果两个节点之间有双向边,他们一定是同一类型态;要出现常返态和非常返态的区别一定是由于单向边的存在;
    2. 进一步的,一个连通子图里面的所有状态一定是同一类型;要出现常返态和非常返态的区别,一定是两个连通子图\(A,B\) 之间只有单向边;
    3. 假如只有 \(A\)\(B\) 的单向边,那么连通子图 \(B\) 中的节点一定都是非常返态。
  2. 零常返态只可能出现在状态个数为无穷的时候。

极限分布/平稳分布是怎样:

  1. 极限分布与初始分布有关,平稳分布只与状态转移链本身的性质有关;
  2. 由于最终一定会收敛到常返态,平稳分布一定是只对常返态非零;
  3. 只考虑常返态的集合,平稳分布就是转移概率矩阵 \(P\) 的左特征向量;

下面讲个故事。

把概率想象成手里的money,初始分布概率就是每个节点手里的本钱;一个(不可约的)连通子图就是一个家族,家族里的每个人之间可以相互借钱来回周转;不同连通子图就是不同的家族。

(不可约的)连通子图内部,一个家族里面的钱再怎么周转也都是内部消化,每个人借出去的总能还回来,所以钱不会消失,总会有一个平衡。

如果有两个连通子图,并且它们之间只有单向边,相当于 \(A\) 总是把一部分 money 白送给 \(B\),再厚的家底也得被掏空了。最后 \(A\)家族就破产了(非常返态,不是闭集)。

如果 \(B\)家族的钱只进不出,也就是没有向外的有向边(闭集),那么他们就会完成资本的积累,最终钱总会聚集到他们那里(常返态)。

如果有两个连通子图,他们相互之间没有任何边相联系,那就是他们互不打扰,两不相欠(各自有一个平稳分布)。100年之后各自有多少钱取决于现在各自有多少钱,不多不少(联合组成无穷个平稳分布,因为初始分布有无穷种情况)。

🔲 ⭐

【高等数值分析】常微分方程数值解

1. 预备理论

求解常微分方程初值问题数值解 \[\begin{align}&\frac{dy}{dx} = f(x,y), \quad a < x < b, |y| < \infty \\&y(a) = y_0\end{align}\] 存在唯一性定理:若 \(f(x,y)\)连续,对 \(y\) 满足 Lipschitz条件,那么初值问题有唯一解。

对上面的常微分方程积分就有 \(y(t_{n+k}) -y(t_{n-j}) = \int_{t_{n-j}}^{t_{n+k}} f(t,y(t))dt\),所以实际上可以转化为数值积分的问题,所以这一节的方法和数值积分很类似。但是这里的问题在于被积函数值是不知道的。

2. 线性多步法

在数值积分里面 \(f(t_k,y(t_k))\)是已知的,主要是在设计求积节点对应的系数。在这里还需要首先估计每个节点的函数值。对前面的式子采样求积就得到\(y(t_{n+k}) - y(t_{n-j}) = h\sum_i \beta_if_{n-i}\),就引出了下面要介绍的线性多步法。

2.1 线性多步法

\(f_n=f(t_n,y_n),y_n=y(t_n)\),线性多步法(线性 \(k\) 步法)的一般表达式为 \[\sum_{i=0}^k \alpha_i y_{n+i} = h\sum_{i=0}^{k} \beta_i f_{n+i}\] 其中 \(\alpha_k=1,\alpha_0^2+\beta_0^2\ne0\)。给一些特例:

  • \(\alpha_0=-1,\beta_0=1,\beta_1=0\)时为显式Euler法 \(y_{n+1}=y_n +hf_n\)
  • \(\alpha_0=-1,\beta_0=0,\beta_1=1\)时为隐式Euler法 \(y_{n+1}=y_n +hf_{n+1}\)
  • \(\alpha_0=-1,\beta_0=\beta_1=1/2\)时为梯形法 \(y_{n+1} = y_n +h/2(f_n+f_{n+1})\)

可以定义“特征多项式”,在后面分析稳定性和收敛性的时候会很有用\[\rho(\xi) = \sum_{i=0}^k \alpha_i \xi^i, \quad \sigma(\xi)=\sum_{i=0}^k\beta_i\xi^i\]如何衡量数值求解方法的精度呢?定义局部截断误差\[T_{n+k} = \sum_{i=0}^k [\alpha_i y(t_{n+i}) - h\beta_i y'(t_{n+i})]\] 对上面的式子做Taylor展开,可以得到形如 \(T_{n+1} = -1/90 h^5y^{(5)}(x_n)+O(h^6)\)(这是Simpson公式的局部截断误差),这个式子表明Simpson公式是4阶的。

对于线性 \(k\)步法,要想设计一种迭代方法使得数值精度是 \(p\)阶的,可以采用待定系数法,保证局部截断误差中 \(y^{(1)}(x),...,y^{(q)}(x)\)前面的系数都是0。

2.2 稳定性与收敛性

2.2.1 相容性

首先引入相容性的概念。若某个多步法的截断误差满足 \(T(t,y(t),h)=o(h)\)那么称其为相容的。若 \(T(t,y(t),h)=O(h^{q+1})\),则称其为 \(q\) 阶的。

实际上,相容性就是说该方法至少是1阶的,也就是当 \(y\)为一次线性函数时,该方法要能够得到准确解。

定理:相容的多步法充要条件是 \(\rho(1)=0,\rho'(1)=\sigma(1)\)。

2.2.2 零稳定性

根条件:若多项式 \(\rho(\xi)\) 的 \(k\) 个根的模长都不大于1,并且模值等于 1 的根都是单根,则称其满足根条件。

Note:根条件考虑的是齐次差分方程的解。当 \(f\equiv 0\)的时候,齐次差分方程的解具有指数形式 \(y(n)=\sum_i P(n) \xi_i^n\),其中 \(P(n)\) 是一个多项式, \(\xi_i\) 就是 \(\rho(\xi)\)的根(这里没考虑重根的情况,不过是类似的),如果 \(|\xi_i|<1\) 就意味着 \(\lim_{n\to\infty}y(n) =0\),那么这跟稳定性有什么关系呢?

定理:线性多步法关于初值稳定的充要条件是 \(\rho(\xi)\) 满足根条件。

Note:这个定理说明了根条件和(零)稳定性二者的等价关系。为什么呢?假设真实的初值为\(y_0=y(t_0)\),真实的解为 \(y(n)=y_0\xi_0^{(n-t_0)}\)。而我们计算的时候由于各种原因拿到的初值是有误差的,也就是\(\hat{y}_0\),那么最后求得的误差就是\(\hat{y}(n)-y(n)=(\hat{y}_0-y_0)\xi_0^{(n-t_0)}\),如果不满足根条件,当\(n\to\infty\) 的时候,\(\hat{y}(n)-y(n) \to\infty\),解就是不稳定的。

2.2.3 收敛性

定义:假设在区间 \([a,b]\) 上等距划分 \(N\) 个区间,\(x_n= a+nh,n=0,1,...,N\),求解初值问题得到的解为 \(y_n,n=0,1,...,N\),最大整体误差为 \(E(h) = \max_{n} |y(x_n) - y_n|\),如果满足\(\lim_{h\to0}E(h)=0\),则称方法是收敛的。如果 \(E(h) \le Ch^p\),则称方法是 \(p\) 阶收敛的。

定理:相容性 + 零稳定性 \(\iff\) 收敛性。

2.2.4 绝对稳定性

前面的零稳定性是跟所选择的方法有关的。而这里的绝对稳定性研究的是微分方程本身的属性。如果方程本身性质很不好,那么可能无论选择什么方法都是不稳定的。

\(f(x,y)=\lambday\),再研究方法的稳定性。以试验方程为例 \[\begin{cases}y' = \lambda y, t\in[a,b] \\y(a) = \tilde{y}_0\end{cases}\] 用Euler法得到 \(y_n =(1+h\lambda)y_{n-1}\),于是两个解的误差满足 \(y_n-z_n = (1+h\lambda)^n (y_0-z_0)\),若\(|1+h\lambda|>1\),则误差总会放大,这是我们不希望的。

要保证稳定性,既与方程本身的性质(也即 \(\lambda\)) 有关,也与所选择的步长 \(h\)有关。从另一个角度而言,方程本身的性质(\(\lambda\))会影响可选择的 \(h\) 的范围。由此引出绝对稳定性的概念。

对前面提到的线性多步法,把 \(f(x,y)=\lambday\) 代回去就有 \[\sum_{i=0}^k (\alpha_i - h\lambda \beta_i) y_{n+i} = 0\]\(\Pi(\xi;z) = \rho(\xi) -z\sigma(\xi)\)稳定多项式

定理:对给定的 \(z\),若 \(\Pi(\xi;z)\) 的 \(k\) 个根的模都小于1,则其是绝对稳定的。(满足这样条件的 \(z\)的集合构成绝对稳定区域。)

3. Runge-Kutta方法

Runge-Kutta方法和线性多步法的主要区别在于,其在 \(x_n,x_{n+1}\)中间又进行了采样、插值。这种方法大概可以表示为 \[\begin{align}& y_{n+1} = y_n + h\sum_{i=1}^{s} b_i k_i \\& k_i = y_n + f(t_n + c_i h, y_n + \sum_{i=1}^s a_{ij} k_j)\end{align}\] 其中 \(0\le c_i\le1\)。实际上\(k_i \approx y(t_n+c_i h)\),就是在\([y_n,y_{n+1}]\)中间又进行了采样插值。

这样提高了精度,同时也会增加计算复杂度。

除此之外,提高精度的方法还有 Richardson 外推方法。

4. 刚性问题

刚性问题也是方程本身的属性,主要是指某些情况下两个特解的尺度相差很大,比如两个指数衰减的过程混合在一起,但是其中一个衰减特别快(\(\exp(-\lambda t)\) 的 \(\lambda\)特别大),另一个则衰减特别慢,那么数值求解的时候很可能只能看到衰减慢的那个过程,另一个则被忽略。这在化学反应中是经常遇到的。

❌