CN103767707A - Blood sample level dependence functional magnetic resonance signal fluctuating frequency clustering analysis method - Google Patents

Abstract

The invention discloses a blood sample level dependence functional magnetic resonance signal fluctuating frequency clustering analysis method which is characterized by comprising the following steps: step 1, pretreatment analysis is performed on resting state data collected from a magnetic resonance machine; step 2: Hilbert-Huang transform is adopted for decomposing the pretreated signals, so as to obtain intrinsic mode functions with different frequencies; the Hilbert-Huang transform comprises two steps of empirical mode decomposition and Hilbert transform; step 3, quantitative indexes of different frequency components are calculated; step 4, clustering of different voxels in a brain is performed by taking the quantitative indexes in the step 3, as a basis for classification and statistical analysis to the different indexes is performed. The method solves the problem of insufficiency of data drive frequency analysis methods in the fMRI research field, an iconography analytical method is provided for basic science and clinical researches and meanwhile new quantitative indexes and physiological markers are provided for clinical disease diagnosis and pathological mechanisms.

Description

The clustering method of blood sample horizontal dependency functional magnetic resonance signal concussion frequency

Technical field

The invention belongs to Medical Image Processing and analysis technical field, specifically, relate to the clustering method of a kind of blood sample horizontal dependency functional magnetic resonance signal concussion frequency.

Background technology

Blood oxygenation level dependent canbedetected functional mri can nondestructively detect and locate brain domain, is a new technique of neuroimaging development.It not only can reflect brain function and physiological change, also can carry out brain function positioning integration neuro anatomy and function information, the reflection neuromechanism relevant to various diseases or Cognitive task directly perceived.

The discharge activities of the axoneuron that the various Physiological Psychology activities of body are all followed.The point activity of neuron collective causes nerve and the biochemical variation of surrounding tissue blood flow.Neuron activity need consume a large amount of energy, and cerebral tissue provides these energy by Repiration, and Repiration needs oxygen consumed, and therefore the increase of cerebration will be followed the increase of local oxygen consumption.In order to provide oxygen to divide to enlivening Nao district, regional flow's perfusion also can increase, and the increase of regional flow's perfusion can bring more HbO2 Oxyhemoglobin.But the increase of local brain tissue blood perfusion is different with the ratio that blood oxygen consumption increases, and blood flow increase obviously exceedes the increase of oxygen consumption, this species diversity causes the venous blood oxygenate hemoglobin of brain active region to increase, and deoxyhemoglobin ratio reduces relatively.Because deoxyhemoglobin is paramagnet, can cause that T2 shortens, in the time of its lowering of concentration, cause T2* and T2 to extend, on T2*WI, signal strengthens, and this is Blood oxygenation level dependent canbedetected.The principle imaging that BOLD-fMRI strengthens according to this magnetic susceptibility contrast.

Than brain electricity (EEG), x-ray tomography imaging (CT) and positron emission tomography (PET), BOLD-fMRI has high, the cold advantage of spatial resolution, can nondestructively observe structure and the metabolic activity at any position, brain deep.But compared to EEG signal, BOLD-fMRI temporal resolution is lower, hematodinamics signal is also considered to low frequency signal conventionally, and BOLD-fMRI research field lacks the frequency analysis technique of the data-driven that is similar to EEG data analysis field at present.

Multinomial research is found, the function corresponding to E.E.G of different frequency, and the propagation between different information Nao district may realize by the E.E.G of different frequency.BOLD-fMRI studies and shows, the blood oxygen level variable signal of different frequency may reflect different brain functioies, and different diseases also affects the BOLD activity of characteristic frequency specifically.These documents show, the frequency of research cerebration is significant.

But in the past clinical and basic scientific research all ignored BOLD signal frequency information and often observation signal amplitude along with the variation of different task.Whether the BOLD signal amplitude of for example visual cortex raises during compared to quiescent condition in the time accepting visual stimulus, also change but seldom there is research to pay close attention to its frequency.Along with the increase of tranquillization state BOLD-fMRI in various diseases and basic physiology research, it is found that the various quantizating index of tranquillization state BOLD signal also can reflect the change of brain function.

Whether have different diseases or task also caused the change of BOLD signal frequency? answer the frequency that first this problem will portray BOLD activity in Different brain region under Normal Subjects state.The research of most BOLD frequency characteristics is all with the band filter that artificially sets in advance parameter, BOLD signal to be processed, and this does not only consider the nonlinear characteristic of BOLD signal non-stationary, and the division of how to confirm frequency range is also a problem.Frequency range is divided to such an extent that incorrectly may cause the unstable of result, and linear filter also makes BOLD distorted signals simultaneously.

Summary of the invention

The technical problem to be solved in the present invention is to overcome above-mentioned defect, the clustering method of a kind of blood sample horizontal dependency functional magnetic resonance signal concussion frequency is provided, solve the problem that current fMRI research field lacks the frequency analysis method of data-driven, for basic science and clinical research provide new imaging Analysis method, simultaneously for the research of clinical disease diagnosis and pathomechanism provides new quantizating index and physiology mark.

For addressing the above problem, the technical solution adopted in the present invention is:

The clustering method of blood sample horizontal dependency functional magnetic resonance signal concussion frequency, is characterized in that: comprise the following steps:

Step 1: the tranquillization state data that collect from magnetic resonance machine are carried out to Preprocessing;

