EEG巨型分析I:跨研究的频谱和振幅特征

EEG巨型分析I:跨研究的频谱和振幅特征

原创 周翊 茗创科技 

导读

通过汇集多项研究的统计结果(元分析),fMRI领域取得了重大成就。最近,fMRI标准化工作的重点是实现跨研究(巨型分析)的fMRI原始数据的联合分析,以期获得更详细的见解。然而,目前尚不清楚在EEG领域的此类分析是否可能或同样富有成效。在这里,研究者使用了来自六个地点的18项研究,并提供了一个大规模EEG巨型分析结果,这些研究代表了几种不同的实验范式。研究结果证明,当跨研究的元数据一致时,通道级和源级EEG巨型分析都是可能的,并且可以提供单一研究中无法提供的见解。该分析使用全自动处理流程来减少线噪声,插值噪声通道,执行鲁棒参考,去除眼电活动,并进一步识别异常信号。研究者定义了几个基于通道振幅和分布的稳健测量,以评估跨研究数据的可比性,并观察各种处理步骤对这些测量的影响。使用基于ICA的偶极源,还观察到整个脑区的总体频率基线振幅的一致差异。例如,观察到后脑区的alpha比前脑区高,颞区的beta值更高。总的来说,本研究应用了巨型分析来评估跨研究的事件相关EEG特征的共性。本分析中使用的连续原始数据和预处理数据可在https://cancta.net上的DataCatalog获得。

前言

事实证明,fMRI研究的数据汇集对于提高统计功效、评估主体间差异以及确定预测的普遍性和再现性非常有价值。fMRI元分析通常结合了基于峰值激活坐标(基于坐标的元分析,CBMA)或激活的3D统计图像(基于图像的元分析,IBMA)的研究。BrainMap项目是建立空间标准化方法、引用公共坐标,并将坐标和图像联结到相关文献的先驱。这些和其他fMRI数据共享工作使得fMRI元分析呈爆炸式增长。这增加了对大脑基本功能的理解,并促进了用于诊断和评估精神疾病治疗的工具和生物标志物的开发。

Poldrack等人(2017)为透明和可重复的fMRI成像制定了全面的策略,并为标准化的自动化处理、基准和大规模共享原始数据以及统计图的需求提供了强有力的论证。BIDS(脑成像数据结构)规范代表了为原始和处理后的神经影像数据以及相关元数据制定数据标准的重要国际努力。fMRI数据的BIDS标准化和支持工具相对成熟。与元分析相比,巨型分析涉及跨记录或研究的原始数据的汇集和联合分析,而不是结合元数据或派生的统计数据。随着用于共享原始数据的开放访问平台的创建,例如OpenfMRI及其后继者OpenNeuro,fMRI巨型分析现在开始出现了。

虽然OpenNeuro目前主要是fMRI数据,但最近对MEG和EEG成像的BIDS扩展大大增加了在OpenNeuro和其他平台上公开EEG数据的兴趣。不幸的是,目前成功的fMRI元分析的标准化统计图没有为EEG/MEG提供等效支持。此外,尽管最近出现了一些关于EEG/MEG数据采集和预处理最佳实践的通用指南,但MEG/EEG预处理尚未标准化,只是刚开始进行系统地基准测试。

巨型分析将记录和研究中的原始数据结合起来,其核心问题是EEG数据应该如何标准化,以提高跨研究信号的可比性。本文通过在六个不同实验地点进行的18项研究中,考察了通道级和源级的空间、时间和频谱信号特性,旨在为这些问题提供初步指导。一篇相关论文应用了更详细的巨型分析和层次统计建模,以评估跨同一语料库中与事件相关的时间和频谱特征的共性。

该语料库包括几组研究,它们使用相同的一般范式,但在协议细节上有所不同。例如,其中五项研究是基于RSVP(快速串行视觉呈现)范式,但具有不同类型的图像、呈现速率以及对标记目标检测的手动响应要求。其他研究则在驾驶模拟中使用诸如车辆扰动下的车道保持等任务。该语料库反映了参与研究小组的兴趣,表明了在公共存储库中可能出现的研究类型和异质性,其中研究间的控制变量可能无法比较,即使在内部使用多因素设计时也是如此。在综合分析中很难系统地计算跨被试的变异性。然而,即使在公共存储库可能遇到的非结构化环境中,也有可以探索和利用的通用元素。

