CN115097533B - 一种基于tls-esprit算法的磁共振测深信号提取方法 - Google Patents

一种基于tls-esprit算法的磁共振测深信号提取方法 Download PDF

Info

Publication number
CN115097533B
CN115097533B CN202210479823.1A CN202210479823A CN115097533B CN 115097533 B CN115097533 B CN 115097533B CN 202210479823 A CN202210479823 A CN 202210479823A CN 115097533 B CN115097533 B CN 115097533B
Authority
CN
China
Prior art keywords
matrix
frequency
mrs signal
signal
mrs
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202210479823.1A
Other languages
English (en)
Other versions
CN115097533A (zh
Inventor
蒋川东
吴昊
孙超
田宝凤
郝尚帅
刘峰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Jilin University
Original Assignee
Jilin University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Jilin University filed Critical Jilin University
Priority to CN202210479823.1A priority Critical patent/CN115097533B/zh
Publication of CN115097533A publication Critical patent/CN115097533A/zh
Application granted granted Critical
Publication of CN115097533B publication Critical patent/CN115097533B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/14—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation operating with electron or nuclear magnetic resonance
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N24/00—Investigating or analyzing materials by the use of nuclear magnetic resonance, electron paramagnetic resonance or other spin effects
    • G01N24/08—Investigating or analyzing materials by the use of nuclear magnetic resonance, electron paramagnetic resonance or other spin effects by using nuclear magnetic resonance
    • G01N24/082—Measurement of solid, liquid or gas content
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/38—Processing data, e.g. for analysis, for interpretation, for correction
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30—Assessment of water resources

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • High Energy & Nuclear Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • Remote Sensing (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Geophysics (AREA)
  • Chemical & Material Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • General Health & Medical Sciences (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)

Abstract

本发明涉及一种磁共振测深信号噪声滤除领域，具体是一种基于TLS‑ESPRIT算法的磁共振测深信号提取方法。将含噪信号进行希尔伯特变换，后求取自相关矩阵，进行特征分解，选取K+1个向量构建矩阵，矩阵分解和变换得到特征值参数，获得信号参数估计值，最后重构信号。该方法可以解决复杂环境噪声干扰下难以有效估计工频谐波基频，导致传统谐波建模方法无法干净去除工频谐波干扰的问题。

Description

一种基于TLS-ESPRIT算法的磁共振测深信号提取方法
技术领域
本发明涉及一种磁共振测深(Magnetic Resonance Sounding,MRS)信号噪声滤除领域，具体是一种基于TLS-ESPRIT算法的磁共振测深信号提取方法。
背景技术
磁共振测深(Magnetic Resonance Sounding,MRS)，是目前地下水探测的先进技术，具有便利、高效、快速且对地下水体无损伤探测的特点。通过对地下水分子中的氢质子进行激发，使氢质子产生能级跃迁，撤掉激发信号后，采集氢质子释能回迁时所感应出的MRS信号，从而获取地下水的相关信息。有效MRS信号是一个衰减震荡信号，其数学模型为：
Figure BDA0003627123840000011
Figure BDA0003627123840000012
该模型中的相关参数主要包括：E0、
Figure BDA0003627123840000013
fL、/>
Figure BDA0003627123840000014
分别表示与地层含水量相关的信号初始振幅、与含水层孔隙度相关的平均横向弛豫时间、与地理位置有关的拉莫尔频率、与地质结构和水层周围物质的导电性相关的初始相位。
由于MRS信号一般为纳伏级别的微弱信号，所以实际采集的初始信号会包含大量的噪声干扰，使得MRS信号被各种复杂噪声所影响，数据质量较差，难以实现有效MRS信号的高精度提取，进而影响数据反演解释结果的准确度。(放在这里可能会起反作用删除了)
发明内容
针对MRS信号中所包含噪声的多样性，以及传统谐波建模算法在谐波基频不稳定或者含有多个基频时较难实现将工频谐波噪声去除干净，难以实现MRS信号有效提取问题，本发明提供了一种基于TLS-ESPRIT算法的磁共振测深信号提取方法，该方法通过对MRS信号中初始振幅、平均横向弛豫时间、拉莫尔频率以及初始相位4个关键参数的有效估计，通过数学模型重构，实现了对MRS信号的有效提取，且算法参数估计的精度高、计算量小、速度快。TLS-ESPRIT算法还适用于激光雷达和扫描雷达等特征参数提取领域。
本发明是这样实现的，
一种基于TLS-ESPRIT算法的磁共振测深信号提取方法，该方法包括如下步骤：
步骤1：对磁共振测深过程中的所采集到的观测磁共振信号，进行希尔伯特变换，得到其复数域值S(n)，并对S(n)按照阵元数M，快拍数L，以d为间隔构建数据矩阵X(n)，求取X(n)的自相关矩阵RSS，自相关矩阵RSS的大小为M×M，L+(M-1)d小于等于信号S(n)的数据长度；其中：
Figure BDA0003627123840000021
步骤2：对自相关矩阵RSS进行特征分解，对得到的特征值按照由大到小的顺序进行排列，与特征值相对应的特征向量也进行对应的排列；
步骤3：对采集到的观测磁共振信号进行傅里叶变换，得到其频谱，确定信号中的工频谐波干扰频点f1、f2......fK，进而确定工频谐波干扰数量K，基于旋转不变子空间算法，根据步骤2中的特征值排列，选取前K+1个对应的特征向量，得到Ex、Ey和Exy，Ex对应1～M-1个阵元的K+1个特征值对应的特征向量，Ey对应2～M个阵元的K+1个特征值对应的特征向量，Exy＝[Ex Ey]；
步骤4：基于总体最小二乘准则，构建数据矩阵Es＝Exy H×Exy，H表示共轭转置,并对数据矩阵Es进行特征分解和矩阵变换，求得特征值λi，i＝1,2,3,…,K+1，计算出各个谐波的频率和MRS信号的拉莫尔频率；
Figure BDA0003627123840000031
其中ang(·)表示取复数的辐角，fs为采样频率，fi为各个谐波的频率和MRS信号的拉莫尔频率估计值fL；
步骤5：将式(2)中最接近于已知发射频率的fi确定为MRS信号的拉莫尔频率fL，并依据fL对应的特征值λL，计算出MRS信号的平均横向弛豫时间T2 *；
Figure BDA0003627123840000032
步骤6：根据步骤4和步骤5获得的工频谐波的频率fi以及MRS信号的拉莫尔频率fL、平均横向弛豫时间T2 *的数值，构造函数矩阵Q，其中函数矩阵Q的前K列是由工频谐波成分构成的Q(l,k)＝exp(j2πfi(k)l/fs),l＝0,1,…,N+K；k＝1,2,…,K；Q的第K+1列是由MRS信号成分构成的
Figure BDA0003627123840000033
Figure BDA0003627123840000034
k＝K+1；l＝0,1,…,N+K；根据采集的数据得到数值矩阵Y＝[S(n)′],n＝1,2,…,N+K+1，′表示转置，构成方程：
S×Q＝Y (4)
其中N为整数，取值小于快拍数L，对H矩阵求逆，得到矩阵S；对矩阵S取幅值abs(S)，得MRS信号的初始幅值E0；对矩阵S取幅角ang(S)，得MRS信号的初始相位
Figure BDA0003627123840000035
步骤7：根据上述获得的MRS信号的初始振幅E0，平均横向弛豫时间T2 *，拉莫尔频率fL和初始相位
Figure BDA0003627123840000036
按照其数学模型/>
Figure BDA0003627123840000037
Figure BDA0003627123840000038
进行MRS信号的重构，实现MRS信号的提取目标。
进一步地，步骤4中对矩阵Es进行特征分解和矩阵变换得到特征值，具体步骤为：
步骤4a：对矩阵Es进行特征分解，得到相应的特征矩阵EV；
步骤4b：将特征矩阵EV分成4个子空间矩阵Ei，计算
Figure BDA0003627123840000041
Figure BDA0003627123840000042
步骤4c：对ψTLS进行特征值分解，得到特征值λi，i＝1,2,3,…,K+1，该特征值包含着各个工频谐波和MRS信号的特征信息。
本发明与现有技术相比，有益效果在于：
本发明方法整体框架是通过对含噪信号进行自相关矩阵的构建，采用旋转不变子空间ESPRIT算法得到特征矩阵，对特征矩阵进行分解能够获得包含关于MRS信号以及工频谐波参数信息的特征值，可以求得MRS信号的拉莫尔频率和平均横向弛豫时间。同时，基于工频谐波和理想MRS信号的数学模型，利用总体最小二乘TLS算法，进一步求得MRS信号的初始振幅和初始相位，实现有效MRS信号的4个关键参数的估计，最后通过数学模型实现对MRS信号的提取。该方法可以实现MRS信号关键参数的高精度估计，进而实现对信号的重构提取。该方法可以解决复杂环境噪声干扰下难以有效估计工频谐波基频，导致传统谐波建模方法无法干净去除工频谐波干扰的问题。该方法的有效应用，对于提高磁共振探测数据的信噪比，拓展磁共振勘探技术在复杂强噪声干扰下有效MRS信号的高精度和高稳定性提取具有重要的理论意义，对于拓展仪器系统的应用范围具有重要价值。
附图说明
图1为本发明基于TLS-ESPRIT算法的磁共振测深信号提取方法的流程框图；
图2为本发明TLS-ESPRIT算法的流程框图；
图3理想MRS信号及频谱，图3A为时间-幅值图，图3B为频率-幅值图；
图4含噪与提取的MRS信号及其频谱，图4A为时间-幅值图，图4B为频率-幅值图；
图5实测与提取的MRS信号及其频谱，图5A为处理前后MRS信号时域图，图5B为对应的频谱。
具体实施方式
为了使本发明的目的、技术方案以及优点更加清晰直观，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
如图1所示，基于TLS-ESPRIT算法的磁共振测深信号提取方法，包括以下步骤：
步骤1：对磁共振测深过程中的所采集到的观测MRS信号，进行希尔伯特变换，得到其复数域值S(n)，并对S(n)按照阵元数M，快拍数L,以d为间隔构建数据矩阵X(n),求取X(n)的自相关矩阵RSS，RSS的大小为M×M，L+(M-1)d小于等于信号S(n)的数据长度；
Figure BDA0003627123840000051
步骤2：对自相关矩阵RSS进行特征分解，对得到的特征值按照由大到小的顺序进行排列，与特征值相对应的特征向量也进行对应的排列；
步骤3：对采集到的观测MRS信号进行傅里叶变换，得到其频谱，确定信号中的工频谐波干扰频点f1、f2......fK，进而确定工频谐波干扰数量K，基于旋转不变子空间(ESPRIT)算法，如图2所示，根据步骤2中的特征值排列，选取前K+1个对应的特征向量，得到Ex、Ey和Exy，Ex对应1～M-1个阵元的K+1个特征值对应的特征向量，Ey对应2～M个阵元的K+1个特征值对应的特征向量，Exy＝[Ex Ey]；
步骤4：基于总体最小二乘准则(TLS)，如图2所示，构建数据矩阵Es＝Exy H×Exy，H表示共轭转置,并对Es进行特征分解和矩阵变换，求得特征值λi(i＝1,2,3,…,K+1)，可计算出各个谐波的频率和MRS信号的拉莫尔频率；
Figure BDA0003627123840000061
其中ang(·)表示取复数的辐角，fs为采样频率，fi为各个谐波的频率和MRS信号的拉莫尔频率估计值fL。
步骤5：将式(2)中最接近于已知发射频率的fi确定为MRS信号的拉莫尔频率fL，并依据fL对应的特征值λL，则可计算出MRS信号的平均横向弛豫时间T2 *；
Figure BDA0003627123840000062
步骤6：根据步骤4和步骤5获得的工频谐波的频率fi以及MRS的拉莫尔频率fL、平均横向弛豫时间T2 *的数值，构造函数矩阵Q，其中Q的前K列是由工频谐波成分构成的Q(l,k)＝exp(j2πfi(k)l/fs),l＝0,1,…,N+K；k＝1,2,…,K；Q的第K+1列是由MRS信号成分构成的
Figure BDA0003627123840000063
Figure BDA0003627123840000064
k＝K+1；l＝0,1,…,N+K；根据采集的数据可以得到数值矩阵Y＝[S(n)′],n＝1,2,…,N+K+1，′表示转置，构成方程：
S×Q＝Y (4)
其中N为整数，取值小于快拍数L。对H矩阵求逆，可以得到S矩阵。对矩阵S取幅值abs(S)，可得MRS信号的初始幅值E0；对矩阵S取幅角ang(S)，可得MRS信号的初始相位
Figure BDA0003627123840000067
步骤7：根据上述获得的MRS信号的初始振幅E0，平均横向弛豫时间T2 *，拉莫尔频率fL和初始相位
Figure BDA0003627123840000065
按照其数学模型/>
Figure BDA0003627123840000066
Figure BDA0003627123840000071
进行MRS信号的重构，实现MRS信号的提取目标。
其中，步骤4中对矩阵Es进行特征分解和矩阵变换得到特征值，具体步骤为：
步骤4a：对矩阵Es进行特征分解，得到相应的特征矩阵EV；
步骤4b：将特征矩阵EV分成4个子空间矩阵Ei，计算
Figure BDA0003627123840000072
Figure BDA0003627123840000073
步骤4c：对ψTLS进行特征值分解，得到特征值λi(i＝1,2,3,…,K+1),该特征值就包含着各个工频谐波和MRS信号的特征信息。
实施例1
本实施例是在MATLAB R2018a编程环境下开展的本发明方法的仿真实验。
基于TLS-ESPRIT算法的磁共振测深信号噪声滤除方法，参照图1，包括以下步骤：
步骤1：利用式
Figure BDA0003627123840000074
构造拉莫尔频率为2325Hz，幅值e0为100nV，弛豫时间/>
Figure BDA0003627123840000075
为0.2s的理想MRS信号，如图3所示。其中，图3A为时间-幅值图，图3B为频率-幅值图。添加随机噪声，设置其均值为0，方差为1，幅值为0-100nV随机变化的白噪声；再添加工频谐波噪声，其幅值Am为0-150nV之间的随机数，工频谐波的基频变化范围bf＝49.9-50.1Hz，取这个范围内的随机数，发射频率fT＝2323Hz，信号采样频率fs＝25kHz，采样时间长度t＝1s。如图4所示给出了含噪MRS信号y(n)的时域和频域图，图4A为时间-幅值图，图4B为频率-幅值图；对该含噪数据进行希尔伯特变换求得复数域值S(n)，并对S(n)按照阵元数M＝1000，快拍数L＝6000,以d＝5为间隔构建数据矩阵X(n)，求取数据矩阵X(n)的自相关矩阵RSS，自相关矩阵RSS的大小为M×M，L+(M-1)d等于10995小于信号S(n)的数据长度25000。
Figure BDA0003627123840000081
步骤2：对自相关矩阵RSS进行特征分解，对得到的特征值按照由大到小的顺序进行排列，与特征值相对应的特征向量也进行对应的排列；
步骤3：对仿真的含噪MRS信号y(n)进行傅里叶变换，得到其频谱，确定信号中的工频谐波干扰频点f1、f2......fK，进而确定工频谐波干扰数量K＝71，基于旋转不变子空间(ESPRIT)算法，如图2所示，根据步骤2中的特征值排列，选取前72个对应的特征向量，得到Ex、Ey和Exy，Ex对应1～999个阵元的72个特征值对应的特征向量，Ey对应2～1000个阵元的72个特征值对应的特征向量，Exy＝[Ex Ey]；
步骤4：基于总体最小二乘准则(TLS)，如图2所示，构建数据矩阵Es＝Exy H×Exy，H表示共轭转置,并对Es进行特征分解得到相应的特征矩阵EV，将特征矩阵EV分成4个子空间矩阵
Figure BDA0003627123840000082
计算/>
Figure BDA0003627123840000083
对ψTLS进行特征值分解，求得特征值λi(i＝1,2,3,…,72)，可计算出各个谐波的频率和MRS信号的拉莫尔频率；
Figure BDA0003627123840000084
其中ang(·)表示取复数的辐角，fs为采样频率，fi为各个谐波的频率和MRS信号的拉莫尔频率估计值fL。
步骤5：根据式(2)，将最接近于已知发射频率fT＝2323Hz的fi确定为MRS信号的拉莫尔频率fL＝2325.0Hz；依据fL对应的特征值λL，根据式(3)计算出MRS信号的平均横向弛豫时间T2 *＝201.5ms；
Figure BDA0003627123840000085
步骤6：根据步骤4和步骤5获得的工频谐波的频率fi以及MRS的拉莫尔频率fL、平均横向弛豫时间T2 *的数值，构造函数矩阵Q，其中Q的前71列是由工频谐波成分构成的Q(l,k)＝exp(j2πfi(k)l/fs),l＝0,1,…,3071；k＝1,2,…,71；Q的第72列是由MRS信号成分构成的
Figure BDA0003627123840000091
k＝72；l＝0,1,…,3071；根据采集的数据可以得到数值矩阵Y＝[S(n)′],n＝1,2,…,3072，′表示转置，构成方程：S×Q＝Y。对Q矩阵求逆，可以得到S矩阵。对矩阵S取幅值abs(S)，可得MRS信号的初始幅值E0＝100.33nV；对矩阵S取幅角arg(S)，可得MRS信号的初始相位/>
Figure BDA0003627123840000092
步骤7：根据上述获得的MRS信号的初始振幅E0，平均横向弛豫时间T2 *，拉莫尔频率fL和初始相位
Figure BDA0003627123840000093
按照其数学模型/>
Figure BDA0003627123840000094
Figure BDA0003627123840000095
进行MRS信号的重构，实现MRS信号的提取目标。
为了验证本发明方法的效果，将去噪后MRS信号进行了信噪比(SNR)估计。经计算，其SNR＝40.78dB，较处理前的SNR提高了68.39dB；此外，对MRS信号进行了包络提取和数据拟合，获得的关键参数误差均满足应用要求。
实施例2
本实施例是以长春市文化广场实地采集的MRS信号作为本发明方法的处理对象。
基于TLS-ESPRIT算法的磁共振测深信号噪声滤除方法，参照图1，包括以下步骤：
步骤1：对磁共振测深(MRS)探水仪采集到的一组观测MRS信号y(n)，如图5所示，图5A为处理前后MRS信号时域波形图，图5B为对应的频谱图，计算其信噪比为-21.29dB；信号采样频率fs＝25kHz，采样时间长度t＝1s。发射频率为2360Hz，平均横向弛豫时间为0.2s，幅值约为100nV的数据进行算法处理实验。对该含噪数据进行希尔伯特变换求得复数域值S(n)，并对S(n)按照阵元数M＝1000，快拍数L＝6000,以d＝5为间隔构建数据矩阵X(n)，求取X(n)的自相关矩阵RSS，RSS的大小为M×M，L+(M-1)d等于10995小于信号S(n)的数据长度25000。
Figure BDA0003627123840000101
步骤2：对自相关矩阵RSS进行特征分解，对得到的特征值按照由大到小的顺序进行排列，与特征值相对应的特征向量也进行对应的排列；
步骤3：对观测MRS信号y(n)进行傅里叶变换，得到其频谱，确定信号中的工频谐波干扰频点f1、f2......fK，进而确定工频谐波干扰数量K＝61，基于旋转不变子空间(ESPRIT)算法，如图2所示，根据步骤2中的特征值排列，选取前62个对应的特征向量，得到Ex、Ey和Exy，Ex对应1～999个阵元的62个特征值对应的特征向量，Ey对应2～1000个阵元的62个特征值对应的特征向量，Exy＝[Ex Ey]；
步骤4：基于总体最小二乘准则(TLS)，如图2所示，构建数据矩阵Es＝Exy H×Exy，H表示共轭转置，并对Es进行特征分解得到相应的特征矩阵EV，将特征矩阵EV分成4个子空间矩阵
Figure BDA0003627123840000102
计算/>
Figure BDA0003627123840000103
对ψTLS进行特征值分解，求得特征值λi(i＝1,2,3,…,72)，可计算出各个谐波的频率和MRS信号的拉莫尔频率；
Figure BDA0003627123840000104
其中ang(·)表示取复数的辐角，fs为采样频率，fi为各个谐波的频率和MRS信号的拉莫尔频率估计值fL。
步骤5：根据式(2)，将最接近于已知发射频率fT＝2360Hz的fi确定为MRS信号的拉莫尔频率fL＝2360.0Hz；依据fL对应的特征值λL，根据式(3)计算出MRS信号的平均横向弛豫时间T2 *＝204.8ms。
Figure BDA0003627123840000111
步骤6：根据步骤4和步骤5获得的工频谐波的频率fi以及MRS的拉莫尔频率fL、平均横向弛豫时间T2 *的数值，构造函数矩阵Q，其中函数矩阵Q的前61列是由工频谐波成分构成的Q(l,k)＝exp(j2πfi(k)l/fs),l＝0,1,…,3061；k＝1,2,…,61；Q的第62列是由MRS信号成分构成的
Figure BDA0003627123840000112
Figure BDA0003627123840000113
k＝62；l＝0,1,…,3061；根据采集的数据可以得到数值矩阵Y＝[S(n)′],n＝1,2,…,3062，′表示转置，构成方程：S×Q＝Y。对Q矩阵求逆，可以得到S矩阵。对矩阵S取幅值abs(S)，可得MRS信号的初始幅值E0＝103.70nV。对矩阵S取幅角arg(S)，可得MRS信号的初始相位/>
Figure BDA0003627123840000114
步骤7：根据上述获得的MRS信号的初始振幅E0，平均横向弛豫时间T2 *，拉莫尔频率fL和初始相位
Figure BDA0003627123840000115
按照其数学模型/>
Figure BDA0003627123840000116
Figure BDA0003627123840000117
进行MRS信号的重构，实现MRS信号的提取目标。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种基于TLS-ESPRIT算法的磁共振测深信号提取方法，其特征在于，该方法包括如下步骤：
步骤1：对磁共振测深过程中的所采集到的观测磁共振信号，进行希尔伯特变换，得到其复数域值S(n)，并对S(n)按照阵元数M，快拍数L，以d为间隔构建数据矩阵X(n)，求取X(n)的自相关矩阵RSS，自相关矩阵RSS的大小为M×M，L+(M-1)d小于等于信号S(n)的数据长度；其中：
Figure FDA0004183003940000011
步骤2：对自相关矩阵RSS进行特征分解，对得到的特征值按照由大到小的顺序进行排列，与特征值相对应的特征向量也进行对应的排列；
步骤3：对采集到的观测磁共振信号进行傅里叶变换，得到其频谱，确定信号中的工频谐波干扰频点f1、f2......fK，进而确定工频谐波干扰数量K，基于旋转不变子空间算法，根据步骤2中的特征值排列，选取前K+1个对应的特征向量，得到Ex、Ey和Exy，Ex对应1～M-1个阵元的K+1个特征值对应的特征向量，Ey对应2～M个阵元的K+1个特征值对应的特征向量，Exy＝[Ex Ey]；
步骤4：基于总体最小二乘准则，构建数据矩阵Es＝Exy H×Exy，H表示共轭转置，并对数据矩阵Es进行特征分解和矩阵变换，求得特征值λi，i＝1，2，3，…，K+1，计算出各个谐波的频率和MRS信号的拉莫尔频率；
Figure FDA0004183003940000012
其中ang(·)表示取复数的辐角，fs为采样频率，fi为各个谐波的频率和MRS信号的拉莫尔频率估计值fL；
步骤5：将式(2)中最接近于已知发射频率的fi确定为MRS信号的拉莫尔频率fL，并依据fL对应的特征值λL，计算出MRS信号的平均横向弛豫时间T2 *；
Figure FDA0004183003940000021
步骤6：根据步骤4和步骤5获得的工频谐波的频率fi以及MRS信号的拉莫尔频率fL、平均横向弛豫时间T2 *的数值，构造函数矩阵Q，其中函数矩阵Q的前K列是由工频谐波成分构成的Q(l，k)＝exp(j2πfi(k)l/fs)，l＝0，1，…，N+K；k＝1，2，…，K；Q的第K+1列是由MRS信号成分构成的
Figure FDA0004183003940000022
Figure FDA0004183003940000023
k＝K+1；l＝0，1，…，N+K；根据采集的数据得到数值矩阵Y＝[S(n)′]，n＝1，2，…，N+K+1，′表示转置，构成方程：
S×Q＝Y (4)
其中N为整数，取值小于快拍数L，对函数矩阵Q求逆，得到矩阵S；对矩阵S取幅值abs(S)，得MRS信号的初始幅值E0；对矩阵S取幅角ang(S)，得MRS信号的初始相位
Figure FDA0004183003940000024
步骤7：根据上述获得的MRS信号的初始振幅E0，平均横向弛豫时间T2 *，拉莫尔频率fL和初始相位
Figure FDA0004183003940000025
按照其数学模型/>
Figure FDA0004183003940000026
Figure FDA0004183003940000027
进行MRS信号的重构，实现MRS信号的提取目标。
2.按照权利要求1所述的方法，其特征在于，步骤4中对矩阵Es进行特征分解和矩阵变换得到特征值，具体步骤为：
步骤4a：对矩阵Es进行特征分解，得到相应的特征矩阵EV；
步骤4b：将特征矩阵EV分成4个子空间矩阵Ei，计算
Figure FDA0004183003940000028
Figure FDA0004183003940000029
步骤4c：对ψTLS进行特征值分解，得到特征值λj，j＝1，2，3，…，K+1，该特征值包含着各个工频谐波和MRS信号的特征信息。
CN202210479823.1A 2022-05-05 2022-05-05 一种基于tls-esprit算法的磁共振测深信号提取方法 Active CN115097533B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210479823.1A CN115097533B (zh) 2022-05-05 2022-05-05 一种基于tls-esprit算法的磁共振测深信号提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210479823.1A CN115097533B (zh) 2022-05-05 2022-05-05 一种基于tls-esprit算法的磁共振测深信号提取方法

Publications (2)

Publication Number Publication Date
CN115097533A CN115097533A (zh) 2022-09-23
CN115097533B true CN115097533B (zh) 2023-06-30

Family

ID=83287085

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210479823.1A Active CN115097533B (zh) 2022-05-05 2022-05-05 一种基于tls-esprit算法的磁共振测深信号提取方法

Country Status (1)

Country Link
CN (1) CN115097533B (zh)

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102053280A (zh) * 2010-11-10 2011-05-11 吉林大学 带有参考线圈的核磁共振地下水探测系统及探测方法

Family Cites Families (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6133734A (en) * 1997-12-30 2000-10-17 Schlumberger Technology Corporation Method and apparatus for evaluating an earth formation using nuclear magnetic resonance techiques
AU2003900418A0 (en) * 2003-01-30 2003-02-13 Qr Sciences Limited Improvements in Signal Processing For Detection Of NQR Signals
US7292189B2 (en) * 2004-09-10 2007-11-06 Worcester Polytechnic Institute Methods and apparatus for high resolution positioning
CN104216021B (zh) * 2014-09-09 2017-03-22 吉林大学 一种基于分步式发射的地下核磁共振探测方法
EP3580586A1 (en) * 2017-02-09 2019-12-18 Services Pétroliers Schlumberger Geophysical deep learning
CN107607998B (zh) * 2017-09-25 2019-02-26 吉林大学 一种核磁共振找水仪磁共振响应信号参数提取方法及系统
CN107783200B (zh) * 2017-11-21 2019-06-07 吉林大学 一种联合emd与tfpf算法的全波磁共振信号随机噪声消减方法
CN109541306A (zh) * 2018-12-06 2019-03-29 华北电力大学 一种基于tls-esprit的谐波间谐波检测方法
CN109828318B (zh) * 2019-01-25 2020-07-17 吉林大学 一种基于变分模态分解的磁共振测深信号噪声滤除方法
CN110989017A (zh) * 2019-12-10 2020-04-10 吉林大学 一种包含变化频率偏量的地面核磁共振反演方法
CN112882111B (zh) * 2021-01-18 2022-05-03 吉林大学 一种基于循环相关的磁共振响应信号参数提取方法及系统
CN114280679A (zh) * 2022-01-10 2022-04-05 吉林大学 一种地面核磁共振信号参数提取方法及系统

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102053280A (zh) * 2010-11-10 2011-05-11 吉林大学 带有参考线圈的核磁共振地下水探测系统及探测方法

Also Published As

Publication number Publication date
CN115097533A (zh) 2022-09-23

Similar Documents

Publication Publication Date Title
Astone et al. A method for detection of known sources of continuous gravitational wave signals in non-stationary data
CN109828318B (zh) 一种基于变分模态分解的磁共振测深信号噪声滤除方法
US5574639A (en) System and method for constructing filters for detecting signals whose frequency content varies with time
Perelli et al. Model-based compressive sensing for damage localization in lamb wave inspection
CN107607998B (zh) 一种核磁共振找水仪磁共振响应信号参数提取方法及系统
CN113887398A (zh) 一种基于变分模态分解和奇异谱分析的gpr信号去噪方法
CN114280679A (zh) 一种地面核磁共振信号参数提取方法及系统
CN115097533B (zh) 一种基于tls-esprit算法的磁共振测深信号提取方法
Wang et al. Simulation of matched field processing localization based on empirical mode decomposition and Karhunen-Loeve expansion in underwater waveguide environment
CN109507292A (zh) 一种信号提取方法
Hou et al. Weak Signal Detection Based on Lifting Wavelet Threshold Denoising and Multi-Layer Autocorrelation Method.
CN110244360A (zh) 基于有效频率波数域去混叠的地震数据分离方法及系统
DENG et al. S-transform spectrum decomposition technique in the application of the extraction of weak seismic signals
Liu et al. Numerical simulation of the Doppler spectrum of a flying target above dynamic oceanic surface by using the FEM-DDM method
Wu et al. The suppression of powerline noise for TEM with coded source based on independent component analysis
CN115826068A (zh) 一种基于自适应高斯滤波的mrs信号包络提取方法
CAO et al. A method of detecting seismic singularities using combined wavelet with fractal
CN109871784B (zh) 遗传算法优化匹配追踪的全波核磁共振信号噪声滤除方法
Rani et al. High-resolution spectral analysis methods for self-potential data inversion
CN112946741B (zh) 基于稀疏重建理论的方位各向异性弱信息提取方法
Jeng et al. A nonlinear method of removing harmonic noise in geophysical data
Dong et al. Research on sea clutter suppression using sparse dictionary learning
Chik et al. Comparing the performance of Fourier decomposition and wavelet decomposition for seismic signal analysis
CN112326017A (zh) 一种基于改进半经典信号分析的微弱信号检测方法
CN116432008A (zh) 基于局部Capon估计的SNMR信号参数估计方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant