定位算法研究记录-第一阶段

定位算法研究记录-第一阶段

基本介绍

背景

在核设施退役过程中，放射源定位是一个重要的技术手段。

本研究关注放射源定位的数学算法。由于数学算法与放射性测量设备相关，本研究中，重点关注以下几个方向，首先是关注\gamma射线探测器，这意味着放射源的定位通过观测\gamma射线来达成，\gamma射线探测器是观测系统的核心部件之一。其次是无准直的\gamma射线探测器，也就是可以同时测量来自各个方向的\gamma射线，或者称为能在空间上进行4\pi​立体角探测的测量系统，也可以称为无视角场(Field of View, FOV)，对应的是有视角场的探测器，这类探测器的特点是能对前方一定视角的区域成像。再次是未知数量放射源，意思是待定位放射源的数量不能在观测开始之前确定，被观测的环境中可能出现任何数量的放射源，如单一放射源(Orphan Source)或者多个放射源，相比于确定数量放射源的情况，这种未知数量放射源的情况往往需要更加复杂的算法支持.

定位的目标介绍

放射源坐标是定位的最本质的目标，但由于不同活度、半衰期、多放射源、放射源移动、地图材料不确定等条件的引入，每一个条件都会影响探测器的观测，致使定位问题复杂化。因此，为了能够获得放射源的坐标，需要加入对活度、半衰期、放射源数量、放射源移动速度等的考虑。

数学上就是指建立有关上述所有因素的模型，在概率学派和贝叶斯学派中，就是指建立一个概率模型，并且是一个将上述所有因素作为函数的概率模型。由于探测器对光子的观测符合泊松分布，概率模型通常使用泊松分布建立，这个泊松分布所对应的概率通常被用来构建与放射源数量的关系。同时，光子在材料中发生的康普顿散射、光电效应、电子对效应等相互作用，综合起来符合指数衰减规律，宏观上，光子数量与放射源到探测器的距离呈平方反比关系(inverse proportional to the square distance from a srouce to a detector)。而这里的距离就被用来构建与放射源活度、放射源坐标、移动速度、地图材料等因素的关系。

值得一提的是，针对放射源移动的研究较少，在少量研究文献中，Morelande(2009)针对多放射源的移动情况，使用贝叶斯估计构造了概率模型，并配合蒙特卡洛方法进行求解。

定位的数学模型介绍

文献资料中报道的定位算法有三边定位(Trilateration)、极大似然估计方法、贝叶斯估计方法等。

三边定位

三边定位(Trilateration)法是一种几何方法，依靠每一个探测器建立一个测量圆，综合考虑所有测量圆的交叉情况，估计出最合适的放射源坐标。Liu(2010)对这个方法进行过介绍，而且这篇论文是针对定位方法的一片综述性文章，有参考价值。

接下来介绍的方法，其基本原理是将放射源定位的目标，如数量、坐标、强度等，视作参数，配合统计学领域的参数估计方法进行求解。典型的方法包括极大似然估计、贝叶斯估计。参数估计视数量、坐标、强度等信息为三个参数，它们具有等同的重要性。

频率学派-极大似然估计

极大似然估计(Maximum Likelihood Estimation, MLE)是一种统计学思想，Morelande(2007)、Liu(2016)等对放射源定位问题使用了MLE方法进行研究。其中，Morelande(2007)提出使用克拉美-罗界限(Cramer-Rao bound, CRB)来评价估计方法的性能。

极大似然估计(Maximum Likelihood Estimation, MLE)中，最大期望法(Maximum Liklihood Expectation Maximumzation, MLEM)对多放射源同时定位有很好的效果，尤其是以放射源数量未知为前提的定位研究中，MLEM算法十分有效。历史上，MLEM方法首先被应用在计算机断层扫描(CT)中，作为呈像算法使用。

贝叶斯学派-贝叶斯估计

贝叶斯估计考虑了方案(假设)的分布，这个分布是在问题开始之前就获得的，也就是术语"先验分布"(Priori)。由于Lindley's paradox的存在，贝叶斯估计经常被频率学派攻击。