在过去的十年中,人们越来越强调在EEG数据分析中使用多层次统计模型(如分层线性模型或混合效应模型)来联合分析被试、组水平和研究的内部和之间的变异性。在部分汇集分析中明确建模记录和研究特定的可变性,可以帮助解决与经典全汇集分析相关的I型错误率膨胀问题。在EEG分析中,这种多层次模型通常应用于事件相关数据的结构化语料库。在本文中,研究者使用全汇集分析研究了数据的一般特征。参照相关论文进行了更详细的多层次统计分析,其中描述了各种实验因素对事件相关时间和频谱特征的影响,同时使用两级层次线性模型显式地建模记录特定的可变性。

本文分为两个主要主题,即考察跨研究的通道级和源级的空间、时间和频谱特征共性。方法部分简要描述了实验数据,概述了用于执行分析的自动化处理流程,并介绍了一类稳健的测量,以探索跨异构研究的EEG的标准化和统计特性。结果部分展示了在通道和源空间中EEG信号的统计差异和相似的结果。总之,本研究提供了对包含在大型、多样化EEG记录语料库中的变异性来源的广泛见解,以及如何解决这些因素以揭示跨频段和脑区的功率谱振幅的一致差异。

方法

实验数据

如表1所示,本文使用了在六个地点进行的18项研究的数据,这些地点与四个机构有关:陆军研究实验室(ARL)、国立交通大学(NCTU)、荷兰应用科学研究所(TNO)和加州大学圣地亚哥分校(UCSD)。所有研究均在被试自愿、充分知情的情况下进行,并获得各自院校机构的评审委员会批准。

表1.数据集概要。


这些研究包括两大类任务的变化:视觉目标检测和有/无运动平台的驾驶、交通、速度控制、分心、听觉反馈或自适应巡航控制(ACC)。视觉目标检测任务包括模拟警卫任务的ID检查(GUARDA和GUARDB)、快速串行视觉呈现(RSVPB、RSVPC、RSVPE、RSVPI、RSVPU)和视觉奇异(VEP)。除ACC外,其他任务都是基于驾驶模拟器中的车道保持变化,ACC是在测试轨道上驾驶实际车辆。该库共包含1173条EEG记录,持续时间为633小时。NCTU、TNO和UCSD实验均在各自机构的同一地点进行,而ARL实验则在三个不同地点进行。ARLH指定在陆军研究实验室的人类研究与工程理事会进行实验,ARLS指定在SAIC实验室进行实验,ARLT指定在Teledyne公司进行实验。ARLA指定在所有三个地点使用相同的实验设置执行任务(LKBase和LKCal)。为了便于下游处理的一致性,整个语料库采用ESS(EEG研究模式)组装。将单个EEG记录转换为EEGLAB格式,并将特定于研究的事件代码映射到分层事件描述符(HED)字符串。

数据组织和预处理

数据处理是完全自动化的,并且分几个阶段进行(见图1)。使用PREP,这是本研究作者开发的一种自动化开源工具,用于执行稳健的平均参考和插值坏导。从每个数据记录中去除线噪声后,PREP使用迭代过程来识别和插值噪声通道,以计算稳健的平均参考。在最终迭代中识别为坏的通道被插值。从每个记录中删除非EEG通道,并根据到标准10-20位置的通道距离为没有标准标签的通道分配10-20通道标签。在基于64导10-20配置选择最多64个通道进行此分析后,使用EEGLAB pop_resample()以128 Hz重采样数据。由于本文中使用的所有18项研究都是以高于128 Hz的采样率获得的,因此重采样实际上是进行降采样。使用EEGLAB的pop_eegfiltnew()函数在1Hz下应用1690点高通滤波来去除低频漂移。


图1.EEG自动大规模分析管道。

使用独立成分分析(ICA)和眼动残差回归来从语料库中去除眼动。在使用Mullen等人(2015)所述的无异常值振幅的数据部分执行Infomax ICA后,应用了修改版的EyeCatch来识别与眼动相关的独立成分(IC)。EyeCatch包含3453个手动识别的眼动ICA头皮地形图数据库,将眼睛IC标记为与数据库中至少一个头皮地形图具有相似的地形(相关性>0.94),或者是适度相似(相关性>0.85)且至少1%的数据帧具有至少100的功率比的IC。为了计算功率比,EyeCatch使用MATLAB cwt()函数执行连续小波变换,并在每个时间点上计算[1,3] Hz范围内的功率除以(3,15] Hz范围内的功率。然后,从通道数据中删除由EyeCatch识别为包含眼动的ICs所跨越的子空间。此外,应用BLINKER以识别在去除眼动伪迹之前与眨眼不同阶段相关的潜伏期。BLINKER在EEG.event structure中插入了五个事件(Left base,Left zero,Blink peak,Right zero,Right base)。BLINKER还提取了连续的“眨眼”信号,该信号遵循每对标记之间眨眼引起的EEG活动。

