Kuramoto模型在Python和MATLAB中的简单实现

事情是这样的，我最近在研究团队编组及内部模式对发挥团队能力的影响，以及如何正确编组让团队能力发挥实现最大化，别问我为什么研究这个，反正稀里糊涂就研究上了。我发现在描述团队编组间及内部同步能力的时候，人们对Kuramoto模型（藏本模型）作了大量的研究,其中包括模型达到完全相位同步的充分条件、耦合强度对于同步的影响、一定条件下振子的收敛速率等。但具体实现一般都在MATLAB中，且网上代码过于复杂（我运行了一遍一堆报错），这里我使用Python和MATLAB对Kuramoto模型简单模拟。

模拟的话还是一遍举个栗子，边分析边测试效果最好，百度学术上有一篇关于Kuramoto模型的简单论文，我们就用它来实现模拟。空间信息支援力量编组模式分析，上连接:http://xueshu.baidu.com/usercenter/paper/show?paperid=11b7c2ef69b3ac7f6e114afeb75f8083&site=xueshu_se

好吧，那我们开始，首先是Python实现：

        N维Kuramoto模型的数学描述如下：

（1）

没错，是讨厌的数学公式，没事，它可以改写成这样：

（2）

好像还是有点长，那我们在改写一下：

（3）

看着好多了，那我就来说说式子中参数的意义，Kij为耦合矩阵，是为了便于描述不同振子间耦合程度不同的情形。最下面那个式子的r就是我们的目标，反应振子间的相关性，这个相关性就可以描述我们想要的编组内部同步能力。

哎呦，这个式子看起来好简单，这里要补充一下知识点：同步能力可不是一下子各组该怎么同步直接确定的了，它是一个从开始到稳定的阶段，也就是说随时间变化，最终反映在各组的同步能力才会确定，那么最后图像是什么样子才算同步能力好呢？

同步能力好，是指随着时间的推移，各组的同步能力r逐渐稳定，波动现象消失或固定在某一个小范围内。需要注意的是这和各组r值之间的差距没有关系，我们要的是一个平稳的状态。那怎么办找r和t的关系呢？

注意看最上面那两个式子，相位（第一项，等号左边那个）上面有个点，这样他可就不简简单单是个相位了，它代表的是相位的变化值，是一个微小的微分值，好吧具体意思就是，那个式子左边展开之后是这样的：

（5）

哎呀，t出现了，其实\theta 与t有关，这里你可能有点绕，因为它们之间的关系是一一对应的，就是说每个时间的t对应了一个\theta ，我下面带入具体数值的时候你就知道了。

组间同步能力与时间t的关系出现了！

也就是说我先用上面的那个公式4计算出来\theta 的值，在带入到公式5，那么t-r关系就可以明确下来了，那现在我们再回过头来看看文章中已经给的例子，看看还有没有未知量。

栗子是这样的栗子：

假设某机构内部有 4 个编组，每个编组包括 5 个节点(其中 1 个节点为领导节点) 。另外，将上级领导作为一个独立的编组，且只包含一个节点。假设在领导机关增加4名信息传递人员。当以独立编组模式编组时，指定1名信息传递人员为指挥者，其指挥关系与其他编组一致; 当分散编组时，信息传递力量节点的关系与所在编组其他节点指挥关系一 致。其中，完全分散编组模式时，各信息传递力量节点之间无信息共享通道; 不完全分散编组时，在各信息传递人员节点之间建立一条信息共享通道。各编组模式及其拓扑结构如下图所示：

独立编组模式


完全分散编组


不完全分散编组


参数数据

参数确定一下有没有未知量：

首先N，数据数目已知，这个有了。

K值是分组内的连接强度，这个是看实际情况，由甲方提供或者自己看着给的，这里就是甲方给的编组图，i与j点的链接强度一目了然，这个有了。

\omega _{i} 是振子i的固有频率，也称自然频率，甲方会给，没法自己估计，这个有了。

\theta ，\theta 怎么办，初始的\theta 会给，自己也能测的出来，但那么多\theta _{j} -\theta_{i} 得多少不知道啊，这里通过翻看文章，我发现其实文章是有一个特殊条件的，不然的话是需要研究耦合因子求三种约束条件解情况的，特殊条件就在这：

假设编组内节点的初始相位差为π/2，且编号最小者为0，随编号增大而增大。

哦，初始相位差知道了，你还告诉了我各个初始相位，那么\theta _{j} -\theta _{i} 的值就在一个范围内的几个固定值里面啊！

好的，没有未知量了，就是找K的时候麻烦点，没办法，这个决定了编组的不同，写脚本算一下吧：

#codingutf-8

##ScriptName:KuramotoSimulation.py 