贝叶斯估计是一类方法，它们以贝叶斯估计为基础，通过引入分类算法、PCA等等方法来提升估计效果。贝叶斯估计的特点是提供并不断更新放射源的概率分布地图。Howse(2001)是比较早期的使用MLE方法的文献记录；Tandon(2015, 2016)在贝叶斯估计的基础上引入了大量算法，提出了所谓Bayesian Aggregation(BA，贝叶斯聚合)的方法，但是Tandon的定位算法可能非常粗糙，需要具体阅读其博士论文。

马尔可夫链蒙特卡罗方法(Markov Chain Monte Carlo, MCMC)，马尔可夫链蒙特卡罗方法)被用来解决贝叶斯估计所构造的数学问题，Morelande(2007)、Bukartas(2019)、Hite(2019a)、Hite(2019b)、王明生(2019)等以这个方法配合贝叶斯估计，达成放射源定位的目标。

粒子滤波(Particle Filter)方法，是以贝叶斯估计为基础的另一种求解方法。Rao(2015)使用粒子滤波进行了放射源定位的研究。

定位的数学求解介绍

定位的数学模型构造了数学问题，针对这些问题的求解方法有许多可供选择，可以分为数值方法和统计方法两类。

数值方法可以分为全局优化方法和局部优化方法。

全局优化方法在Stefanescu(2017)的4.1节中列举了，包括模拟退火(Simulating Annealing, SA)算法，粒子群(Particle swarm)算法，遗传算法(Genetic algorithm)。

局部优化方法主要指Implicit Filtering。

数值方法-最小二乘估计

最小二乘估计法，依靠测量数据与放射源坐标、强度的指数衰减关系，构造矩阵。矩阵方程型如AX=B，其中，A是衰减系数矩阵，每一个矩阵元素代表光子从对应位置的放射源发射，直至被对应位置的探测器接收，这一个路程中所发生衰减的比例，主要由光子在材料中的指数衰减过程贡献；B是测量样本矩阵，列向量，为已知数据；X是每一个坐标像素中的放射源强度，为未知数据，是带求解量。X可以通过最小二乘估计进行求解。

更加侧重数学意义的表现，估计不能这么解。

局部优化方法-Implicit Filtering 非线性优化

非线性优化方法是一种最优解寻找的算法，是一种通用方法。我计划将放射源定位问题进行抽象，并构造为一个优化问题，将放射源的位置、强度等抽象为带求解量，找到若干合适的评价标准构造代价函数，并使用非线性优化方法进行求解。非线性优化方法的分支很多，一阶梯度下降法是最基本的方法，二阶梯度下降方法提升了计算精度，牛顿-高斯法、Levenberg-Marquardt法更加适合实际应用。

缺点是可能得到局部最优解。

全局优化方法-蒙特卡罗方法

蒙特卡罗(Monte Carlo, MC)方法，基本靠猜，只要计算时间够，总能碰到最优解。

全局优化方法-模拟退火

数值方法可能被局部最优解干扰，导致估计结果并非全局最优解。模拟退火方法是一个可以找到全局最优解的数值方法，称为全局优化方法，Global Optimization Technique。

Stefanescu(2017)以局部退火方法确定了包含全局最优解的粗略范围，再配合其他方法进行精确求解。具体而言，Stefanescu(2017)使用的是一种适应性的模拟退火方法(An Adaptive SA Algorithm)。

定位问题解的评价

放射源的定位研究可以分为两个部分，第一部分是放射源分布方案的获取，第二部分是方案的评价。前导研究重点介绍方案的评价方法。

统计学上，方案就是似然函数，或者说以似然函数来描述方案，在进行极大似然估计和贝叶斯估计时，都默认似然函数已经获得。由此，将似然函数的数学构造从统计学方法中提取出来，成为定位研究的前导之一。

探测器响应模型

放射源定位的研究目标中，放射源的数量、位置和强度是提及频率最高的，为了不失一般性，本研究同样以这三个变量定位研究目标。

\gamma射线从放射源发出，到被探测器接收并被响应，这个过程可以分为三部分，并统一使用统计学知识进行描述。第一部分是\gamma射线在飞行过程中的衰减过程，第二部分是\gamma进入探测器的有效探测面，第三部分是\gamma被探测器响应。

