AX=XB 求解

slam 中的手眼标定最后可以归结为求解如下问题
已知 $\boldsymbol A_i, \boldsymbol B_i \in SO(3)$ ，需要求解 $\boldsymbol X \in SO(3)$ ，满足
$\min \sum_{i = 1}^n d(\boldsymbol A_i \boldsymbol X - \boldsymbol X \boldsymbol B_i)$

其中 $SO(3)$ 即为旋转矩阵，网上有多种求解方法，这里只描述其中几种，首先介绍下一些基本的变换

1、基本变换

hat运算符 $\land$

把三维实向量 $\boldsymbol\omega$ 映射到 $3 \times 3$ 反对称阵 $\boldsymbol W$
$\boldsymbol \omega^{\land} =\begin{bmatrix} \omega_1\\ \omega_2\\ \omega_3 \end{bmatrix}^{\land} =\begin{bmatrix} 0 & -\omega_3 & \omega_2 \\ \omega_3 & 0 & -\omega_1 \\ -\omega_2 & \omega_1 & 0 \end{bmatrix} =\boldsymbol W$

性质1： $\boldsymbol a^\land \boldsymbol b = - \boldsymbol b^\land \boldsymbol a$ ，直接推导即可证明

性质2： $\boldsymbol \omega^\land \boldsymbol \omega^\land \boldsymbol \omega^\land = - \| \boldsymbol \omega \|^2 \boldsymbol \omega^\land$ ，直接推导即可证明

性质3： $\boldsymbol \omega \boldsymbol \omega^T - \boldsymbol \omega^\land \boldsymbol \omega^\land = \| \boldsymbol \omega \|^2 \boldsymbol I$ ，直接推导即可证明

性质4： $(\boldsymbol R \boldsymbol v)^\land = \boldsymbol R \boldsymbol v^\land \boldsymbol R^T$ ，证明如下
$\begin{aligned} &\forall \boldsymbol u \in \mathbb R^3, (\boldsymbol R \boldsymbol v) \times (\boldsymbol R \boldsymbol u) = \boldsymbol R(\boldsymbol v \times \boldsymbol u) \\ \Rightarrow &\forall \boldsymbol u \in \mathbb R^3, (\boldsymbol R \boldsymbol v)^\land (\boldsymbol R \boldsymbol u) = \boldsymbol R(\boldsymbol v^\land \boldsymbol u) \\ \Rightarrow &\forall \boldsymbol u \in \mathbb R^3, (\boldsymbol R \boldsymbol v)^\land \boldsymbol R \boldsymbol u = \boldsymbol R \boldsymbol v^\land \boldsymbol u \\ \Rightarrow &(\boldsymbol R \boldsymbol v)^\land \boldsymbol R = \boldsymbol R \boldsymbol v^\land \\ \Rightarrow &(\boldsymbol R \boldsymbol v)^\land = \boldsymbol R \boldsymbol v^\land \boldsymbol R^T \end{aligned}$

vee运算符 $\lor$

把 $3 \times 3$ 反对称阵 $\boldsymbol W$ 映射到三维实向量 $\boldsymbol\omega$
$\boldsymbol W^{\lor} =\begin{bmatrix} 0 & -\omega_3 & \omega_2 \\ \omega_3 & 0 & -\omega_1 \\ -\omega_2 & \omega_1 & 0 \end{bmatrix}^{\lor} =\begin{bmatrix} \omega_1\\ \omega_2\\ \omega_3 \end{bmatrix} =\boldsymbol\omega$

指数映射

把 $3 \times 3$ 反对称阵 $\boldsymbol W$ 映射到旋转矩阵 $\boldsymbol R$
$e^{\boldsymbol W} = \boldsymbol I + \frac{\sin \| \boldsymbol W^\lor \|}{\| \boldsymbol W^\lor \|} \boldsymbol W + \frac{1 - \cos \| \boldsymbol W^\lor \|}{\| \boldsymbol W^\lor \|^2} \boldsymbol W^2 = \boldsymbol R$