import matplotlib.pyplot as plt 

from pylab import 

from sympy import 

from matplotlib.ticker import MultipleLocator, FormatStrFormatter 

import math 

import numpy as np 

N = 31      #总节点数 

c=[0,0,math.pi  2,math.pi,3  math.pi 2,0,0,0,math.pi  2,math.pi,3  math.pi 2,0,0,0,math.pi  2,math.pi,3  math.pi 2,0,0,0,math.pi  2,math.pi,3  math.pi 2,0,0,0,math.pi  2,math.pi,3  math.pi 2,0,0] 

w = [4,3,3,2,2,1,4,3,3,2,2,1,4,3,3,2,2,1,4,3,3,2,2,1,4,3,3,2,2,1,5] 

k1 = [ 

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2], 

[1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,1,0.9,0.9,0.9,0.9,0.9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2], 

[0,0,0,0,0,0,0.9,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0.9,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0.9,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0.9,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0.9,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,1,0.8,0.8,0.8,0.8,0.8,0,0,0,0,0,0,0,0,0,0,0,0,2], 

[0,0,0,0,0,0,0,0,0,0,0,0,0.8,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0.8,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0.8,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0.8,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0.8,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0.7,0.7,0.7,0.7,0.7,0,0,0,0,0,0,2], 

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.7,1,0,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.7,0,1,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.7,0,0,1,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.7,0,0,0,1,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.7,0,0,0,0,1,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0.6,0.6,0.6,0.6,0.6,2], 

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.6,1,0,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.6,0,1,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.6,0,0,1,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.6,0,0,0,1,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.6,0,0,0,0,1,0], 

[2,0,0,0,0,0,2,0,0,0,0,0,2,0,0,0,0,0,2,0,0,0,0,0,2,0,0,0,0,0,1]]   

n = [i + 1 for i in range(22)]                #目标划分,24个值,1-24 

t = [j  for j in range(1000)] 

ci = 0 

C = [] 

C.append(c) 

for d in range(1100) 

    for i in n 

        for j in range(22) 

            cj = c[j + 1] - c[i] 

            ci += k1[i][j + 1]  math.sin(cj) 

        cii = ci + w[i] 

        h = [0.01  cii + z for z in c] 

        C.append(h) 

        c = h 

def r_function(u) 

    y1 = 0 

    y2 = 0 

    y12 = 0 

    y22 = 0 

    r = [] 

    for x in range(1000) 

        for ul in u 

            y1 += math.cos(C[x][ul]) 

            y2 += math.sin(C[x][ul]) 

        y12 = y1  2 

        y22 = y2  2 

        r.append((float(1)  N )  ((y12 + y22)  0.5)) 

    return r 

r1 = r_function([1,2,3,4,5,6]) 

r2 = r_function([7,8,9,10,11,12]) 

r3 = r_function([13,14,15,16,17,18]) 

r4 = r_function([19,20,21,22,23,24]) 

r5 = r_function([25,26,27,28,29,30]) 

#r6 = r_function([31]) 


ax = subplot(111) #注意一般都在ax中设置,不再plot中设置 

ymajorLocator  = MultipleLocator(0.1) #将y轴主刻度标签设置为0.5的倍数 

ax.yaxis.set_major_locator(ymajorLocator) 

plt.plot(t, r1, marker='o', color='green', label='1') 

plt.plot(t, r2, marker='D',color='red', label='2') 

plt.plot(t, r3, marker='+',color='skyblue', label='3') 

plt.plot(t, r4, marker='h', color='blue', label='4') 

plt.plot(t, r5, marker='_',color='yellow', label='5') 

#plt.plot(t, r6, color='red', label='6') 


plt.legend() # 显示图例 

plt.xlabel('iteration times') 

plt.ylabel('r') 

plt.show() 

独立编组结果如图：

k1独立编组

好的，从图像我们来看看Kuramoto模型在描述这个编组的时候，5组最终稳定，我们说这个团队编组还算科学，但我们改变一下K的值，换成分散编组：

k2 = [

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[1,1,0.9,0.8,0.7,0.6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2], 

[1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0.9,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0.8,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0.7,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0.6,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,1,1,0.9,0.8,0.7,0.6,0,0,0,0,0,0,0,0,0,0,0,2], 

[0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0.9,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0.8,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0.7,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0.6,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,1,0.9,0.8,0.7,0.6,0,0,0,0,0,0,2], 

[0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0.9,0,1,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0.8,0,0,1,0,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0.7,0,0,0,1,0,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0.6,0,0,0,0,1,0,0,0,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0.9,0.8,0.7,0.6,2], 

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.9,0,1,0,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.8,0,0,1,0,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.7,0,0,0,1,0,0], 