Step 2: adopt Hilbert-Huang transform to decompose pretreated signal, obtain the intrinsic mode function of different frequency; Hilbert-Huang transform comprises empirical mode decomposition and Hilbert transform two steps;

Step 3: the quantizating index that calculates different frequency composition;

Step 4: take the quantizating index of step 3 as classification foundation, by different voxel clusters in brain, the above different indexs of statistical analysis.

Improve described step 1 as one) in, be BOLD-fMRI signal data to the tranquillization state data that collect from magnetic resonance machine, the specific implementation step of BOLD-fMRI signal data being carried out to Preprocessing is:

1), utilize SPM software based on MATLAB platform, the tranquillization state data that collect are carried out to time difference rectification;

2), the time is corrected to data later and carry out the moving alignment of head;

3), the data after correct moving alignment are carried out Spatial normalization;

4), the data after Spatial normalization are gone to linear drift and time domain standardization.

Improve described step 2 as one) in, adopt Hilbert-Huang transform to decompose pretreated signal, obtain the composition of different frequency, i.e. intrinsic mode function, specific implementation step is:

1), adopt empirical mode decomposition signal is decomposed, obtain some intrinsic mode functions;

2), the intrinsic mode function obtaining after empirical mode decomposition is carried out to Hilbert transform.

Improve as one, adopt Hilbert-Huang transform to decompose pretreated signal, wherein in step 1), adopt empirical mode decomposition to decompose signal, the specific implementation step that obtains some intrinsic mode functions is:

Step 1.1), find all extreme points of primary signal x (t);

Step 1.2), between all maximum points, interpolation obtains the coenvelope eup (t) of x (t), between all minimum points, interpolation obtains the lower envelope elow (t) of x (t);

Step 1.3), calculate eup (t) and average time series emean (t)=[elow (t)+eup (t)]/2 of elow (t);

Step 1.4), x (t) deducts average time series and obtains concussion pattern a: r (t)=x (t)-emean (t);

Step 1.5), whether the disconnected concussion pattern r (t) of judgement meet feature or the monotonous curve of intrinsic mode function, if meet IMF feature, should " disconnected concussion pattern " become an intrinsic mode function, IMF _i(t)=r (t); Monotonous curve no r (t) becomes residual error, and EMD process finishes, and repeats above step until obtain an IMF otherwise establish x (t)=r (t); IMF is characterized as lower envelope symmetry, and the number of times of IMF curve zero crossing equates with extreme point number or differs at most 1.

Improve as one, described employing Hilbert-Huang transform decomposes pretreated signal, wherein step 2) in the intrinsic mode function obtaining after empirical mode decomposition is carried out to Hilbert transform specific implementation step be:

Step 2.1), to establish certain IMF be that its Hilbert transform Y (t) of X (t) is

wherein P is Cauchy's principal value $P{\∫}_a^{b}f\left(x\right)\mathrm{dx}={\lim }_{e\→\∞}({\∫}_a^{x-e}f\left(x\right)\mathrm{dx}+{\∫}_{x+e}^{b}f\left(x\right)\mathrm{dx});$

Step 2.2), by X (t) and Y(t) merge into a complex signal Z(t), make Z (t)=X (t)+iY (t)=a (t) ^{ei θ (t)}, wherein a (t)=[X2 (t)+Y2 (t)] ^1/2the envelope of complex signal, Q(t)=tan ^-1(Y/X) be the instantaneous phase of complex signal; The instantaneous frequency of hilbert spectrum is defined as

Step 2.3), by all intrinsic mode functions are done to Hilbert transform, the primary signal x (t) after empirical mode decomposition can be expressed as $x\left(t\right)={\mathrm{\Σ}}_{i=1}^K{a}_i\left(t\right)\exp \left\{i{\∫w}_j\left(t\right)\mathrm{dt}\right\}.$

Improve step 3 as one) in, calculate different frequency composition, i.e. the quantizating index of intrinsic mode function, specific implementation step is:

Step 3.1), calculate the energy of each intrinsic mode function;

Step 3.2), calculate the Hilbert weighted frequency of each intrinsic mode function and the Hilbert weighted average frequency of primary signal;

Step 3.3), calculate other indexs, as locally coherence.

Improve described step 3.1 as one) in, the circular that calculates the energy of each intrinsic mode function is,

Improve described step 3.2 as one) in, the specific implementation step of calculating the Hilbert weighted frequency of each intrinsic mode function and the Hilbert weighted average frequency of primary signal is:

Step 3.2.1), to establish primary signal x (t) length be N time point, empirical mode decomposition obtains K intrinsic mode function, the Hilbert weighted frequency of j intrinsic mode function is

wherein a _jand w (i) _j(i) envelope and the instantaneous frequency of the moment i in being respectively;

Step 3.2.2), the Hilbert weighted average frequency of primary signal x (t) by all intrinsic mode functions HWFj}, j=1 ..., K, obtains through the weighting of each intrinsic mode function envelope norm, $\mathrm{HWMF}=\frac{{\mathrm{\Σ}}_{i=1}^N\left|\right|a_j\left|\right|{\mathrm{HWF}}_j}{{\mathrm{\Σ}}_{i=1}^N\left|\right|a_j\left|\right|}.$

Improve step 3.3 as one) in, specific implementation step is:

The all voxels of full brain are calculated to its certain intrinsic mode function and the ReHo of the identical intrinsic mode function of adjacent voxels around one by one,