性质1：假设 $\theta = || \boldsymbol W^\lor||=\sqrt{\boldsymbol \omega_1^2 + \boldsymbol \omega_2^2 + \boldsymbol \omega_3^2}$ 为向量模长，即 $\boldsymbol \omega = \theta \hat{\boldsymbol \omega}$
$\begin{aligned} e^{\boldsymbol W} &= \sum_{n = 0}^\infty \frac{1}{n!} \boldsymbol W^n = \sum_{n = 0}^\infty \frac{1}{n!} (\boldsymbol \omega^\land)^n = \sum_{n = 0}^\infty \frac{1}{n!} (\theta \hat{\boldsymbol \omega}^\land)^n \\ &= \boldsymbol I + \theta \hat{\boldsymbol \omega}^\land + \frac{1}{2!}\theta^2 \hat{\boldsymbol \omega}^\land \hat{\boldsymbol \omega}^\land + \frac{1}{3!}\theta^3 \hat{\boldsymbol \omega}^\land \hat{\boldsymbol \omega}^\land \hat{\boldsymbol \omega}^\land + \frac{1}{4!}\theta^4 \hat{\boldsymbol \omega}^\land \hat{\boldsymbol \omega}^\land \hat{\boldsymbol \omega}^\land \hat{\boldsymbol \omega}^\land + \cdots \\ &= \boldsymbol I + \theta \hat{\boldsymbol \omega}^\land + \frac{1}{2!}\theta^2 \hat{\boldsymbol \omega}^\land \hat{\boldsymbol \omega}^\land - \frac{1}{3!}\theta^3 \hat{\boldsymbol \omega}^\land - \frac{1}{4!}\theta^4 \hat{\boldsymbol \omega}^\land \hat{\boldsymbol \omega}^\land + \cdots \\ &= \boldsymbol I + \hat{\boldsymbol \omega}^\land \hat{\boldsymbol \omega}^\land + \bigg(\theta - \frac{1}{3!} \theta^3 + \cdots \bigg) \hat{\boldsymbol \omega}^\land - \bigg(1 - \frac{1}{2!} \theta^2 + \frac{1}{4!} \theta^4 - \cdots \bigg) \hat{\boldsymbol \omega}^\land \hat{\boldsymbol \omega}^\land \\ &= \boldsymbol I + \hat{\boldsymbol \omega}^\land \hat{\boldsymbol \omega}^\land +\sin \theta \hat{\boldsymbol \omega}^\land - \cos \theta \hat{\boldsymbol \omega}^\land \hat{\boldsymbol \omega}^\land\\ &= \boldsymbol I + \sin \theta \hat{\boldsymbol \omega}^\land + (1 - \cos \theta) \hat{\boldsymbol \omega}^\land \hat{\boldsymbol \omega}^\land \end{aligned}$

这也正是罗德里格斯公式，可以从旋转几何出发进行推导

性质2： $\forall \boldsymbol X \in SO(3), e^{\boldsymbol X \boldsymbol W \boldsymbol X^T} = \boldsymbol X e^\boldsymbol W \boldsymbol X^T$ ，证明如下
$\begin{aligned} e^{\boldsymbol X \boldsymbol W \boldsymbol X^T} &= \sum_{n=0}^\infty \frac{1}{n!} (\boldsymbol X \boldsymbol W \boldsymbol X^T)^n \\ &= \sum_{n=0}^\infty \frac{1}{n!} \underbrace{\boldsymbol X \boldsymbol W \boldsymbol X^T \boldsymbol X \boldsymbol W\boldsymbol X^T \cdots \boldsymbol X \boldsymbol W\boldsymbol X^T}_n \\ &= \sum_{n=0}^\infty \frac{1}{n!} \boldsymbol X \boldsymbol W^n \boldsymbol X^T \\ &= \boldsymbol X \bigg(\sum_{n=0}^\infty \frac{1}{n!} \boldsymbol W^n \bigg) \boldsymbol X^T \\ &= \boldsymbol X e^{\boldsymbol W} \boldsymbol X^T \end{aligned}$

对数映射

把旋转矩阵 $\boldsymbol R$ 映射到 $3 \times 3$ 反对称阵 $\boldsymbol W$
$\log \boldsymbol R = \frac{\varphi (\boldsymbol R - \boldsymbol R^T)}{2 \sin \varphi} =\boldsymbol W$
其中 $\varphi = \arccos \bigg( \frac{tr(\boldsymbol R) - 1}{2} \bigg)$