通道分析

通道分析的目标是评估跨研究通道数据的统计和可变性。一些分析步骤通过去除中值并除以每个通道的稳健标准差来对数据进行标准化。将这种标准化称为“稳健的z分数”数据转换。稳健标准差定义为1.4826倍的绝对中位差(MAD)。本文将EEG通道的稳健标准差称为“稳健通道振幅”。

头皮位置和记录通道振幅的变化:为了研究EEG信号在通道位置和记录中是否始终具有不同的振幅,通过26×1正向量开发了每个记录振幅的简化表示,称为“记录振幅矢量”。首先选择所有记录中共有的26个通道来计算记录振幅矢量:Fp1、Fp2、F3、Fz、F4、F7、F8、FC3、FCz、FC4、FT7、FT8、C3、Cz、C4、TP7、TP8、CP3、CPz、CP4、P3、Pz、P4、O1、Oz和O2。然后,计算每个通道在20Hz低通滤波通道信号后的稳健标准差。为了在多个R记录中进行分析,这里形成了“振幅矩阵”,即这些列向量的26×R矩阵。(注意,这种低通滤波和通道下选仅用于通道分析目的,不适用于源计算。)

通过记录特定的常数来表征变化:EEG记录联合分析的一个已知困难是,耳机技术和传感器放置以及被试头发、头皮和颅骨特征引入的可变性会导致数据的统计特性不同。解决此问题的一种简单方法是用特定于记录的比例因子对可变性进行建模,并通过将其通道数据除以特定于记录的常数来“标准化”每个记录。本研究考察了几种计算特定记录的标准化因子,包括均值、Huber均值和记录振幅矢量的欧氏(L2)范数。还研究了直接基于通道信号而不是记录振幅矢量的测量值,但这些方法的可比性较差,因此这里没有报告。虽然所有这些方法都捕获了记录振幅的集中趋势,但均值和欧氏范数对异常值很敏感。Huber均值是一种迭代技术,用于在异常值存在的情况下稳健地逼近均值。

标准化因子是在振幅矩阵(跨通道)上按列计算的,如图2中概括的Huber均值。无论采用何种标准化方法,得到的振幅矩阵都按列除以记录标准化因子的1×R向量,以生成每种标准化方法的归一化振幅矩阵。


图2.计算Huber均值的步骤。其他归一化方法也采用了类似的步骤。

跨头皮位置和频率的通道基线幅值变化:为了研究不同频率的EEG信号是否在通道位置上始终具有不同的振幅,研究者分析了去除眼动,并通过Huber记录特定常数进行归一化的连续EEG。选择用于研究通道信号振幅的26个相同通道子集,并计算每个通道的连续小波变换(使用MATLAB cwt()函数),以获得时变振幅谱图。为了关注每个通道中不同频率的相对幅值,将这些基线频谱幅值矢量归一化,使欧式范数为1(1)。这里将结果称为“归一化通道频谱图”。

为了研究在每个频率上观察到的通道幅值差异的统计显著性,对每个(频率、记录)对的通道频谱幅值的稳健z分数应用Wilcoxon符号秩检验(MATLAB signrank()函数)。这些统计检验的零假设是,对于每个(频率、记录)对,每个通道频谱振幅可以由一个常数加上一个从零均值化分布中抽取的随机值来建模。

源分析

使用来自ICA的源对脑源信号变化进行评估,这是一种广泛使用的EEG源分析技术。使用在预处理过程中计算的ICA成分(非眼动)作为分析的起点。