wherein W is exactly ReHo value mistake! Do not find Reference source., scope 0 to 1,

wherein r _ijj voxel mistake! Do not find Reference source.In the order of i time point;

r _iaverage, n mistake! Do not find Reference source.It is the time point number of this intrinsic mode function, total number of this voxel of k and adjacent voxels, if adjacent voxels is defined as the k=27 of voxel contacting with this voxel face, limit or fixed point, if be defined as face or the tactile k=19 of voxel of edge joint, if be defined as the k=7 of voxel of face contact.

Improve as one, in described step 4, take above-mentioned quantizating index as classification foundation, by different voxel clusters in brain, and the specific implementation step of the above different indexs of statistical analysis is:

Step 4.1), the indices obtaining using step 3 or wherein some as feature, carries out K-means cluster to all voxels of full brain;

Step 4.2), utilize the mind map of the indices that SPM software obtains step 3, obtain final result not carrying out on the same group statistical between object of study.

Improve step 4.1 as one), the indices obtaining using step 3 or wherein some as feature, the concrete steps of all voxels of full brain being carried out to K-means cluster are:

Step 4.1.1), from all voxels of full brain, choose at random the characteristic index barycenter of K voxel;

Step 4.1.2), measure the distance of its characteristic index to each barycenter to remaining each voxel, and it is grouped into the class of nearest barycenter, make each voxel be grouped into a certain apoplexy due to endogenous wind until travel through all voxels of full brain;

Step 4.1.3), recalculate the barycenter (average of such characteristic index) of each class having obtained;

Step 4.1.4), iteration 4.1.2～4.1.3 walks until new barycenter equates or be less than assign thresholds with former barycenter, and algorithm finishes.

Owing to having adopted technique scheme, compared with prior art, the present invention includes four steps: first the BOLD-fMRI data that collect from magnetic resonance machine are carried out to Preprocessing; Then adopt advanced data-driven method---Hilbert-Huang transform, by pretreated BOLD signal adaptive be decomposed into the IMF of different frequency range; Calculate subsequently that energy that HHT decomposes the each IMF obtaining, Hilbert weighted frequency, locally coherence etc. are all kinds of to be estimated or quantizating index; Finally using these estimate or quantizating index as classification foundation by the voxel cluster of different parts in brain, which or estimate or quantizating index there are differences by above-mentioned in brain district between statistical method comparison different crowd or different conditions, thereby find patient's central neuropathy region, or reflect the change of brain function activity under normal person's different conditions.The method that we provide for find various nervus centralis diseases focus, drug target, probe into sacred disease pathogenesis new method be provided, the support of methodology, iconography is also provided for the basic research of cognitive neuroscience, neuro physiology, has there is important scientific research and clinical value.

The invention solves the problem that current fMRI research field lacks the frequency analysis method of data-driven, for basic science and clinical research provide new imaging Analysis method, simultaneously for the research of clinical disease diagnosis and pathomechanism provides new quantizating index and physiology mark.

Below in conjunction with the drawings and specific embodiments, the invention will be further described simultaneously.

Accompanying drawing explanation

Fig. 1 is the realization flow schematic diagram of the FCA method of an embodiment of the present invention;

Fig. 2 is the EMD realization flow figure of an embodiment of the present invention;

Fig. 3 uses the normal tested tranquillization state fMRI data acquired results of FCA methods analyst in an embodiment of the present invention.

The specific embodiment

Embodiment:

The clustering method of blood sample horizontal dependency functional magnetic resonance signal concussion frequency, comprises the following steps:

Step 1, carries out Preprocessing to the tranquillization state data that collect from magnetic resonance machine.

Step 2, adopts Hilbert-Huang transform, and Hilbert-Huang transform comprises empirical mode decomposition (EMD) and Hilbert transform (HT) two steps, and pretreated signal is decomposed, and obtains the intrinsic mode function of different frequency.

Step 3, the various quantizating index of calculating different frequency intrinsic mode function (IMF), as energy, Hilbert weighted frequency (HWF), locally coherence () etc.

Step 4, take above-mentioned quantizating index as classification foundation, by different voxel clusters in brain.The above different indexs of statistical analysis.

In the present embodiment, step 1) specific implementation step that the BOLD-fMRI signal data to collecting from magnetic resonance machine collecting from magnetic resonance machine is carried out to Preprocessing is:

Step 1.1, utilizes the SPM software based on MATLAB platform, and the tranquillization state data that collect are carried out to time difference rectification.

Step 1.2, corrects data later to the time and carries out the moving alignment of head.

Step 1.3, the data after correct moving alignment are carried out Spatial normalization.

Step 1.4, goes linear drift and time domain standardization to the data after Spatial normalization.

In the present embodiment, to magnetic resonance acquisition to data carry out time difference while correcting, the difference of correcting exactly in 1 volume acquisition time between layers.

The data that time corrected are carried out to the moving alignment of head, exactly each two field picture of a sequence is all carried out to registration with the first two field picture of this sequence, be registrated under the same coordinate system, moving to correct head, after alignment, the data after correct moving rectification are carried out hand inspection again, spend and just exclude, will not analyze if translation and rotation exceed respectively 1 millimeter and 1;

Data after correct moving alignment are carried out Spatial normalization, exactly tested brain registration are normalized in the brain template of standard; When data after Spatial normalization are carried out to space smoothing, the data after the gaussian kernel function smoothing processing standardization of 6 millimeters of full width at half maximum of employing.