第一部分与介于放射源和探测器之间的介质有关，使用指数衰减模型进行描述。

设定，探测器总数为N_D，探测器的编号以D为标志，坐标使用三元向量描述，第i个探测器的坐标为\boldsymbol{r}_{Di}。

设定，放射源的总数为N_S，放射源的编号以S为标志，第i个放射源的坐标为\boldsymbol{r}_{Si}、放射性强度为I_{Si}。

设定，\gamma射线在材料中的衰减过程为指数衰减，衰减系数为\mu({\boldsymbol{r})}，总宏观界面(the total macroscopic cross sections of material)为\Sigma(\boldsymbol{r})，\boldsymbol{r}表示某一指定位置，这样设置是考虑到在放射源与探测器之间，存在各种介质，比如空气、混凝土墙壁、树木、土壤等等，不同介质对\gamma射线的衰减效应不同。dis(\boldsymbol{r})表示向量\boldsymbol{r}的距离，也就是\boldsymbol{r}​的二范数。
I_{attenuated,Si,Di}= I_{Si} \cdot exp({- \int_{Si}^{Di} \mu({\boldsymbol{r}}) \cdot \rho({\boldsymbol{r}}}) \cdot \rm{d}(dis(\boldsymbol{r}))) \\ I_{attenuated,Si,Di}= I_{Si} \cdot exp({- \int_{Si}^{Di} \Sigma({\boldsymbol{r}}}) \cdot \rm{d}(dis(\boldsymbol{r})))
第二部分与探测器有效探测面积有关。默认放射源发射粒子的方向是各向同性，在三维空间中，也就是4\pi方向发射。探测器与放射源距离为dis(Di,Si)，探测器的有效探测面积为A_{Di}，设定探测器的有效探测面的法线平行于探测器与放射源的连线\boldsymbol{r_{Si}}-\boldsymbol{r_{Di}}，则所有释放出来的光子中，可以进入探测器有效探测面的数量I_{Si,Di}​表示如下。
I_{attenuatedAndGeo,Si,Di} = I_{attenuated,Si,Di} \cdot \frac{A_{Di}}{4\pi \left\| \boldsymbol{r_{Si}}-\boldsymbol{r_{Di}} \right\|_2^2}
第三部分与探测器的固有探测效率有关。并不是所有进入探测器的光子都能被探测器探测，对于一个能量的光子，探测器能够探测到的光子数量比例是一个固定值，称为固有探测效率，以\epsilon_{Di}表示。被第i个探测器测量到的光子数量I_{Di}​表示如下。
I_{Di,Si,Expected} = I_{attenuatedAndGeo,Si,Di} \cdot \epsilon_{Di}
综合上面三个部分，从第i个放射源释放的\gamma​光子，能够被第i个探测器测量到的数量I_{Di}​表示如下。
I_{Di,Si,Expected} = I_{Si} \cdot exp({- \int_{Si}^{Di} \Sigma({\boldsymbol{r}}}) \cdot \rm{d}(dis(\boldsymbol{r})))\cdot \frac{A_{Di}}{4\pi \left\| \boldsymbol{r_{Si}}-\boldsymbol{r_{Di}} \right\|_2^2} \cdot \epsilon_{Di}

泊松分布和光子的观测概率模型

本研究中，一个探测器对一个放射源的观测过程，是这个探测器对若干光子的响应过程。每一个光子被探测器响应遵循二项分布，对于大量光子入射探测器的情况，探测器对某一数量的光子产生响应的概率过程遵循泊松分布，以Poisson表示。本节以单放射源单探测器观测的情况为例，对概率模型的建立进行介绍。

设定只有一个探测器和一个放射源存在，探测器选用第Di个，放射源选用第Si个。来自第Di个放射源，并且被第Si个探测器响应的光子数量为I_{Di,Si}，那么真实情况下，这个探测器的响应数量I_{Di,Si,Observed}不一定正好是I_{Di,Si}，而是一个以I_{Di,Si}​为期望的泊松分布。

仍然遵守本节内前文的设定，若第Di个探测器对第Si个放射源的响应数量，也就是观测到的\gamma​光子数量为I_{Di,Si,Observed}​，而理论上应该观测到的数量为I_{Di,Si,Expected}​，这种事件的发生概率可以使用泊松分布进行计算。
P(\boldsymbol r_{Si},I_{Si}|Obs_{Di}) \equiv Poisson(I_{Di,Si,Observed},I_{Di,Si,Expected})