脑源在空间和频率上的变化:通过将EEGLAB dipfit插件应用于在预处理中未被识别为眼动伪迹的ICs头皮图,推断出最大独立EEG活动源的等效偶极子位置。为了进一步选择最可能与脑源相关的ICs,本研究只考虑了符合以下所有标准的偶极子:(a)由偶极子解释的头皮图至少占85%的方差,(b)EEGLAB MARA插件没有将成分头皮图识别为伪迹,并且(c)偶极子位于脑体积内,由EEGLAB fieldtrip_lite插件的ft_sourcedepth()函数识别。注意,每个偶极子只与一个记录相关联。与通道空间分析一样,研究者计算了与每个偶极子相关的IC激活的连续小波变换(使用 MATLAB cwt()函数),以获得时变振幅谱图。为了关注不同频率的相对幅值,并方便交叉记录比较,对每个频谱图进行缩放,使其总体基线频谱振幅矢量的欧氏范数为1。这里将结果称为IC的“归一化频谱图”。

不同脑区的背景频谱功率特征:通道级或源级电生理信号的功率谱可以建模为周期性和非周期性成分的组合,分别反映窄带振荡(频峰)和约1/f背景波动。虽然已经广泛研究了各种频段中周期性/振荡活动的特性,并将其与一系列认知和行为状态以及神经系统疾病联系起来,但功率谱的非周期性成分却很少受到关注。为了表征非周期谱的空间分布,研究者计算了每个IC的所有时间点的Huber平均功率谱。使用Huber均值而不是中值,以便更准确地估计平均功率,同时忽略由伪迹引起的异常值。然后,使用fooof python工具箱进行建模。

结果

通道分析

为了研究将每个EEG记录除以单个记录特定的数字是否提高了记录之间的可比性,研究者考察了不同归一化方法和不同处理阶段的通道振幅离散。图3显示了不同类型的处理对通道振幅离散的影响。


图3.预处理和通道归一化对通道振幅离散的影响。每个头皮图显示了每个通道的通道离散度中值。括号内的数字对应于头皮图的总体中值离散度。所有归一化方法的结果是在去除眼动伪迹后得到的。

有趣的是,虽然去除眼动伪迹改变了离散的空间分布,但它并没有显著减少离差量。此外,在去除眼动伪迹之前的最大离差并不会出现在额通道中。注意,离散是一个无量纲量,按常数缩放记录不会改变离散度。与非归一化记录相比,使用三种方法(均值、L2或Huber)中的任何一种进行归一化都能显著降低通道离散度。对于没有异常值的数据集,Huber均值和均值本质上是相同的,但是当存在异常值时,Huber均值法产生的离散度略低。对于所有三种归一化方法的离散空间分布是相似的。

图4显示了去除动伪迹前(黑)、去除眼动伪迹后(蓝)以及去除眼动伪迹后进行Huber归一化(绿)的通道离散分布情况。Huber归一化显著降低了所有实验任务的中值离散度。对于大多数实验,归一化还减少了通道离散异常值的数量。注意,LKBase和LKCal是基本的车道保持研究,包括在三个不同的实验地点使用两种不同的耳机类型获得的记录。


图4.在去除动伪迹前(黑)、去除眼动伪迹后(蓝)以及去除眼动伪迹后进行Huber归一化(绿)的通道离散分布。

眼动伪迹的去除和记录特定常数的缩放都会影响通道幅值的相对分布,如图5所示。图5的第一行显示了在复合通道振幅矩阵上的计算结果,对应于128Hz降采样的连续EEG数据(在去除眼动伪迹前)。中间行显示了去除眼动伪迹后在复合通道振幅矩阵上的计算结果。底行显示了在去除眼动伪迹,并按记录特定值(Huber均值)缩放每个记录后,基于复合通道振幅矩阵的计算结果。


图5

图5的第一列显示了各个通道振幅矩阵的头皮地形图:分别在去除眼动伪迹之前、去除眼动伪迹之后和Huber归一化之后。使用EEGLAB topoplot()函数和方法部分中列出的26个常用通道。去除眼动伪迹后,稳健通道振幅中值降低(3.76 vs 4.97),稳健通道振幅中值的空间分布由集中在额叶区域变为双叶形式,幅值分布更均匀,主要分布在额叶和枕叶区域。注意,归一化不会显著改变地形图的空间分布,因为归一化只对应于对中值有贡献的点的重新加权。然而,归一化确实大大降低了总体中值稳健通道振幅。26个常用通道的稳健通道振幅的统计检验(Wilcoxon符号秩)表明,在去除眼动伪迹之前,所有通道的稳健振幅都与每个记录中值存在显著差异。在去除眼动伪迹和归一化后,Fp2是唯一与记录中值没有显著差异的通道。