In the present embodiment, the data after correct moving alignment are carried out Spatial normalization, and the specific implementation step tested brain registration being normalized in the brain template of standard is:

Data after correct moving alignment are utilized the affine transformation the average image of 12 parameters to be registrated to MNI standard form to lay equal stress on and be cut to the voxel of 3x3x3 cubic millimeter, and by MNI coordinate transform to Talairach coordinate system.

Go linear drift and the standardized specific implementation step of time domain to be to the data after Spatial normalization:

By linear regression, the composition of linear growth or reduction is removed, the time series of each voxel is deducted to the time domain average of self the time domain standard deviation divided by self.

In the present embodiment, step 2) in utilization is write based on MATLAB platform program software, adopt Hilbert-Huang transform to decompose pretreated signal, the specific implementation step that obtains the composition of different frequency is:

Step 2.1, adopts EMD to decompose signal, obtains some IMF;

Step 2.2, the IMF obtaining after EMD is decomposed carries out Hilbert conversion.

In the present embodiment, step 2.1) adopt EMD to decompose signal, the specific implementation step that obtains some IMF is:

Step 2.1.1, finds all extreme points of primary signal x (t);

Step 2.1.2, between all maximum points, interpolation obtains the coenvelope eup (t) of x (t), and between all minimum points, interpolation obtains the lower envelope elow (t) of x (t);

Step 2.1.3, the average time series of calculating eup (t) and elow (t)

emean(t)=[elow(t)+eup(t)]/2;

Step 2.1.4, x (t) deducts average time series and obtains one " concussion pattern "

r(t)=x(t)-emean(t);

Step 2.1.5, whether judgement " disconnected concussion pattern " r (t) meets feature or the monotonous curve of IMF, should " disconnected concussion pattern " become an IMF, IMF if meet IMF feature _i(t)=r (t), if monotonous curve no r (t) become residual error, EMD process finishes, and repeats above step until obtain an IMF otherwise establish x (t)=r (t).IMF is characterized as lower envelope symmetry, and the number of times of IMF curve zero crossing equates with extreme point number or differs at most 1.

In the present embodiment, step 2) in EMD is decomposed after the IMF that the obtains specific implementation step of carrying out Hilbert conversion be:

Step 2.2.1, establishing certain IMF is X (t) mistake! Do not find Reference source.Its Hilbert conversion Y (t) mistake! Do not find Reference source.For mistake! Do not find Reference source.

wherein P is Cauchy's principal value $P{\∫}_a^{b}f\left(x\right)\mathrm{dx}={\lim }_{e\→\∞}({\∫}_a^{x-e}f\left(x\right)\mathrm{dx}+{\∫}_{x+e}^{b}f\left(x\right)\mathrm{dx})$ A mistake! Do not find Reference source.；

Step 2.2.2, merges into a complex signal Z (t) by X (t) and Y (t), order ^{a mistake! Do not find Reference source.}z (t)=X (t)+iY (t)=a (t) ^{ei θ (t)}, wherein a (t)=[X2 (t)+Y2 (t)] ^1/2the envelope of complex signal, Q (t)=tan ^-1(Y/X) be the instantaneous phase of complex signal.The instantaneous frequency of hilbert spectrum is defined as

Step 2.2.3, by all IMF are done to Hilbert transform, primary signal x (t) mistake after EMD decomposes! Do not find Reference source.(ignoring residual error) can be expressed as

In embodiments of the present invention, step 3) the middle various quantizating index that calculate different frequency intrinsic mode function, as the specific implementation step of energy, Hilbert weighted frequency, locally coherence etc. is:

Step 3.1, calculates the energy of each IMF;

Step 3.2, calculates the Hilbert weighted frequency (HWF) of each IMF and the Hilbert weighted average frequency (HWMF) of primary signal;

Step 3.3, calculates other indexs, as locally coherence (ReHo) etc.

In embodiments of the present invention, the circular that calculates the energy of each IMF is:

In embodiments of the present invention, the specific implementation step of calculating the Hilbert weighted frequency (HWF) of each IMF and the Hilbert weighted average frequency (HWMF) of primary signal is:

Step 3.2.1, establishes primary signal x (t) length for a time point, and EMD decomposition obtains K IMF, and the HWF of j IMF is

wherein a _jand w (i) _j(i) be respectively envelope and the instantaneous frequency of the moment i in step 2.2.2;

Step 3.2.2, the HWMF of primary signal x (t) is by { the HWF of all IMF _j, j=1 ..., K, obtains through the weighting of each IMF envelope norm,

In the present embodiment, calculate other indexs, as the specific implementation step of locally coherence (ReHo) etc. is:

The all voxels of full brain are calculated to its certain IMF and the ReHo of the identical IMF of adjacent voxels around one by one,

wherein W is exactly ReHo value mistake! Do not find Reference source., scope 0 to 1,

wherein r _ijj voxel mistake! Do not find Reference source.In the order of i time point; r _iaverage, n mistake! Do not find Reference source.Be the time point number of this IMF, total number of this voxel of k and adjacent voxels, if adjacent voxels is defined as the k=27 of voxel contacting with this voxel face, limit or fixed point, if be defined as face or the tactile k=19 of voxel of edge joint, if be defined as the k=7 of voxel of face contact.

In the present embodiment, step 4) in take above-mentioned quantizating index as classification foundation, by different voxel clusters in brain, and the specific implementation step of the above different indexs of statistical analysis is:

Step 4.1, the indices obtaining using step 3 or wherein some as feature, carries out K-means cluster to all voxels of full brain, also can use other clustering algorithms.

Step 4.2, utilizes the mind map of the indices that SPM obtains step 3 to carry out statistical test, obtains final result not carrying out on the same group statistical between object of study.

In the present embodiment, the indices obtaining using step 3 or wherein some as feature, the concrete steps of all voxels of full brain being carried out to K-means cluster are:

Step 4.1.1 chooses at random the characteristic index barycenter of K voxel from all voxels of full brain;

Step 4.1.2, measures the distance of its characteristic index to each barycenter to remaining each voxel, and it is grouped into the class of nearest barycenter, makes each voxel be grouped into a certain apoplexy due to endogenous wind until travel through all voxels of full brain;

Step 4.1.3, recalculates the barycenter (average of such characteristic index) of each class having obtained;

Step 4.1.4, iteration 4.1.2～4.1.3 walks until new barycenter equates or be less than assign thresholds with former barycenter, and algorithm finishes.

Below in conjunction with drawings and the specific embodiments, application principle of the present invention is further described.

The embodiment of the present invention with analyze normal tested under tranquillization state the image of magnetic resonance imaging gained be treated to example, the FCA method of attempting to propose by the present invention drives BOLD-fMRI signal data the composition that is decomposed into different frequency range, and then the various features of these frequency contents of quantitative analysis, by the voxel cluster of feature similarity, between research Different brain region, which frequency range to carry out information interaction and integration by.

The technical thought that realizes this object is: normal tested tranquillization state data are first carried out to pretreatment; Adopt again Hilbert-Huang transform (HHT), by pretreated BOLD signal adaptive be decomposed into the signal of different frequency range; Calculate subsequently that energy that HHT decomposes the each frequency band signals obtaining, Hilbert weighted frequency (HWF), locally coherence (ReHo) etc. are all kinds of to be estimated or quantizating index; Finally using these estimate or quantizating index as classification foundation by the voxel of different parts in brain (voxels) cluster.。

Concrete technical scheme comprises as follows:

(1) data pre-treatment step:

(1a) the SPM software based on MATLAB (Matrix Lab) platform, carries out time difference rectification to the tranquillization state data that collect;

(1b) time is corrected to data later and carry out the moving alignment of head;

(1c) data based on after the moving alignment of head, carry out Spatial normalization;

(1d) data after Spatial normalization are gone linear drift and time domain standardization;

(2) seasonal effect in time series step after HHT decomposition pretreatment:

(2a) adopt EMD to decompose signal, obtain some IMF;

(2b) IMF obtaining after EMD decomposition is carried out to Hilbert conversion;

(3) step of the various quantizating index of calculating different frequency composition (IMF):

(3a) calculate the energy of each IMF;

(3b) calculate the Hilbert weighted frequency (HWF) of each IMF and the Hilbert weighted average frequency (HWMF) of primary signal;

(3c) profit is calculated other indexs, as locally coherence (ReHo) etc.;

(4) take above-mentioned quantizating index as foundation is by different voxel clusters in brain, and the specific implementation step of the above different indexs of statistical analysis is:

(4a) indices obtaining using step 3 or wherein some as feature, carries out K-means cluster (also can use other clustering algorithms) to all voxels of full brain;

(4b) utilize the mind map of the indices that SPM obtains step 3, obtain final result not carrying out on the same group statistical between object of study.

With reference to Fig. 1, the present invention includes: to the data pretreatment collecting from magnetic resonance machine; Adopt again HHT by pretreated BOLD signal adaptive be decomposed into the signal of different frequency range; Calculating subsequently energy that HHT decomposes the each frequency band signals obtaining, Hilbert weighted frequency, locally coherence etc. estimates or quantizating index; Finally using these estimate or quantizating index as classification foundation by the voxel cluster of different parts in brain.Concrete steps are as follows:

One, data preprocessing phase

Step 1: to magnetic resonance acquisition to data carry out time difference rectification, it is exactly the difference of correcting in 1 volume acquisition time between layers that time difference is corrected, and then guarantees between each layer it is all to obtain from the identical time;

Step 2: because cerebral function imaging duration of experiment is long, the head movement that the tested physiologic factor such as breathing, blood flow pulsation causes is unavoidable, so the data that the time corrected are carried out to the moving alignment of head, namely each two field picture of a sequence is all carried out to registration with the first two field picture of this sequence, be registrated under the same coordinate system, moving to correct head, after alignment, the data after correct moving rectification are carried out hand inspection again, spend and just exclude, will not analyze if translation and rotation exceed respectively 1 millimeter and 1.

Step 3: the data space standardization after the moving alignment of head, in test, exist multiple tested, there is certain difference in the brain size shape between tested and tested, for follow-up statistical analysis, must carry out the normalization of brain shape, tested brain registration be normalized in the brain template of standard;

3a) data after correct moving alignment utilize the affine transformation the average image of 12 parameters to be registrated to MNI standard form;

3b) be heavily cut to the voxel of 3x3x3 cubic millimeter, and MNI (Montreal Neurological Institute) coordinate transform is arrived to Talairach coordinate system;

Step 4: by linear regression, the composition of linear growth or reduction is removed, the time series of each voxel is deducted to the time domain average of self the time domain standard deviation divided by self, to eliminate the difference of different voxel baseline values.

Two, the HHT stage