P(\boldsymbol r_{Si},I_{Si}|Obs_{Di}) = \frac{{I_{Di,Si,Expected}}^{I_{Di,Si,Observed}}}{I_{Di,Si,Observed}!} \cdot e^{-I_{Di,Si,Expected}}

单放射源多探测器观测的概率模型

若环境中仍然只有一个放射源，但将探测器的数量增加，并且分布于环境中的不同位置。假设经过一次观测，则每一台探测器都有一个属于自己的观测数值，也就是对\gamma​光子的响应次数。设定，放射源编号Si；探测器编号从D0开始，至DN；编号为Di的探测器，所观测到的响应次数为I_{Di,Si,Observed}​，理论上这台探测器的响应次数为I_{Di,Si,Expected}​。对于这种设定下的观测结果，可以将每一个探测器的观测事件发生概率进行累积，并且以乘法的方式累积，用以表示本次观测的总概率P(\boldsymbol r_{Si},I_{Si}|Obs_{All})​，其计算方法如下。
P(\boldsymbol r_{Si},I_{Si}|Obs_{All}) = \prod_{Di=D0}^{DN} {P(\boldsymbol r_{Si},I_{Si}|Obs_{Di})}

P(\boldsymbol r_{Si},I_{Si}|Obs_{All}) = \prod_{Di=D0}^{DN} {\frac{{I_{Di,Si,Expected}}^{I_{Di,Si,Observed}}}{I_{Di,Si,Observed}!} \cdot e^{-I_{Di,Si,Expected}}}

多放射源多探测器观测的概率模型

若环境中有多个放射源，使用多个探测器进行观测，概率模型将变得更加复杂一点。设定，放射源编号从S0开始，至SN；第Si个放射源的坐标为\boldsymbol r_{Si}​、强度为I_{Si}​；探测器编号从D0开始，至DN；编号为Di的探测器，所观测到的响应次数为I_{Di,SAll,Observed}​，理论上这台探测器的响应次数为I_{Di,SAll,Expected}​。
I_{Di,SAll,Expected} = \sum_{Si=S0}^{SN} {I_{Di,Si,Expected}}
通过对照可以发现，这里的探测器响应次数并不是针对某一个放射源的，而是对所有放射源的累积结果，这是因为探测器并不能区分当前观测到的\gamma​光子是来自于哪一个放射源。其实，在一定程度上可以区分，比如采用具有能量分辨能力的探测器，如闪烁体型探测器、半导体型探测器，这属于本研究的一个细分领域，且不与本研究冲突，暂不讨论。

假设在一次观测之后，分析得到了一种放射源分布和探测器观测的方案，编号为Plani。这种方案包含了两个领域的内容，一个是放射源领域，具体包括方案内每一个放射源的坐标和强度；另一个是探测器领域，具体包括一个探测器的响应次数分别对应到每一个放射源上的贡献值，也就是每一个放射源在该探测器上的响应次数。例如，对编号为Di的探测器而言，其观测值是I_{Di,SAll,Expected}​，其中来自编号为Si放射源的响应次数为I_{Di,Si,Observed}​。

那么，编号为Plani的方案发生的概率表示为P(\boldsymbol r_{SAll},I_{SAll}|Obs_{All,Plani})​。
P(\boldsymbol r_{SAll},I_{SAll}|Obs_{All,Plani}) = \prod_{Di=D0}^{DN} \prod_{Si=S0}^{SN} {\frac{{I_{Di,Si,Expected}}^{I_{Di,Si,Observed}}}{I_{Di,Si,Observed}!} \cdot e^{-I_{Di,Si,Expected}}}

似然函数

似然函数是概率统计下贝叶斯估计中的一个名词，结合本文的研究背景，似然函数是以探测器的所有观测结果Obs为条件，第i个备选方案Plani为结果的概率函数，以P(\boldsymbol r_{SAll},I_{SAll}|Obs_{All,Plani})或者{\cal{L}}_{Obs,Plani}({\boldsymbol r}_{SAll},I_{SAll})表示，其中的\boldsymbol r表示所有放射源坐标的集合，I同理，表示所有放射源强度的集合。