2、Navy手眼标定算法

论文 "Robot sensor calibration: solving AX=XB on the Euclidean group"
考虑到 $\boldsymbol A_i, \boldsymbol B_i$ 均为旋转矩阵，因此存在三维实向量 $\boldsymbol \omega_{A_i} = (\log \boldsymbol A_i)^\lor, \boldsymbol \omega_{B_i} = (\log \boldsymbol B_i)^\lor$ 使得
$\boldsymbol A_i = e^{\boldsymbol \omega_{A_i}^\land}, \boldsymbol B_i = e^{\boldsymbol \omega_{B_i}^\land}$

利用前面的性质可以得到
$\begin{aligned} \boldsymbol A_i \boldsymbol X &= \boldsymbol X \boldsymbol B_i \\ \boldsymbol A_i &= \boldsymbol X \boldsymbol B_i \boldsymbol X^T \\ e^{\boldsymbol \omega_{A_i}^\land} &= \boldsymbol X e^{\boldsymbol \omega_{B_i}^\land} \boldsymbol X^T \\ e^{\boldsymbol \omega_{A_i}^\land} &= e^{\boldsymbol X \boldsymbol \omega_{B_i}^\land \boldsymbol X^T} \\ e^{\boldsymbol \omega_{A_i}^\land} &= e^{(\boldsymbol X \boldsymbol \omega_{B_i})^\land} \end{aligned}$

因此求解 $\boldsymbol A_i \boldsymbol X = \boldsymbol X \boldsymbol B_i$ 等价于求解 $\boldsymbol \omega_{A_i} = \boldsymbol X \boldsymbol \omega_{B_i}$ 因此原问题将会变为如下问题
$\min \sum_{i = 1}^n w_i d(\boldsymbol X \boldsymbol \omega_{\boldsymbol B_i} - \boldsymbol \omega_{A_i})$

如果定义距离 $d$ 为欧式距离，则该问题将会变成经典问题
$\min \sum_{i = 1}^n w_i \| \boldsymbol X \boldsymbol \omega_{\boldsymbol B_i} - \boldsymbol \omega_{A_i} \|^2$

具体推导在另一篇文章《点云精配准》求解旋转矩阵 $\boldsymbol R$ 中有详细描述，这里直接给出求解步骤：

找到一系列旋转矩阵 $\boldsymbol A_i, \boldsymbol B_i$
计算对应的三维向量 $\boldsymbol \omega_{A_i} = (\log \boldsymbol A_i)^\lor, \boldsymbol \omega_{B_i} = (\log \boldsymbol B_i)^\lor$
对协方差矩阵 $\boldsymbol \Sigma$ 进行SVD分解 $\boldsymbol \Sigma = \boldsymbol U \boldsymbol \Lambda \boldsymbol V^T$ ，其中协方差矩阵 $\boldsymbol \Sigma$ 为
$\boldsymbol \Sigma = \begin{bmatrix} \boldsymbol \omega_{\boldsymbol B_1} & \cdots & \boldsymbol \omega_{\boldsymbol B_n} \end{bmatrix} \begin{bmatrix} w_1 & 0 & \cdots & 0 \\ 0 & w_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & w_n \end{bmatrix} \begin{bmatrix} \boldsymbol \omega_{\boldsymbol A_1}^T \\ \vdots \\ \boldsymbol \omega_{\boldsymbol A_n}^T \end{bmatrix}$
求解得到旋转矩阵 $\boldsymbol X^* = \boldsymbol V \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & |\boldsymbol V \boldsymbol U^T| \end{bmatrix}\boldsymbol U^T$

带有初值的近似解

如果有一个不错的初值 $\boldsymbol X_0 \in SO(3)$ ，问题变成需要求解得到一个 $\boldsymbol Y$ 使得
$\boldsymbol Y \boldsymbol X_0 \boldsymbol \omega_{\boldsymbol B_i} = \boldsymbol \omega_{A_i}, \quad i = 1,\cdots, n$