Step 5: with reference to Fig. 2, this stage is utilized the program based on MATLAB platform of oneself writing, adopts EMD to decompose pretreated signal, and flow process is as follows:

5a) find all extreme points of primary signal x (t);

5b) between all maximum points, interpolation obtains the coenvelope eup (t) of x (t), and between all minimum points, interpolation obtains the lower envelope elow (t) of x (t);

5c) average time series emean (t)=[elow (t)+eup (t)]/2 of calculating eup (t) and elow (t);

5d) x (t) deducts average time series and obtains " concussion pattern " r (t)=x (t)-emean (t);

Whether 5e) judgement " disconnected concussion pattern " r (t) meets feature or the monotonous curve of IMF, should " disconnected concussion pattern " become an IMF, IMF if meet IMF feature _i(t)=r (t), if monotonous curve no r (t) become residual error, EMD process finishes, and repeats above step until obtain an IMF otherwise establish x (t)=r (t).IMF is characterized as lower envelope symmetry, and the number of times of IMF curve zero crossing equates with extreme point number or differs at most 1.

Step 6: EMD is decomposed to the IMF obtaining and carry out Hilbert transform:

6a) establishing certain IMF is X (t) mistake! Do not find Reference source.Its Hilbert conversion Y (t) mistake! Do not find Reference source.For mistake! Do not find Reference source.

wherein P is Cauchy's principal value $P{\∫}_a^{b}f\left(x\right)\mathrm{dx}={\lim }_{e\→\∞}({\∫}_a^{x-e}f\left(x\right)\mathrm{dx}+{\∫}_{x+e}^{b}f\left(x\right)\mathrm{dx})$ A mistake! Do not find Reference source.；

6b) X (t) and Y (t) are merged into a complex signal Z (t), order ^{a mistake! Do not find Reference source.}z (t)=X (t)+iY (t)=a (t) ^{ei θ (t)}, wherein a (t)=[X2 (t)+Y2 (t)] ^1/2the envelope of complex signal, Q (t)=tan ^-1(Y/X) be the instantaneous phase of complex signal.The instantaneous frequency of hilbert spectrum is defined as

6c) by all IMF are done to Hilbert transform, primary signal x (t) mistake after EMD decomposes! Do not find Reference source.(ignoring residual error) can be expressed as

Three, calculate the various quantizating index stages of each IMF

Step 7: calculate each IMF energy:

Step 8: calculate the Hilbert weighted frequency (HWF) of each IMF and the Hilbert weighted average frequency (HWMF) of primary signal,

8a) establish primary signal x (t) length for a time point, EMD decomposition obtains K IMF, and the HWF of j IMF is wherein a _jand w (i) _j(i) be respectively envelope and the instantaneous frequency of the moment i in step 2.2.2.

8b) HWMF of primary signal x (t) is by { the HWF of all IMF _j, j=1 ..., j, obtains through the weighting of each IMF envelope norm, $\mathrm{HWMF}=\frac{{\mathrm{\Σ}}_{i=1}^N\left|\right|a_j\left|\right|{\mathrm{HWF}}_j}{{\mathrm{\Σ}}_{i=1}^N\left|\right|a_j\left|\right|}.$

Step 9: calculate certain voxel and the locally coherence of the identical IMF of voxel around, all voxels of full brain are calculated to its certain IMF and the ReHo of the identical IMF of adjacent voxels around one by one,

wherein W is exactly ReHo value mistake! Do not find Reference source., scope 0 to 1,

wherein r _ijj voxel mistake! Do not find Reference source.In the order of i time point;

r _iaverage, n mistake! Do not find Reference source.Be the time point number of this IMF, total number of this voxel of k and adjacent voxels, if adjacent voxels is defined as the k=27 of voxel contacting with this voxel face, limit or fixed point, if be defined as face or the tactile k=19 of voxel of edge joint, if be defined as the k=7 of voxel of face contact.

Four, cluster and the statistics stage

Step 10: the indices obtaining using step 3 or wherein some as feature, all voxels of full brain are carried out to K-means cluster, also can use other clustering algorithms;

10a) from all voxels of full brain, choose at random the characteristic index barycenter of K voxel;

10b) measure the distance of its characteristic index to each barycenter to remaining each voxel, and it is grouped into the class of nearest barycenter, make each voxel be grouped into a certain apoplexy due to endogenous wind until travel through all voxels of full brain;

10c) recalculate the barycenter (average of such characteristic index) of each class having obtained;

10d) iteration 4.1.2～4.1.3 walks until new barycenter equates or be less than assign thresholds with former barycenter, and algorithm finishes.

Step 11: utilize the mind map of the indices that SPM obtains step 3 to add up, obtain final result by two sample t-checks not carrying out on the same group statistical between object of study.

By above 11 steps, realize without priori BOLD signal decomposition and the rounded analysis to BOLD signal different frequency characteristic.

Fig. 3 is the result of carrying out cluster according to voxel frequency characteristic, different classes represents by different colors, wherein the spatial distribution of some class is consistent with the functional network that research had been found that in the past, the English abbreviation of these representational classes marks in the drawings, and the average HWMF value of all voxels of these representative classes also marks in the drawings.Through full brain BOLD signal adaptive is decomposed and quantitative analysis, find the Different brain region frequency that takes on a different character, different Nao district or functional node are integrated together and are formed large-scale brain network by specific frequency.Above-mentionedly be found to be neuroscience basic research Brain function integration and separation, parallel task processing provides new approaches, also for changing, the central nervous system of clinical research various diseases provides reference baseline and iconography new method, obtain good result of study, to prove that the method is feasible.