所谓"放射源坐标的集合"，实际上就是每一个放射源的坐标的集合，每一个放射源指当前方案Plani内包括的所有放射源。数学上，可以将一个放射源的坐标\boldsymbol r_{Si}视作一个1 \times 3矩阵，将所有坐标纵向排列，形成一个SN \times 3矩阵。

设定，放射源的备选方案包括放射源总数量SN，每一个放射源的位置\boldsymbol r_{Si}和强度I_{Si}​。
\boldsymbol r \equiv {\boldsymbol r}_{All} = \{ {\boldsymbol r_{S0}}, \dots, {\boldsymbol r_{Si}}, \dots, {\boldsymbol r_{SN-1}} \}

I \equiv I_{All} = \{ {I_{S0}}, \dots, {I_{Si}}, \dots, {I_{SN-1}} \}

参考上文可以发现，这里的似然函数就是编号为Plani的方案的发生概率，P(\boldsymbol r_{SAll},I_{SAll}|Obs_{All,Plani})​。
{\cal{L}}_{Obs,Plani}(\boldsymbol r,I) \equiv P(\boldsymbol r_{SAll},I_{SAll}|Obs_{All,Plani})
在一个放射源定位问题的研究中，会出现若干备选方案。这里为每一个方案设定编号，从Plan0开始，至Plan(N-1)。最优解存在于这些方案之中。根据本节开头的描述，每一个方都对应一个似然函数，第i个方案的似然函数表示如下。
{\cal{L}}_{Obs,Plani}(\boldsymbol r,I) \equiv P(\boldsymbol r_{SAll},I_{SAll}|Obs_{All,Plani}) = \prod_{Di=D0}^{DN-1} \prod_{Si=S0}^{SN-1} {\frac{{I_{Di,Si,Expected}}^{I_{Di,Si,Observed}}}{I_{Di,Si,Observed}!} \cdot e^{-I_{Di,Si,Expected}}}

方案的评价

方案，在本研究背景下指放射源的分布方案。每一种方案的优劣，可以理解为一种适配程度，是针对探测器的观测结果的，因此方案的优劣可以表述成当前这种方案对观测结果的适配程度。所谓适配程度，可以理解为某一种分布下的概率值，这个分布由当前方案构造出，形式不唯一，由研究者设计和选择，频率学派的似然函数和贝叶斯学派的后验概率都可以作为这里的分布使用。以频率学派的似然函数为例，对单放射源单探测器情况而言，其分布就是它们所构造出的泊松分布，而对多放射源多探测器情况而言，也就是各个探测器所对应的泊松分布之积。

定位问题细化

前文中所说的方案，就是本研究的目标，众多方案中的最优方案就是定位问题的解。但是，在前导研究中，只说明了在已获得方案的基础上，如何评价每一个方案，并未说明方案的获得方法。这是因为本研究将定位问题划分为方案的获取和方案的评价，前导研究中主要针对第二个部分进行介绍，而第一部分，也就是方案的获取，将在后文介绍。

定位问题的构造和求解

极大似然估计(MLE)

数学问题的构造

极大似然估计通过将似然(Likelihood)函数最大化的方法，寻找最合适的解。在放射源定位中，似然函数简单说就是探测器接收到放射源的概率。但是，极大似然估计只是在构造数学问题，并不提供求解方法。因此，还需要配合数学求解方法，才能进行极大似然估计。求解的方法有很多，MLEM方法是其中研究比较多的，但为了简明表达，先从极大似然的定义本身出发，介绍一种最本质的求解方法。

数学问题的求解-原理性求解

MLE最本质的求解方法，可以从离散和连续两种角度出发，其中离散的角度最简单，因此这里从离散角度考虑MLE问题。由于坐标和活度是一个放射源的基本信息，而环境中放射源数量未知，我们可以不失一般性地设定放射源定位问题的目标是确定放射源的数量、坐标和活度。现在假设已经有若干种备选方案诞生，并且其中一种方案就是真实情况。那么由于备选方案有多种，每一种方案都有一个似然函数，通过比较所有这些似然函数的值，也就是概率，找到最大概率，再找到对应这个最大概率的方案，就认为这个方案是最合适的解。