注意到 $\boldsymbol Y = \boldsymbol X^* \boldsymbol X_0^T$ ，因此 $\boldsymbol Y \in SO(3)$ 可以用欧拉角表示，同时注意到 $\boldsymbol X_0$ 是一个不错的初值，因此有 $yall, pitch, roll \to 0$ ，进而 $\boldsymbol Y$ 可以近似表示成
$\begin{aligned} \boldsymbol Y =& \boldsymbol R_z(yall) \boldsymbol R_y(pitch) \boldsymbol R_x(roll) \\ =& \begin{bmatrix} \cos y & -\sin y & 0\\ \sin y & \cos y & 0\\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} \cos p & 0 & \sin p\\ 0 & 1 & 0\\ -\sin p & 0 & \cos p \end{bmatrix} \begin{bmatrix} 1 & 0 & 0\\ 0 & \cos r & -\sin r\\ 0 & \sin r & \cos r\\ \end{bmatrix} \\ \approx& \begin{bmatrix} 1 & -y & 0\\ y & 1 & 0\\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & p\\ 0 & 1 & 0\\ -p & 0 &1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & -r\\ 0 & r & 1\\ \end{bmatrix} \\ \approx& \begin{bmatrix} 1 & -y & p\\ y & 1 & -r\\ -p & r & 1\\ \end{bmatrix} \\ =& \boldsymbol I + \begin{bmatrix} r\\ p\\ y \end{bmatrix}^\land \end{aligned}$

带入前式可以得到
$\boldsymbol \omega_{A_i} = \begin{bmatrix} r\\ p\\ y \end{bmatrix}^\land \boldsymbol X_0 \boldsymbol \omega_{\boldsymbol B_i} + \boldsymbol X_0 \boldsymbol \omega_{\boldsymbol B_i} = -(\boldsymbol X_0 \boldsymbol \omega_{\boldsymbol B_i})^\land \begin{bmatrix} r\\ p\\ y \end{bmatrix} + \boldsymbol X_0 \boldsymbol \omega_{\boldsymbol B_i}$

因此原问题变成
$(\boldsymbol X_0 \boldsymbol \omega_{\boldsymbol B_i})^\land \begin{bmatrix} r\\ p\\ y \end{bmatrix} = \boldsymbol X_0 \boldsymbol \omega_{\boldsymbol B_i} - \boldsymbol \omega_{A_i}, \quad i = 1,\cdots,n$

将 $n$ 个方程可以写成矩阵形式
$\begin{bmatrix} (\boldsymbol X_0 \boldsymbol \omega_{\boldsymbol B_1})^\land \\ \vdots \\ (\boldsymbol X_0 \boldsymbol \omega_{\boldsymbol B_n})^\land \end{bmatrix} \begin{bmatrix} r\\ p\\ y \end{bmatrix} = \begin{bmatrix} \boldsymbol X_0 \boldsymbol \omega_{\boldsymbol B_1} - \boldsymbol \omega_{A_1} \\ \vdots \\ \boldsymbol X_0 \boldsymbol \omega_{\boldsymbol B_n} - \boldsymbol \omega_{A_n} \end{bmatrix}$

此问题为经典 $\boldsymbol A \boldsymbol x = \boldsymbol b$ 问题，其解为
$\boldsymbol x = (\boldsymbol A^T \boldsymbol A)^{-1} \boldsymbol A^T \boldsymbol b = \boldsymbol V \Lambda^+ \boldsymbol U^T \boldsymbol b$

其中 $\boldsymbol A = \boldsymbol U \Lambda \boldsymbol V^T$ ， $\Lambda^+$ 为伪逆
$\Lambda^+ = \begin{bmatrix} \frac{1}{\lambda_1}\\ & \ddots\\ && \frac{1}{\lambda_m}\\ &&& 0\\ &&&& \ddots\\ &&&&& 0 \end{bmatrix}$

误差分析

