贾会宾
华中师范大学心理学院2012年6月7日
其中的KeyStep是数据处理的必要的步骤;ExploratoryStep是为了更好的了解EEGLAB的特点的探索性非必要步骤。
第一章:将数据载入EEGLAB
KeyStep1:StartMATLAB
KeyStep2(Optional):Switchtothedatadirectory(folder)--------也就是把directory转换为数据所在的文件夹。可以通过MATLAB界面上方的“CurrentFolder”或者在commandline通过“cd”。这并不是必须的步骤,在实际的操作中有没有这个步骤没有影响。KeyStep3:StartEEGLAB
在MATLABcommandline键入“eeglab”,并点击“Enter”键。EEGLAB的主窗口将会弹出,如下图所示
KeyStep4:载入数据
目前EEGLAB支持绝大多数的数据类型。
通过File>Importdata可以查看支持的数据类型,当然也包括BP的“***.vhdr”文件。EEGLAB自身的文件格式是“****.set”文件。
在下面,我将使用EEGLABtutorial中自带的“Eeglab_data.set”。载入方式是:File>Loadexistingdataset。
载入后,结果如下:
从这个界面可以看到该数据集一些信息。Channelsperframe为32,也就是channelnumber是32;Framesperepoch为30504,表示的是一个epoch的样本点数目(由于在这个原始数据中还没有分段,所以默认有一个epoch,在这里30504表示的是每个channel样本点总数);Epochs为1,表示的是epoch的数目;Events为1,event表示的该数据集中事件(即诸如刺激呈现、被试反应等称为event)的数目;Samplingrate(Hz)为128;Epochstart(sec)为0,Epochend(sec)为238.305,所以整个数据长度为238.305s;Reference为参考电极;Channellocation为No,这是因为我们还没有对各个channel的头皮坐标进行定位;ICAweights为No,ICAweights表示的是进行ICA分析后得到的weights,此处有两个值“No”和“Yes”。
ExploratoryStep:Scrollingdata
通过如下方式,查看各个channel的波形
通过该界面下的Setting下拉菜单及界面下方的那些选项可以实现一些操作。例如Timerangetodisplay调整该界面显示的时间长度,上图为5s;numberofchannelstodisplay调整一个界面显示的channel数目;并可以选择一个时间段,并“REJECT”。本界面的操作很简单的,与BP相似。
第二章对各个channel进行头皮定位
为了以2-D或3-D的方式绘制EEGscalpmaps,以及对ICA分析后的成分进行源定位,需要数据集包括各个channel的头皮坐标的相关信息。KEYSTEP5:Loadthechannellocations
需要选择Edit>channellocations,得如下界面
在本例中,eeglab开发者给出了该数据集的channellocation文件,需要选择左下方的Readlocations。在弹出的对话框中,选择对应的位置文件。
得到如下结果:
通过上面的页面可以查看各个channel的坐标,通过“Plot2-D”可以查看定位后结果(当然Plot3-D(xyz)可以绘制三维的。
就BP而言,我们没有现成的channellocation文件,那我们需要选择“Lookuplocs”,EEGLAB可以帮助我们自动识别channel位置。在弹出的窗口,点击OK就好了。为了比较一下自动识别后的channellocation是否正确,可以选择“Plot2-D”,并与BP自带的电极位置文件进行比较。我试过,貌似是可以的。还有一种方法(听人说的,没试过):我们依据BP自带的channel坐标的文件,在一个txt文档中写上各个channel的坐标,并将文件的拓展名改为locs。
第三章绘制channel频谱图
最好在进行数据下面的数据处理前,查看各个channel数据的波形,并reject那些不好的数据段。
在本部分,我们将查看数据集的频谱图。这一部分可以为我们的后续操作提供参考依据,并不是必须的步骤。
ExploratoryStep:PlotChannelSpectraandMaps
需要选择Plot>Channelspectraandmaps,将弹出如下窗口
保持默认就好了,点击OK
得到如下结果,其中的各个曲线表示各个的channel,32条;在上方绘制了61022Hz的头皮功率分布。
注意:我们在上述界面默认选择了15%的数据进行分析,因此我们每次得到的结果将有所不同。当取值为100时,当然是无误的结果。
当然,对于已经epoch的数据,我们同样也可以进行如上操作。
第四章预处理工具
1.改变取样率
Tools>Changesamplingrate没有改变)
并不是必须的,但是降低取样率可以节省空间(本例
2.滤波
为了消除线性趋势,我们一般需要高通滤波KEYSTEP6:Removelineartrends
我们需要通过如下选项:Tools>Filterthedata。Filterthedata有两个子菜单BasicFIRfilter和shortIIRfilter。它们是两种不同的滤波方法。
前者界面如下:
后者界面如下:
我们可以在一个还没有epoch的连续数据进行滤波,也可以在epoch后的数据进行滤波。对连续的数据进行滤波可以排除epochboundaries的滤波伪迹。
在本例中,我们选择Tools>Filterthedata>BasicFIRfilter,Loweredge选择1(Hz)。弹出如下结果
在本例中,我们选择Overwriteitinmemory(当然不是必须),这样的话得到结果将覆盖住原始数据集。如果不选择的话,将生成一个新的数据集。
如果我们希望进行带通滤波,需要分别进行高通和低通,在FIR滤波器中直接带通滤波可能产生一些问题(究其原因,开发者没说……)但是IIR滤波器貌似没有这些问题。3.重新设置参考电极
在这个数据集的记录过程中,研究者使用的是将乳突参考。
为了重新设置参考电极,我们需要选择Tools>Re-reference。如下页面将弹出
在这里,我们不更改参考电极。
第五章ExtractingDataepochs
KEYSTEP7:Extractdataepochs
需要选择Tools>Extractepoch,弹出如下窗口
点击“…”弹出如下窗口
在这里我们选择的事件类型是“square”,在本实验中表示的将“square”出现时刻作为每个epoch的时间零点。
在这里,我们不需要把每个epoch的limit进行更改。这里每个epoch时间长度为3s。通常我们将每个epoch设置的时间长些,因为这样可以在低于10Hz的频段范围内进行时频分解。
KEYSTEP8:Removebaselinevalues完成上述操作后,会自动弹出如下窗口
在这里基线为[-10000],点击OK。
第六章数据叠加平均
1.PlottingtheERPdataonasingleaxiswithscalpmapsExploratoryStep:Plottingall-channelERPs
在这一步,我们将绘制所有epochs的叠加平均(ERP)和某一特定潜伏期的ERPscalpmaps。我们需要选择
Plot>ChannelERP>Withscalpmaps
弹出如下窗口
保持默认值就好,点击Ok
在上图结果中,每条曲线对应于各个channel,上方的地形图是430ms时刻的平均电压地形图分布(当我们使用默认值的时候,eeglab将默认绘制ERP方差最大时刻的地形图,在本例是430ms)。
2.PlottingERPtracesinatopographicarray
ExploratoryStep:PlottingERPsinaTopographicMap
在这里我们希望按照各个channel的头皮分布分别绘制各个channel的ERP。需要选择Plot>ChannelERPs>Inscalp/rect.array
弹出如下窗口
点击Ok,弹出如下窗口
当我们双击其中的任何一个channel的时候,会弹出相应的channelERP。3.Plottingaseriesof2-DERPscalpmaps
在这里,我们将绘制一系列的2-Dscalpmaps,其中每个图表明的是一个特定的潜伏期的电压分布。
ExploratoryStep:Plottingaseriesof2-DERPScalpMaps选择Plot>ERPmapseries>In2-D,弹出如下窗口
假设我们感兴趣的是潜伏期0100200300400500ms的电压分布。我们需要在上面的对话框第一行键入0:100:500。点击Ok
第七章选择数据的epochs,并比较
1.Selectingeventsandepochsfortwoconditions
为了比较一个被试两种条件下的ERP,我们需要首先为两种条件各创建一个dataset。在这个实验中,一半的目标刺激呈现在位置1,一半的目标刺激呈现在位置2。ExploratoryStep:SelectingEventsandEpochsforTwoConditions需要选择Edit>Selectepochsorevents弹出如下窗口
点击“Stimulusposition”
我们分别键入1和2,可创建两个dataset,并分别将其命名为“Square,Pos.1”和“Square,Pos.2”。
在EEGLAB中,另外一个可用于数据集的子集选择的函数是pop_select.m,这个函数可以通过如下方式调用Edit>Selectdata弹出如下窗口
在上述窗口,键入相应参数。例如
通过上述操作,我们将提取每个epoch的-0.5至1s,去除了第2、3、4个epoch。2.ComputingGrandMeanERPs
ExploratoryStep:ComputingGrandMeanERPs
我们需要选择Plot>Sum/CompareERPs,弹出如下窗口
我们在第一行中间的文本框,输入3和4(它们分别表示上述两个生成的数据集。3和4分别是它们在EEGLABGUI界面上的编号)。右方的avg表示总平均ERP,std表示standarddeviation,allERPs表示输出所有数据集的ERP。
在highlightsignificantregions键入0.05,表示t-testsignificanceprobability(p)threshold。
点击OK
双击FPz
3.ComparingERPsintwoconditions
ExploratoryStep:ComparingERPsinTwoConditions
我们需要选择Plot>Sum/CompareERPs弹出如下窗口
键入参数
表示比较上述两个数据集,低通滤波30Hz,并绘制标题,“Position1-2”点击OK
双击FPz
第八章绘制ERPimages
为了对ERP效应有一个更好理解,EEGLAB有一个比较有特色的功能,即绘制ERPimage。这个ERPimage是一个2-Dimage,其中的横轴是每个epoch的时刻值,纵轴是各个epoch的编号,而该图像中的每一点表示相应的epoch的相应时刻的电压值。至于纵轴上的各个epoch的顺序,EEGLAB默认是按照它们在实验中出现的顺序进行排序。当然,研究者可以依据自己的个人兴趣,对各个epoch的纵轴顺序重新排序(例如,依据subjectreactiontime,alpha-phaseatstimulusonset)------理论上来说,排序的方式可以有很多种。当然每种排序所能提供的信息不尽相同,需要针对自己的研究来具体选择。
同时需要注意的是,有些时候如果不仔细考虑会很容易对结果进行错误解释。例如,利用某一频段的相位进行排序,可能会掩盖相同数据的不同频段的振荡活动。1.选择需要绘制的channel
假设我们想选择alpha频段(大约为10Hz)功率最大的channel。需要选择Plot>Channelspectraandmaps
可以看出alpha频段的能量主要集中在枕叶。
鉴于此,我们选择POz(27号)电极
2.使用pop_erpimage()绘制ERPimagesExploratoryStep:ViewingaChannelERP选择Plot>ChannelERPimage弹出如下窗口
选择channel27,smoothing1(表示的是在临近的epochs进行平滑绘图的的时候,所取的值),其他默认,点击Ok弹出如下结果
上图的最上方为选择的电极的头皮位置
中间的图即为ERPimage,下图为该电极的ERP。ExploratoryStep:PlottingaSmoothedERP
在这里我们将smoothing设置为10,它表示在相邻的10个epoch平滑(我的理解),弹出如下窗口
可以看出,上述结果可以更好的显示alpha频段的振荡活动。3.对ERPimage中的trial进行排序
在上面的ERPimage中,从下到上,EEGLAB默认按照trial在实验过程中出现的顺序排列。当然也可以按照其他的变量进行排序。在下面我们使用的变量的是被试的反应时间,具体来说是事件‘rt’的潜伏期(latency)。
ExploratoryStep:SortingTrialsinanERPImage选择Plot>ChannelERPimage
双击,弹出如下窗口
选择latency,点击OK
双击,弹出如下窗口
选择rt,点击Ok。
在“Eventtimerange”下方的文本框,键入-200800(不是很理解这一选项)。点击Ok。弹出如下窗口
上图中的黑线表示事件“rt”的“latency”。
如果不对“Eventtimerange”下方的文本框进行设置的话,弹出如下结果
看不出两种结果之间的不同….
4.使用频谱选项绘制ERPimage
ExploratoryStep:SortingTrialsinanERPbyPhaseValue我们键入的参数如下图
点击OK,弹出如下结果
(ITC)ExploratoryStep:Inter-TrialCoherenceCoherence(ITC)
为了在ERPimagefigure中绘制ITC,我们需要在pop_erpimage中设置如下参数
这样我们便可以分析[911]Hz,并进行显著性检验(显著性水平0.01)。点击Ok
上图的最上方子图的ERPimage。
倒数第二个子图是ERSP(EventRelatedSpectralPower),表示的是功率的平均变化(单位为dB)。在这个子图中,曲线没有超出蓝色区域,表示在选择的频率10.13Hz(见右下角)相对基线水平(25.93dB),功率没有显著变化。
最下方的子图的ITC。它表示的相对刺激呈现而言,相位同步化的程度。10.13Hz表示选择分析的频段。可以看出在300Hz附近相位同步化显著增强。
第九章使用ICA分解数据
KEYSTEP9:CalculateICAComponents需要选择Tools>RunICA,弹出
这里默认的算法是runica,按照这个默认的就好,点击Ok。耐心等待,速度一般较慢。1.绘制2-DComponentScalpMaps
我们需要选择Plot>Componentmaps>In2-D。弹出如下窗口
鉴于我们想绘制1:12个成分,我们需要将“Componentnumbers”后面的文本框改为1:12。点击Ok
凭借这个我们可以了解各个成分的头皮分布。EEGLAB有一个很有用的插件ADJUST,它可以半自动的去除伪迹成分。ADJUST需要自己安装,不是EEGLAB默认携带的。步骤如下,选择Tools>ADJUST。弹出如下窗口
鉴于我们已经运行ICA了,选择RunADJUST。弹出如下结果
从上图可以看出,ADJUST认为第7、10、27成分是伪迹。我们双击出如下结果
,弹
单凭有初步经验的人也可以看出,此成分明显是一个眨眼伪迹。ADJUST的开发者说,ADJUST的准确率相对专家可达到95%。所以相对我这种没有太多经验的人来说,姑且就认为ADJUST的判断正确吧。
下一步是删去第7、10、27成分,我们需要选择Tools>Removecomponents
在右方的文本框,键入71027,点击Ok。
然后选择“Accept”。
第十章WorkingwithICAComponents
1.绘制componentspectraandmaps
在本部分,我们希望了解哪些成分对特定的频段贡献最大。我们需要选择plot>Componentspectraandmaps。弹出如下窗口:
表示是需要随机选择的数据的百分比。百分比小的话,会节省时间;
百分比越大,结果当然越准确。在这里我们选择100。其余选项默认,点击Ok。
在下面,为了估计各个成分对某一个特定的channel的贡献,我们将介绍一个更为准确的方法(至于为什么是一个更为准确的做法,开发者只是说了“fortechnicalreasons”)。
假设我们需要考查27号channel。我们需要
右方的文本框,键入
27。
uncheck。点击OK。
,右方的
同时commandline出现如下结果
“Component1percentvarianceaccountedfor:1.16Component2percentvarianceaccountedfor:4.67Component3percentvarianceaccountedfor:13.79Component4percentvarianceaccountedfor:32.87Component5percentvarianceaccountedfor:-0.96Component6percentvarianceaccountedfor:4.91Component7percentvarianceaccountedfor:0.18Component8percentvarianceaccountedfor:5.05Component9percentvarianceaccountedfor:27.82Component10percentvarianceaccountedfor:1.08Component11percentvarianceaccountedfor:40.35Component12percentvarianceaccountedfor:3.10Component13percentvarianceaccountedfor:2.56Component14percentvarianceaccountedfor:Component15percentvarianceaccountedfor:Component16percentvarianceaccountedfor:
3.4.032.86
Component17percentvarianceaccountedfor:1.Component18percentvarianceaccountedfor:-1.68Component19percentvarianceaccountedfor:0.02Component20percentvarianceaccountedfor:1.11Component21percentvarianceaccountedfor:0.12Component22percentvarianceaccountedfor:-0.11Component23percentvarianceaccountedfor:-0.22Component24percentvarianceaccountedfor:0.38Component25percentvarianceaccountedfor:-0.03
Component26percentvarianceaccountedfor:-0.41Component27percentvarianceaccountedfor:0.12Component28percentvarianceaccountedfor:Component29percentvarianceaccountedfor:2.绘制成分ERPs
0.000.06”
我们需要选择Plot>ComponentERPs>Inrectangulararray。
点击Ok。
为了查看各个成分的ERP,双击相应的成分。
为了比较多个dataset的成分的ERP,我们可以采用与channelERPs相似的方法,选择Plot>Sum/Comparecomp.ERPs。事实上,对话框也很相似,如下:
不再详述………
3.绘制成分ERP贡献
为了得到成分的ERP对数据的ERP的贡献,我们需要选择Plot>ComponentERPs>withcomponentmaps。
点击Ok就好了。
上图中的黑色包络线是各个时间点的所有channel的最大值和最小值,彩线是各个成分的的ERP。
现在我们希望了解对200-500ms数据ERP贡献最大的成分。
右方的文本框改为200500。
在commandline返回如下结果:
“Comparingmaximumprojectionsforcomponents:IC1maximummeanpowerofback-projection:80.1385IC2maximummeanpowerofback-projection:1.73828IC3maximummeanpowerofback-projection:1.69376IC4maximummeanpowerofback-projection:2.55062IC5maximummeanpowerofback-projection:6.99448IC6maximummeanpowerofback-projection:9.00769IC7maximummeanpowerofback-projection:3.06103IC8maximummeanpowerofback-projection:6.6IC9maximummeanpowerofback-projection:2.88371IC10maximummeanpowerofback-projection:4.23078IC11maximummeanpowerofback-projection:3.42227IC12maximummeanpowerofback-projection:4.721IC13maximummeanpowerofback-projection:1.90804IC14maximummeanpowerofback-projection:1.8131IC15maximummeanpowerofback-projection:0.853632IC16maximummeanpowerofback-projection:0.683555IC17maximummeanpowerofback-projection:5.9612IC18maximummeanpowerofback-projection:0.86812IC19maximummeanpowerofback-projection:0.73882IC20maximummeanpowerofback-projection:0.61765IC21maximummeanpowerofback-projection:0.287868IC22maximummeanpowerofback-projection:0.277423IC23maximummeanpowerofback-projection:0.173851IC24maximummeanpowerofback-projection:0.0984656IC25maximummeanpowerofback-projection:0.2357IC26maximummeanpowerofback-projection:0.0806276IC27maximummeanpowerofback-projection:0.100677IC28maximummeanpowerofback-projection:0.0798501IC29maximummeanpowerofback-projection:0.0388474intheinterval203msto500ms.
Plottingenvelopesof7componentprojections.Topomapswillshowcomponents:176
8
1015withmaxvarattimes(ms):2812
313
352
398
430epochframes:319320323328334
338
Componentsortvarininterval:80.141.741.692.556.999.013.06Summedcomponent'ppaf'ininterval[203.125500]ms:88.98%Plotlimits(sec,sec,uV,uV)[-1,1.99219,-18.7508,22.2963]”
12453341
第十一章Time/Frequencydecomposition
1.DecomposingchanneldataKeyStep10:ERSPandITC
为了检测ERSP(event-relatedspectralperturbation)和ITC(inter-trialcoherence),我们需要选择Plot>Timefrequencytransforms>Channeltime-frequency。弹出如下对话框
Channelnumber选择14(Cz),Bootstrapsignificancelevel选择0.01,其它默认,点击Ok。
得到的结果分为两个子图,上方是ERSP子图,下方是ITC子图。
ERSP子图:该子图左方的panel是基线的平均功率谱,而各个时间点的ERSP包络线。而topimage是相对基线,每个时间点、每个频率的频谱功率改变(event-relatedchangesinspectralpower(frompre-stimulusbaseline)ateachtimeduringtheepochandateachfrequency(<50Hz))。
ITC子图:我们在前面遇到过ITC。AsignificantITCindicatesthattheEEGactivityatagiventimeandfrequencyinsingletrialsbecomesphase-locked(notphase-randomwithrespecttothetime-lockingexperimentalevent).可以看出,在本图中,大约在10Hz左右出现了significantinter-trialcoherence。
2.Computingcomponenttime/frequencytransforms
在本部分,我们将对成分进行时频分解,因为各个成分可能更直接反应大脑的EEGsource。
我们需要选择Plot>Time/frequencytransforms>Componenttime-frequency。
Componentnumber选择10,Subepochtimelimits选择-5001000,选择UseFFT,Bootstrapsignificancelevel选择0.01。点击Ok
备注:之所以选择FFT是因为,相对wavelets而言,它可以计算更低的频率。3.计算成分的cross-coherences
为了决定两个成分活动之间同步化的程度,我们可以绘制它们之间的event-relatedcross-coherence。我们需要选择Plot>Time-frequencytransforms>Componentcross-coherence。
我们的参数设置如下:
点击Ok。
本结果的上方子图是两个成分的coherencemagnitude,而下方的子图表示的是:在上方各个coherencemagnitude显著的区域,两个成分的相位差。
当然我们同样可以绘制各个channel之间的coherence,需要选择Plot>Time-frequencytransforms>Channelcross-coherence。
补充:
一、使用DIPFIT.2对成分进行等价偶极子定位
在使用DIPFIT之前需要载入channellocation,并进行ICA分解。在这里我们将使用另一个数据集eeglab_dipole.set。
进行偶极子定位时,需要三个步骤,分别是(1)Settingmodelandpreferences;(2)Gridscanning;(3)Non-linearinteractivefitting。1.SettingupDIPFITmodelandpreferences
我们需要选择Tools>LocatedipolesusingDIPFT>Headmodelandsettings。
Headmodel选项为我们提供了四个可供选择的headmodel。在这里我们选择BEM。现在说一下Co-registerchan.locs.Withheadmodel选项。如果我们使用的电极位置是10-20系统的,我们不需要co-register。需要选择“NoCo-Reg”。否则的话,我们需要选择ManualCo-Reg。
注意:在eeglab的wikitutorial中,作者说针对样本数据集eeglab_dipole.set我们不需要co-register。但是我发现这样做的话,会得到错误的结果(经过与UCSD的老师联系,这确实是一个错误!),所以在这里我选择manualco-reg。点击ManualCo-Reg。
选择。弹出如下对话框
点击OK。返回coregister对话框,点击OK。返回Dipfitsetting对话框,点击Ok。
2.Initialfitting
这一步的fitting的结果不精确,较为粗糙;但是可以为后一步的精确fitting提供一个起点。我们需要选择Tools>LocatedipolesusingDIPFT>Coarsefit(gridscan)。
点击Ok。
3.Non-linearinteractivefitting
我们需要选择Tools>LocatedipolesusingDIPFT>Autofit(coarsefit,finefit,plot)
默认fit所有的71个成分,但是这通常会很费时间。为了简单起见,我们选择前10个独
立成分。
后面的文本框选择1:10。点击OK。
为了查看fitting后的结果,eeglab提供了两种方法:3-D的和2-D的。(1)Tools>LocatedipolesusingDIPFIT>Plotcomponentdipoles
点击OK。
在上图中可以很容易的从各个视角查看各个偶极子。也可以通过plotone/plotall和Keep|Next查看哪些偶极子。很简单的,不再详述。(2)Plot>Componentmaps>In2-D
在这里我们绘制前20个成分。,打上勾。点击Ok。
在上图中可以看到,第三个成分的偶极子位于中线位置,scalpmap也几乎是对称分布的,所以在头皮左右两侧各绘制一个偶极子似乎更为合适。在这里选择Tools>LocatedipolesusingDIPFIT>Finefit(iterative)。
修改参数设置,
选择。程序运行结束后,选择。
成功了!Manualdipolefit对话框点击Ok。
二、EEG构架数组
各个数据集在上述各个操作阶段产生的数据存储在一个EEG构架数组。在命令行,键入EEG>>EEGEEG=
setname:'ap82'
filename:'eeglab_dipole.set'
filepath:'C:\\Users\\Administrator\\Desktop\\'subject:''group:''condition:''session:[]
comments:[10x61char]
nbchan:69(channel的数目)
trials:5
pnts:1000(每个trial的取样点)srate:250(取样率)
xmin:-3(每个trial的开始时间点)xmax:0.9960(每个trial的结束时间点)times:[1x1000double]
data:[69x1000x5single](各个channel的数据)icaact:[71x1000x5single]icawinv:[69x71double]icasphere:[71x69double]
icaweights:[71x71double](进行成分分析时得到的weight)icachansind:[1x69double]
chanlocs:[1x69struct]urchanlocs:[]
chaninfo:[1x1struct]
ref:'common'event:[1x9struct]urevent:[1x9struct]
eventdescription:{[]
[]
''''''''''''''''''''''''}
(各个成分的数据)
epoch:[1x5struct]
epochdescription:{}
reject:[1x1struct]stats:[1x1struct]specdata:[]specicaact:[]splinefile:''
icasplinefile:''
dipfit:[1x1struct](上面dipfit过程中产生的数据)history:[1x26char]saved:'no'etc:[]datfile:[]
因篇幅问题不能全部显示,请点此查看更多更全内容
Copyright © 2019- sceh.cn 版权所有 湘ICP备2023017654号-4
违法及侵权请联系:TEL:199 1889 7713 E-MAIL:2724546146@qq.com
本站由北京市万商天勤律师事务所王兴未律师提供法律服务