从贝叶斯的角度(Beyasian perspective)来看，极大似然估计就是一个先验为均匀分布的贝叶斯估计。

数学问题的求解-最大期望法MLEM

最大期望法，简称MLEM方法，以极大似然估计方法构造的数学问题为基础，通过一种迭代的方法进行求解，这种迭代方法被称为最大期望(EM)法。

贝叶斯估计

数学问题的求解-马尔可夫链蒙特卡洛方法(MCMC)

通过使用贝叶斯公式来构造有关放射源的概率方程。王明生(2019)以放射源数量已知为前提，一定程度地简化了数学问题。

放射源定位范例-单放射源极大似然估计

参考文献

[Bukartas, et. al., 2019] A. Bukartas, R. Finck, J. Wallin, C.L. Rääf, A Bayesian method to localize lost gamma sources, Applied Radiation and Isotopes, Volume 145, 2019, Pages 142-147, ISSN 0969-8043, https://doi.org/10.1016/j.apradiso.2018.11.008. (http://www.sciencedirect.com/science/article/pii/S0969804318304019)

[Hite, et al., 2019, a] Jason Hite, John Mattingly, Bayesian Metropolis methods for source localization in an urban environment, Radiation Physics and Chemistry, Volume 155, 2019, Pages 271-274, ISSN 0969-806X, https://doi.org/10.1016/j.radphyschem.2018.06.024. (http://www.sciencedirect.com/science/article/pii/S0969806X17307867)

[Hite, et al., 2019, b] Jason Hite, John Mattingly, Dan Archer, Michael Willis, Andrew Rowe, Kayleigh Bray, Jake Carter, James Ghawaly, Localization of a radioactive source in an urban environment using Bayesian Metropolis methods, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, Volume 915, 2019, Pages 82-93, ISSN 0168-9002, https://doi.org/10.1016/j.nima.2018.09.032. (http://www.sciencedirect.com/science/article/pii/S016890021831177X)

[Howse, et al., 2001] Howse J W , Ticknor L O , Muske K R . Least squares estimation techniques for position tracking of radioactive sources[J]. Automatica, 2001, 37(11):1727-1737.

[Liu, et al., 2010] Liu Y , Yang Z , Wang X , et al. Location, Localization, and Localizability[J]. Journal of Computer Science and Technology, 2010, 25(2):274-297.

[Liu, 2016] Liu Zheng, Mobile radiation sensor networks for source detection in a fluctuating background using geo-tagged count rate data[D], University of Illinois at Urbana-Champaign, 2016. http://hdl.handle.net/2142/90850

[Morelande, et al., 2007] Morelande M R , Kreucher C M , Kastella K . A Bayesian Approach to Multiple Target Detection and Tracking[J]. IEEE Transactions on Signal Processing, 2007, 55(5):1589-1604.

[Morelande, et al., 2009] Morelande, Ristic. Radiological Source Detection and Localisation Using Bayesian Techniques[J]. IEEE Transactions on Signal Processing, 2009, 57(11):4220-4231.

[Rao, et al., 2015] Rao N S V , Sen S , Prins N J , et al. Network algorithms for detection of radiation sources[J]. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 2015, 784:326-331.

[Tandon, 2015] Tandon P . Bayesian Aggregation of Evidence for Detection and Characterization of Patterns in Multiple Noisy Observations[M]. ACM, 2015.

[Tandon, 2016] Tandon P , Huggins P , Maclachlan R , et al. Detection of radioactive sources in urban scenes using Bayesian Aggregation of data from mobile spectrometers[J]. Information Systems, 2016, 57:195-206.

[Stefanescu, et al., 2017] Stefanescu, Razvan, Schmidt K , Hite J , et al. Hybrid optimization and Bayesian inference techniques for a non-smooth radiation detection problem[J]. International Journal for Numerical Methods in Engineering, 2017.

[王明生，2019] 王明生，肖宇峰，刘冉，刘成，杨川，基于自适应 M-H 采样的放射源定位算法，测控技术，2019

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

推荐阅读更多精彩内容