图5的中间列绘制了通道i与通道j(其中i≠j)的稳健幅值,分别在去除眼动伪迹前(顶行)、去除眼动伪迹后(中间行)和Huber归一化后(底行)。顶行的图相对于45°角线有些不对称,反映了额叶通道的振幅优势。去除眼动伪迹后,这种关系变得更加对称,具有可见的线性趋势,表明通道振幅存在整体记录水平的均匀缩放。在将通道数据除以记录的Huber平均归一化因子后,通道对幅值之间的线性趋势似乎基本上消除了(底图)。通过对每组通道幅值对拟合线性回归模型得到的调整后的R平方中值,在去除眼动伪迹前为0.465,去除眼动伪迹后为0.481,归一化后为0.015。使用去除眼动伪迹后的数据并分别对不同通道对执行线性拟合发现,在100%的线性拟合中,斜率因子是非零的。这种线性关系解释了通道对幅值的大约一半的可变性,在Huber均值归一化后大部分消失。图5最右边一列显示了记录振幅矢量在二维平面中的t-SNE投影。虽然归一化后投影所占的面积更趋于圆形,但Huber归一化前后的点高度重叠。这强调了交叉研究分析中一致的眼动伪迹去除的重要性。

有些通道在某些频率上始终具有较高的振幅。图6显示了不同频率下的归一化通道幅值。4和8Hz之间的幅值模式类似于图5中去除眼动后的中值通道振幅图,这表明整体通道振幅模式受到4-8Hz范围内幅值差异的高度影响。


图6

源分析

为了研究与特定频段和脑区相关的振幅是否存在系统变化,研究者将测量投影分析(MPA)应用于从脑内IC偶极子计算的频谱图偏差,并使用AAL图谱将结果映射到大脑区域,如图7所示。最初有56574个偶极子,其中26175个在大脑中。


图7

图7的左列显示了不同频率范围内具有高(红色)或低(绿色)局部频谱幅值偏差均值的脑区3D视图:delta([2,4] Hz),theta([4,8] Hz)、alpha([8,14] Hz)和beta([14,31] Hz)。通过12mm标准差3D高斯平滑得到的平均局部振幅偏差在指定频段内至少为5%的脑区显示为红色。振幅偏差最大为-5%的区域用绿色表示。

这些3D视图在alpha振幅上显示出明显的前/后二分情况:大脑后部区域(红色)通常比前部区域(绿色)具有更多的alpha振幅。感觉运动、额中部和颞叶区域的beta活动似乎有所增加。此外,枕下区和小脑区域的活动增加,这可能反映了颈部肌肉的活动。Beta活动沿中央中线结构减弱,包括前扣带回、额上区、上顶叶、后扣带回和内侧枕叶。还可以观察到前扣带回和额上回的theta活动。通过随机置换检验计算基线水平。

研究者调查了这些差异对于整个数据语料库是否具有统计学意义,如图8所示。该图的顶行显示了覆盖有偶极子密度图像的横向脑部视图。采用高斯空间平滑计算偶极子密度,标准差为12 mm。虽然背侧和中线大脑区域具有更高的偶极子密度,但偶极子是分布在整个大脑体积中的。


图8

结论

EEG巨型分析是一个新兴的领域,本研究是对这种分析方法的性能和潜力的初步展示。这些分析的自动化性质和大量数据的汇集将允许对与认知现象相关的大脑动力学进行更系统的研究,同时评估从单一研究或范式中获得的结果的普遍适用性。本文介绍的指标可用于确定预处理中的特定选择是否会改变EEG数据集的一般信号特征,或者某个特定记录是否是该研究的异常值。本研究应用多层次统计建模的巨型分析来评估跨研究的事件相关EEG特征的共性。现阶段和未来的研究可以扩展这些分析,以便对EEG变化的源、分布式皮层源活动以及大规模的脑连接进行更细致的探索。本文使用的原始数据可在https://cancta.net上的DataCatalog上获得。

参考文献:Bigdely-Shamlo, N. , Touryan, J. , Ojeda, A. , Kothe, C. , & Robbins, K.. Automated EEG mega-analysis I: Spectral and amplitude characteristics across studies. NeuroImage, 207.

小伙伴们点个“在看”,加

(星标)关注茗创科技,将第一时间收到精彩内容推送哦~


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

推荐阅读更多精彩内容