Compare with existing research, the present invention has advantages of as follows:

1 research in the past, by artificial setpoint frequency scope research BOLD frequency characteristic, is used band filter to obtain the composition of different frequency range simultaneously, the randomness that this causes frequency range to be chosen, and do not consider the nonlinear feature of BOLD signal non-stationary.The present invention with the method for data-driven without artificially setting any parameter and the condition including frequency range, simultaneously adaptively by BOLD signal decomposition to different frequency range, considered the nonlinear feature of BOLD signal non-stationary.

2 only analyze that original BOLD data obtain the Global Information of full frequency band and the characteristic information that cannot obtain each frequency range compared to research in the past, the present invention is by HHT and multiplely estimate and quantizating index (energy, HWF, ReHo) combination, from many aspects quota portray the feature of each voxel different frequency BOLD activity (different I MF), also can use and estimate and quantizating index simultaneously, the characteristic of portraying different I MF as entropy, small-world network etc., thus the present invention can excavate more useful informations of BOLD signal hiding.

The frequency clustering method that the embodiment of the present invention provides, first to normal tested tranquillization state data are carried out to pretreatment, adopt again Hilbert-Huang transform, by pretreated BOLD signal adaptive be decomposed into the signal of different frequency range, calculate subsequently that energy that HHT decomposes the each frequency band signals obtaining, Hilbert weighted frequency, locally coherence etc. are all kinds of to be estimated or quantizating index, finally using these estimate or quantizating index as classification foundation by the voxel cluster of different parts in brain.Found the Different brain region frequency that takes on a different character, different Nao district or functional node are integrated together and are formed large-scale brain network by specific frequency.These are found to be neuroscience basic research Brain function integration and separation, parallel task processing provides new approaches, also provide reference baseline and iconography new method for the central nervous system of clinical research various diseases changes, and have obtained good result of study.

BOLD-fMRI frequency provided by the invention clustering method, for find various nervus centralis diseases focus, drug target, probe into sacred disease pathogenesis new method be provided, the support of methodology, iconography is also provided for the basic research of cognitive neuroscience, neuro physiology, has there is important scientific research and clinical value.

The foregoing is only preferred embodiment of the present invention, not in order to limit the present invention, all any modifications of doing within the spirit and principles in the present invention, be equal to and replace and improvement etc., within all should being included in protection scope of the present invention.

Claims (10)

1. blood sample horizontal dependency functional magnetic resonance signal concussion frequency clustering method, is characterized in that: comprise the following steps:

Step 1: the tranquillization state data that collect from magnetic resonance machine are carried out to Preprocessing;

Step 2: adopt Hilbert-Huang transform to decompose pretreated signal, obtain the intrinsic mode function of different frequency; Hilbert-Huang transform comprises empirical mode decomposition and Hilbert transform two steps;

Step 3: the quantizating index that calculates different frequency composition;

Step 4: take the quantizating index of step 3 as classification foundation, by different voxel clusters in brain, the above different indexs of statistical analysis.

2. blood sample horizontal dependency functional magnetic resonance signal concussion frequency according to claim 1 clustering method, it is characterized in that: in described step 1, be BOLD-fMRI signal data to the tranquillization state data that collect from magnetic resonance machine, the specific implementation step of BOLD-fMRI signal data being carried out to Preprocessing is:

1), utilize SPM software based on MATLAB platform, the tranquillization state data that collect are carried out to time difference rectification;

2), the time is corrected to data later and carry out the moving alignment of head;

3), the data after correct moving alignment are carried out Spatial normalization;

4), the data after Spatial normalization are gone to linear drift and time domain standardization.

3. blood sample horizontal dependency functional magnetic resonance signal concussion frequency according to claim 1 clustering method, it is characterized in that: in described step 2, adopt Hilbert-Huang transform to decompose pretreated signal, obtain the composition of different frequency, be intrinsic mode function, specific implementation step is:

1), adopt empirical mode decomposition signal is decomposed, obtain some intrinsic mode functions;

2), the intrinsic mode function obtaining after empirical mode decomposition is carried out to Hilbert transform.

4. blood sample horizontal dependency functional magnetic resonance signal concussion frequency according to claim 3 clustering method, it is characterized in that: adopt Hilbert-Huang transform to decompose pretreated signal, wherein in step 1), adopt empirical mode decomposition to decompose signal, the specific implementation step that obtains some intrinsic mode functions is:

Step 1.1), find all extreme points of primary signal x (t);

Step 1.2), between all maximum points, interpolation obtains the coenvelope eup (t) of x (t), between all minimum points, interpolation obtains the lower envelope elow (t) of x (t);

Step 1.3), calculate eup (t) and average time series emean (t)=[elow (t)+eup (t)]/2 of elow (t);

Step 1.4), x (t) deducts average time series and obtains concussion pattern a: r (t)=x (t)-emean (t);

Step 1.5), whether the disconnected concussion pattern r (t) of judgement meet feature or the monotonous curve of intrinsic mode function, if meet IMF feature, should " disconnected concussion pattern " become an intrinsic mode function, IMF _i(t)=r (t); Monotonous curve no r (t) becomes residual error, and EMD process finishes, and repeats above step until obtain an IMF otherwise establish x (t)=r (t); IMF is characterized as lower envelope symmetry, and the number of times of IMF curve zero crossing equates with extreme point number or differs at most 1.