对于优化方程 $f(x, z_1,\cdots,z_n)$ ，其中 $x$ 为需要求解的变量， $z_i$ 为观测得到的变量，如果通过求导来求得最优解 $x^*$ ，那么有 $\frac{\partial f(x^*, z_1,\cdots,z_n)}{\partial x} = F(x^*, z_1,\cdots,z_n) = 0$ ，该方程通过计算可以得到 $x^* = g(z_1,\cdots,z_n)$ ，假设 $z_i$ 均在真值 $z_i^0$ 附近，将函数 $g$ 在真值处进行泰勒展开可以得到
$x^*=g(z_1,\cdots,z_n) \approx g(z_1^0,\cdots,z_n^0)+\sum_{i=1}^n\frac{\partial g}{\partial z_i}(z_i - z_i^0)$

如果观测量 $z_i$ 相互之间独立，则有
$cov(x) \approx \sum_{i = 1}^n \frac{\partial g}{\partial z_i} cov(z_i) \frac{\partial g}{\partial z_i}^T$

根据隐函数定理
$\begin{aligned} \frac{\partial g}{\partial z_i} =& -\bigg( \frac{\partial F}{\partial x} \bigg)^{-1}\frac{\partial F}{\partial z_i} \\ =& -\bigg( \frac{\partial^2 f}{\partial x^2} \bigg)^{-1} \frac{\partial^2 f}{\partial z_i \partial x} \end{aligned}$

由此可以得到最优解 $x^*$ 的协方差矩阵为
$cov(x) \approx \sum_{i = 1}^n \bigg( \frac{\partial^2 f}{\partial x^2} \bigg)^{-1} \frac{\partial^2 f}{\partial z_i \partial x} cov(z_i) \frac{\partial^2 f}{\partial z_i \partial x}^T \bigg( \frac{\partial^2 f}{\partial x^2} \bigg)^{-1}$

当前
$\begin{aligned} f(\boldsymbol{\omega}_{\boldsymbol{X}}, \boldsymbol{\omega}_{\boldsymbol{B}_i}, \boldsymbol{\omega}_{\boldsymbol{A}_i}) =& \sum_{i = 1}^n w_i \| e^{\boldsymbol{\omega}_{\boldsymbol{X}}^\land} \boldsymbol{\omega}_{\boldsymbol{B}_i}- \boldsymbol{\omega}_{\boldsymbol{A}_i} \|^2\\ =& \sum_{i = 1}^n w_i (\boldsymbol{\omega}_{\boldsymbol{B}_i}^T {e^{\boldsymbol{\omega}_{\boldsymbol{X}}^\land}}^T e^{\boldsymbol{\omega}_{\boldsymbol{X}}^\land} \boldsymbol{\omega}_{\boldsymbol{B}_i} - 2\boldsymbol{\omega}_{\boldsymbol{A}_i}^T e^{\boldsymbol{\omega}_{\boldsymbol{X}}^\land} \boldsymbol{\omega}_{\boldsymbol{B}_i} + \boldsymbol{\omega}_{\boldsymbol{A}_i}^T \boldsymbol{\omega}_{\boldsymbol{A}_i} )\\ =& \sum_{i = 1}^n w_i (\boldsymbol{\omega}_{\boldsymbol{B}_i}^T\boldsymbol{\omega}_{\boldsymbol{B}_i} - 2\boldsymbol{\omega}_{\boldsymbol{A}_i}^T e^{\boldsymbol{\omega}_{\boldsymbol{X}}^\land} \boldsymbol{\omega}_{\boldsymbol{B}_i} + \boldsymbol{\omega}_{\boldsymbol{A}_i}^T \boldsymbol{\omega}_{\boldsymbol{A}_i} ) \end{aligned}$

采用右扰动，推导得到一阶导
$\frac{\partial f}{\partial \boldsymbol{\omega_X}} = 2\sum_{i=1}^n w_i (\boldsymbol{X}^T \boldsymbol{\omega_{A_i}})^\land \boldsymbol{\omega_{B_i}} = -2\sum_{i=1}^n w_i \boldsymbol{\omega_{B_i}}^\land \boldsymbol{X}^T \boldsymbol{\omega_{A_i}}$