[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.6,0,0,0,0,1,0], 

[2,0,0,0,0,0,2,0,0,0,0,0,2,0,0,0,0,0,2,0,0,0,0,0,1] 


k2独立编组

这两幅图像都在开始阶段大幅波动，而后在一定范围内趋于稳定，那么到底哪个分组模式最符合实际，最能突出编组能力呢？

这里还有一个公式，来解决这个问题，编组同步能力的量化：

（6）

M_{s} 就可以描述某编组的同步效果，r_{s} 是达到稳定状态后序参量的均值，β∈(0,1)是调节因子。我们可以用M_{s} 比较编组内部的好坏。那编组间能力的好坏怎样比较呢？

这个Kuramoto模型同样有所考虑，它有一个描述整个系统编组能力的公式：

（7）

其中，P是编组的数量，M_{i} 是第i个编组的同步能力，\omega _{i} 是编组在整个系统中的权重，\psi _{e} 是各编组平均 相位的均值，\Delta \psi 是各编组平均相位的标准差。具体的计算不是这篇文章的重点，就不在计算M_{s} 和M的值来比较上述例子独立编组和分散编组的好坏了，本篇文章主要是讲下Kuramoto模型的解决思路，尤其是上面解决\theta 值的方法可以套用在其他Kuramoto模型中，做一个目标估计绰绰有余的。

下面是解决Kuramoto模型常用的MATLAB编程方法，具体思路与上述基本一致，这里不再赘述，K的值我们给另一种编组模式：不完全分散编组模式，也是现在实际上最长用的编组方式，直接上代码：

clc;

clear all; 

k=[0  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0.5; 

  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   

  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   

  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   

  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   

  0  0  0  0  0  0  1  1  1  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0.5; 

  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   

  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   

  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   

  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   

  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  0  0  0  0  0  0  0  0  0  0.5; 

  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   

  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   

  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   

  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0;   

  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  1  1  1  0  0  0  0  0.5; 

  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0;   

  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0;   

  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0;   

  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  0  0  0  0  0;   

  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  2  2  2  0.5; 

  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  2  0  0  0  0;   

  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  2  0  0  0  0;   

  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  2  0  0  0  0;   

  0.5 0 0  0  0  0.5 0 0  0  0 0.5 0  0  0  0  0.5 0 0  0  0  0.5 0 0  0  0;   

  ]     

N = 25; 

Step= 10000; 

Theta=zeros(Step,N); 

Omega =zeros(Step,N); 

DeltaT=0.01 

Theta(1,:)=[0 pi/2 pi 3*pi/2 0 pi/2 pi 3*pi/2 0 pi/2 pi 3*pi/2 0 pi/2 pi 3*pi/2 0 pi/2 pi 3*pi/2 0 pi/2 pi 3*pi/2 0]; 

Omega(1,:)=2*pi*[2  3  3  4      4 2    3  3      4 4    2  3      3 4    4  2      3 3    4  4      2 2    3  4      1]; 

%  求解微分方程 

for n = 1:Step 

    for i = 1:N 

            Sigma = 0; 

        for j = 1:N 

            %%判断k对同步效果的影响 

= Sigma + k(i,j) * sin( Theta(n,i) - Theta(n,j) ); 

        end 

            Omega(n+1,i) = Omega(n,i) - Sigma; 

            Theta(n+1,i) = Omega(n+1,i)*DeltaT + Theta(n,i); 

    end 

end 

%  求解序参量r(t) 

groupN = 5; 

Step = 1000; 

for i = 1:Step 

    sigma1 = 0; 

    sigma2 = 0; 

    sumTheta = 0; 

    for j =6:10  %表示相应的群组 

        sigma1=sigma1+cos(Theta(i,j)); 

        sigma2=sigma2+sin(Theta(i,j)); 

    end 

    r(i)= sqrt( sigma1^2+sigma2^2 )/groupN; 

end 

x= 0.01:DeltaT:DeltaT*Step; 

plot(x,r); 

MATLAB版不完全分散编组结果如下图：


不完全分散编组（MATLAB）

文章能看到最后真的很不容易，代码及对数据的延伸都已上传至GitHub和Gitee上，求star：

GitHub：https://github.com/wangwei39120157028/Spatial_Information_Support_Force_Grouping_Mode_Analysis

Gitee：https://gitee.com/wwy2018/Spatial_Information_Support_Force_Grouping_Mode_Analysis

欢迎再到我的GitHub和Gitee上看看我的关于逆向工程的项目，点赞求星。

GitHub：https://github.com/wangwei39120157028/IDAPythonScripts

Gitee：https://gitee.com/wwy2018/IDAPythonScripts

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

推荐阅读更多精彩内容