5. blood sample horizontal dependency functional magnetic resonance signal concussion frequency according to claim 3 clustering method, it is characterized in that, described employing Hilbert-Huang transform decomposes pretreated signal, wherein step 2) in the intrinsic mode function obtaining after empirical mode decomposition is carried out to Hilbert transform specific implementation step be:

Step 2.1), to establish certain IMF be that its Hilbert transform Y (t) of X (t) is

wherein P is Cauchy's principal value $P{\∫}_a^{b}f\left(x\right)\mathrm{dx}={\lim }_{e\→\∞}({\∫}_a^{x-e}f\left(x\right)\mathrm{dx}+{\∫}_{x+e}^{b}f\left(x\right)\mathrm{dx});$

Step 2.2), by X (t) and Y(t) merge into a complex signal Z(t), make Z (t)=X (t)+iY (t)=a (t) ^{ei θ (t)}, wherein a (t)=[X2 (t)+Y2 (t)] ^1/2the envelope of complex signal, Q(t)=tan ^-1(Y/X) be the instantaneous phase of complex signal; The instantaneous frequency of hilbert spectrum is defined as

Step 2.3), by all intrinsic mode functions are done to Hilbert transform, the primary signal x (t) after empirical mode decomposition can be expressed as $x\left(t\right)={\mathrm{\Σ}}_{i=1}^K{a}_i\left(t\right)\exp \{i\∫w_j\left(t\right)\mathrm{dt}\}.$

6. blood sample horizontal dependency functional magnetic resonance signal concussion frequency according to claim 1 clustering method, is characterized in that, in step 3, calculates different frequency composition, i.e. the quantizating index of intrinsic mode function, and specific implementation step is:

Step 3.1), calculate the energy of each intrinsic mode function;

Step 3.2), calculate the Hilbert weighted frequency of each intrinsic mode function and the Hilbert weighted average frequency of primary signal;

Step 3.3), calculate other indexs, as locally coherence.

7. blood sample horizontal dependency functional magnetic resonance signal according to claim 6 concussion frequency clustering method, is characterized in that: in described step 3.1, the circular that calculates the energy of each intrinsic mode function is,

8. blood sample horizontal dependency functional magnetic resonance signal concussion frequency according to claim 5 clustering method, it is characterized in that: described step 3.2) in, the specific implementation step of calculating the Hilbert weighted frequency of each intrinsic mode function and the Hilbert weighted average frequency of primary signal is:

Step 3.2.1), to establish primary signal x (t) length be N time point, empirical mode decomposition obtains K intrinsic mode function, the Hilbert weighted frequency of j intrinsic mode function is

wherein a _jand w (i) _j(i) envelope and the instantaneous frequency of the moment i in being respectively;

Step 3.2.2), the Hilbert weighted average frequency of primary signal x (t) by all intrinsic mode functions HWFj}, j=1 ..., K, obtains through the weighting of each intrinsic mode function envelope norm, $\mathrm{HWMF}=\frac{{\mathrm{\Σ}}_{i=1}^N\left|\right|a_j\left|\right|{\mathrm{HWF}}_j}{{\mathrm{\Σ}}_{i=1}^N\left|\right|a_j\left|\right|}.$

9. blood sample horizontal dependency functional magnetic resonance signal concussion frequency according to claim 6 clustering method, is characterized in that step 3.3) in, specific implementation step is:

The all voxels of full brain are calculated to its certain intrinsic mode function and the ReHo of the identical intrinsic mode function of adjacent voxels around one by one, wherein W is exactly ReHo value mistake! Do not find Reference source., scope 0 to 1, wherein r _ijj voxel mistake! Do not find Reference source.In the order of i time point;

r _iaverage, n mistake! Do not find Reference source.It is the time point number of this intrinsic mode function, total number of this voxel of k and adjacent voxels, if adjacent voxels is defined as the k=27 of voxel contacting with this voxel face, limit or fixed point, if be defined as face or the tactile k=19 of voxel of edge joint, if be defined as the k=7 of voxel of face contact.

10. blood sample horizontal dependency functional magnetic resonance signal concussion frequency according to claim 1 clustering method, it is characterized in that: in described step 4, take above-mentioned quantizating index as classification foundation, by different voxel clusters in brain, and the specific implementation step of the above different indexs of statistical analysis is:

Step 4.1), the indices obtaining using step 3 or wherein some as feature, carries out K-means cluster to all voxels of full brain;

Step 4.2), utilize the mind map of the indices that SPM software obtains step 3, obtain final result not carrying out on the same group statistical between object of study;

Wherein, step 4.1) in the indices that obtains using step 3 or wherein some as feature, the concrete steps of all voxels of full brain being carried out to K-means cluster are:

Step 4.1.1), from all voxels of full brain, choose at random the characteristic index barycenter of K voxel;

Step 4.1.2), measure the distance of its characteristic index to each barycenter to remaining each voxel, and it is grouped into the class of nearest barycenter, make each voxel be grouped into a certain apoplexy due to endogenous wind until travel through all voxels of full brain;

Step 4.1.3), recalculate the barycenter (average of such characteristic index) of each class having obtained;

Step 4.1.4), iteration 4.1.2～4.1.3 walks until new barycenter equates or be less than assign thresholds with former barycenter, and algorithm finishes.

Images