继续推导得到二阶导
$\begin{aligned} \frac{\partial^2 f}{\partial \boldsymbol{\omega_X}^2} =& -2\sum_{i=1}^n w_i (\boldsymbol{X}^T \boldsymbol{\omega_{A_i}})^\land \boldsymbol{\omega_{B_i}}^\land = -2\sum_{i=1}^n w_i (\boldsymbol{\omega_{B_i}}^\land)^2 \\ \frac{\partial^2 f}{\partial \boldsymbol{\omega_{B_i}}\partial \boldsymbol{\omega_X}} =& -2w_i (\boldsymbol{X}^T \boldsymbol{\omega_{A_i}})^\land = -2w_i \boldsymbol{\omega_{B_i}}^\land \\ \frac{\partial^2 f}{\partial \boldsymbol{\omega_{A_i}}\partial \boldsymbol{\omega_X}} =& 2w_i \boldsymbol{X} \boldsymbol{\omega_{B_i}}^\land \end{aligned}$

采用左扰动，推导得到一阶导
$\frac{\partial f}{\partial \boldsymbol{\omega_X}} = 2\sum_{i=1}^n w_i \boldsymbol\omega_{A_i}^\land \boldsymbol{X} \boldsymbol{\omega_{B_i}} = -2\sum_{i=1}^n w_i (\boldsymbol{X} \boldsymbol{\omega_{B_i}})^\land \boldsymbol{\omega_{A_i}}$

继续推导得到二阶导
$\begin{aligned} \frac{\partial^2 f}{\partial \boldsymbol{\omega_X}^2} =& -2\sum_{i=1}^n w_i (\boldsymbol{X}\boldsymbol{\omega_{B_i}})^\land \boldsymbol{\omega_{A_i}}^\land = -2\sum_{i=1}^n w_i (\boldsymbol{\omega_{A_i}}^\land)^2 \\ \frac{\partial^2 f}{\partial \boldsymbol{\omega_{B_i}}\partial \boldsymbol{\omega_X}} =& -2w_i \boldsymbol{X}^T \boldsymbol{\omega_{A_i}}^\land\\ \frac{\partial^2 f}{\partial \boldsymbol{\omega_{A_i}}\partial \boldsymbol{\omega_X}} =& 2w_i (\boldsymbol{X} \boldsymbol{\omega_{B_i}})^\land = 2w_i \boldsymbol{\omega_{A_i}}^\land \end{aligned}$

因此有
$\begin{aligned} \Delta\boldsymbol{\omega_X} =& -\sum_{i=1}^n \bigg( \frac{\partial^2 f}{\partial \boldsymbol{\omega_X}^2} \bigg)^{-1} \frac{\partial^2 f}{\partial \boldsymbol{\omega_{B_i}}\partial \boldsymbol{\omega_X}} \Delta\boldsymbol{\omega_{B_i}} \\ &- \sum_{i=1}^n \bigg( \frac{\partial^2 f}{\partial \boldsymbol{\omega_X}^2} \bigg)^{-1} \frac{\partial^2 f}{\partial \boldsymbol{\omega_{A_i}}\partial \boldsymbol{\omega_X}} \Delta\boldsymbol{\omega_{A_i}} \\ =& -\bigg( \sum_{i=1}^n w_i(\boldsymbol{\omega_{B_i}}^\land)^2 \bigg)^{-1} \sum_{i=1}^n w_i \boldsymbol{\omega_{B_i}}^\land \Delta\boldsymbol{\omega_{B_i}}\\ &+\bigg( \sum_{i=1}^n w_i(\boldsymbol{\omega_{A_i}}^\land)^2 \bigg)^{-1} \sum_{i=1}^n w_i \boldsymbol{\omega_{A_i}}^\land \Delta\boldsymbol{\omega_{A_i}} \end{aligned}$

最后编辑于：2023.05.12 16:41:41

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 205,132评论 6赞 478
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 87,802评论 2赞 381
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 151,566评论 0赞 338
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 54,858评论 1赞 277
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 63,867评论 5赞 368
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 48,695评论 1赞 282
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 38,064评论 3赞 399
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 36,705评论 0赞 258
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 42,915评论 1赞 300
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,677评论 2赞 323
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,796评论 1赞 333
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 33,432评论 4赞 322
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 39,041评论 3赞 307
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,992评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 31,223评论 1赞 260
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 45,185评论 2赞 352
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 42,535评论 2赞 343