CN113945981B - 一种浅海海底节点二次定位方法 - Google Patents

一种浅海海底节点二次定位方法 Download PDF

Info

Publication number
CN113945981B
CN113945981B CN202111210053.2A CN202111210053A CN113945981B CN 113945981 B CN113945981 B CN 113945981B CN 202111210053 A CN202111210053 A CN 202111210053A CN 113945981 B CN113945981 B CN 113945981B
Authority
CN
China
Prior art keywords
kth
wave
submarine
shot
arrival
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.)
Expired - Fee Related
Application number
CN202111210053.2A
Other languages
English (en)
Other versions
CN113945981A (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.)
Changan University
Original Assignee
Changan 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 Changan University filed Critical Changan University
Priority to CN202111210053.2A priority Critical patent/CN113945981B/zh
Publication of CN113945981A publication Critical patent/CN113945981A/zh
Application granted granted Critical
Publication of CN113945981B publication Critical patent/CN113945981B/zh
Expired - Fee Related 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
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/38—Seismology; Seismic or acoustic prospecting or detecting specially adapted for water-covered areas
    • G01V1/3817—Positioning of seismic devices
    • G01V1/3835—Positioning of seismic devices measuring position, e.g. by GPS or acoustically
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/38—Seismology; Seismic or acoustic prospecting or detecting specially adapted for water-covered areas
    • G01V1/3843—Deployment of seismic devices, e.g. of streamers
    • G01V1/3852—Deployment of seismic devices, e.g. of streamers to the seabed

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Acoustics & Sound (AREA)
  • Oceanography (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种浅海海底节点二次定位方法，该方法包括以下步骤：一、炮点和海底节点的布设及地震波数据的采集；二、初至波初至时间的拾取及分离；三、参与定位的有效炮点的选择；四、海底节点二次定位目标函数的获取；五、海底节点二次定位位置的获取。本发明方法步骤简单，设计合理，利用直达波初至和折射波初至联合进行二次定位，准确且快速获取海底节点的位置，能够解决浅海仅利用直达波初至进行二次定位时存在的初至数量少、方位分布不均且受炮端干扰的问题，实用性强。

Description

一种浅海海底节点二次定位方法
技术领域
本发明属于海底节点二次定位技术领域，尤其是涉及一种浅海海底节点二次定位方法。
背景技术
目前，海底节点(OBN—Ocean Bottom Nodes)在海洋地震勘探中显示出极大优势，不仅能获取宽方位、大偏移距、含有转换横波的高质量地震数据，在施工中也不受调查船设备配置和电缆连接的限制，尤其灵活适用于船只和作业平台密集的海域。OBN在实际作业中根据GPS及其他定位系统按照预设测线位置布设，称为一次定位。在OBN布设过程中，由于受到人工抛掷误差、船速变化、沉降速度、洋流潮汐引起的漂移等影响，实际投放位置往往偏离设计位置，即使布设位置准确，在勘探工期内，由于渔船、气候变化、海洋生物等也会影响海底节点位置。忽略位置的变化，仍使用一次定位的位置处理地震波数据将会直接影响CMP面元的准确划分、动校正的正确使用等，对地震波数据的成像造成较大误差，所以要通过二次定位获取海底节点的真实位置。
目前主要的二次定位方法有声波定位和地震波定位。声波定位是使用声学设备及系统，能快速得到定位结果，效率较高，但设备昂贵、定位成本较高。地震波定位手段灵活，能满足施工的精度要求，同时能实现炮检双向定位，比较典型的方法有：基于波场的波场模拟法、波场延拓法，基于直达波初至的近正四面体法、搜索法、曲面扫描拟合法等。基于波场的定位方法需要从波场出发，做模拟或外推延拓，该方法计算量较大、计算速度慢、定位效率低。初至波定位大多是基于直达波初至，直达波定位广泛应用于深水域，在浅海地区存在许多问题。浅海域海底坚硬、地形平坦时容易形成稳定的折射面，折射波较为发育，地震波数据记录中的初至主要是折射波，直达波初至数量少、方位分布不均且因炮端的干扰初至难以准确拾取。仍使用直达波初至定位会降低定位的精度，甚至定位失败，因此，可以引入宽方位、远偏移距、数据量丰富的折射波初至参与定位。
发明内容
本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种浅海海底节点二次定位方法，其方法步骤简单，设计合理，利用直达波初至和折射波初至联合进行二次定位，准确且快速获取海底节点的位置，能够解决浅海仅利用直达波初至进行二次定位时存在的初至数量少、方位分布不均且受炮端干扰的问题，实用性强。
为解决上述技术问题，本发明采用的技术方案是：一种浅海海底节点二次定位方法，其特征在于，该方法包括以下步骤：
步骤一、炮点和海底节点的布设及地震波数据的采集：
步骤101、在海面上布设多条炮线，每条炮线上设置多个炮点；其中，炮线的数量为L，L条炮线按照震源船放炮的顺序依次标记为第1条炮线，…，第l条炮线，...，第L条炮线，l和L均为正整数，1≤l≤L，每条炮线上炮点的数量为N个，N个炮点位于同一直线上，N个炮点的编号按照震源船的运行方向由前至后依次标记为第1个炮点，...，第n个炮点，...，第N个炮点，n为炮点编号，n和N为正整数，且1≤n≤N；
步骤102、在海底设计多个海底检波点，并在每个海底检波点处投放海底节点；其中，海底检波点和海底节点的数量均为K，K个海底检波点位于同一直线上，K个海底检波点按照震源船的运行方向由前至后依次标记为第1个海底检波点，...，第k个海底检波点，...，第K个海底检波点，相应的K个海底节点依次标记为第1个海底节点，...，第k个海底节点，...，第K个海底节点，k和K均为正整数，且1≤k≤K；
步骤103、建立O-XYZ坐标系；其中，原点为O，X轴和Y轴形成的平面XOY处于海平面上，Z轴垂直平面XOY且向下为正；
步骤104、在O-XYZ坐标系下，设定第k个海底检波点的设计位置P(Xk,Yk,Zk)和第l条炮线中第n个炮点的坐标位置Q(xl,n,yl,n,zl,n)；其中，Xk,Yk,Zk分别表示第k个海底检波点在X轴方向、Y轴方向、Z轴方向的坐标分量；xl,n,yl,n,zl,n分别表示第l条炮线中第n个炮点在X轴方向、Y轴方向、Z轴方向的坐标分量；
步骤105、采用震源船依次在L条炮线上N个炮点激发地震波至海底的过程中，K个海底节点分别采集到L条炮线上N个炮点激发的地震波数据；
步骤二、初至波初至时间的拾取及分离：
步骤201、根据K个海底节点分别采集L条炮线上N个炮点激发的地震波数据，拾取第l个炮线上第n个炮点激发地震波至海底时第k个海底节点采集到的初至波初至时间Tk(l,n)；
步骤202、根据公式
Figure BDA0003308593480000031
得到出现折射波初至的临界距离d；其中，V1c表示地震波在海水中的传播速度测量值，V2c表示地震波在海底介质中的传播速度测量值，h表示海水深度；
步骤203、根据第k个海底节点采集到的初至波初至时间Tk(l,n)以及出现折射波初至的临界距离d，从步骤105中的地震波数据中获取第k个海底节点采集到的初至波为直达波所对应的多个炮点和第k个海底节点采集的初至波为折射波时所对应的多个炮点，并将第k个海底节点采集到的初至波为直达波所对应的多个炮点记作初至波为直达波集合，第k个海底节点采集的初至波为折射波时所对应的多个炮点记作初至波为折射波集合；
步骤三、参与定位的有效炮点的选择：
步骤301、将每一条炮线上相邻两个炮点之间的间距记作Ld；
步骤302、获取第k个海底节点的待选择炮点；
步骤303、剔除异常炮点获取第k个海底节点的有效炮点；
步骤304、根据步骤203中初至波为直达波集合和初至波为折射波集合，从步骤303中第k个海底节点的有效炮点中得到第k个海底节点采集的初至波为直达波的有效炮点和第k个海底节点采集的初至波为折射波的有效炮点；
步骤305、第k个海底节点采集的初至波为直达波的有效炮点按照炮线顺序和炮点顺序排序并依次标记为第1个直达波有效炮点，...，第p个直达波有效炮点，...，第P个直达波有效炮点；其中，p和P均为正整数，且1≤p≤P；
将第k个海底节点采集的初至波为折射波的有效炮点按照炮线顺序和炮点顺序排序并依次标记为第1个折射波有效炮点，...，第p′个折射波有效炮点，...，第P′个折射波有效炮点；其中，p′和P′均为正整数，且1≤p′≤P′；
步骤四、海底节点二次定位目标函数的获取：
步骤401、设定第k个海底节点的实际坐标为P′(Xk′,Yk′,Zk′)，根据公式
Figure BDA0003308593480000041
得到第p个直达波有效炮点的直达波到达第k个海底节点的理论初至时间td(p,k)；其中，xd(p),yd(p),zd(p)分别表示第p个直达波有效炮点在X轴方向、Y轴方向、Z轴方向的坐标分量；
步骤402、根据公式
Figure BDA0003308593480000042
得到P个直达波有效炮点的直达波到达第k个海底检波点的理论和实际初至时间的误差Qd(k)；其中，Td(p,k)表示第k个海底节点采集到第p个直达波有效炮点的初至波的实际初至时间；
步骤403、根据公式
Figure BDA0003308593480000043
得到第p′个折射波有效炮点的折射波到达第k个海底节点的理论初至时间tz(p′,k)；其中，xz(p′),yz(p′),zz(p′)分别表示第p′个折射波有效炮点在X轴方向、Y轴方向、Z轴方向的坐标分量，V1表示地震波在海水中的传播速度，V2表示地震波在海底介质中的传播速度；
步骤404、根据公式
Figure BDA0003308593480000051
得到P′个折射波有效炮点的折射波到达第k个海底检波点的理论和实际初至时间的误差Qz(k)；其中，Tz(p′,k)表示第k个海底节点采集到第p′个折射波有效炮点的初至波的实际初至时间；
步骤405、根据公式Q(k)＝w1×Qd(k)+w2×Qz(k)，得到第k个海底节点采集到的初至波的理论和实际初至时间的目标函数Q(k)；w1表示直达波加权因子，w2表示折射波加权因子；
步骤五、海底节点二次定位位置的获取：
采用模拟退火算法，得到目标函数Q(k)的实际坐标的最优解，并将实际坐标的最优解中P*(Xk*,Yk*,Zk*)记作第k个海底节点二次定位位置。
上述的一种浅海海底节点二次定位方法，其特征在于：步骤302中获取第k个海底节点的待选择炮点，具体过程如下：
步骤3021、以第k个海底节点所对应的第k个海底检波点的设计位置在XOY平面的投影为圆心，以半径[d+10Ld,d+30Ld]画圆，得到第k个海底节点待选择的炮点区域，并将第k个海底节点待选择的炮点区域记作第k个待选择炮点区域；
步骤3022、将第k个待选择炮点区域中的各个炮点记作第k个海底节点的待选择炮点；
步骤3023、将第k个海底节点的待选择炮点按照炮线顺序和炮点顺序排序分别记作第k1个待选择炮点，...，第ka个待选择炮点，...，第kA个待选择炮点；其中，a和A均为正整数，且1≤a≤A，A表示第k个海底节点的待选择炮点的总数，k1，ka和kA均为正整数。
上述的一种浅海海底节点二次定位方法，其特征在于：步骤303中剔除异常炮点获取第k个海底节点的有效炮点，具体过程如下：
步骤3031、将A个待选择炮点和第k个海底检波点的水平间距分别记作A个炮检距；根据步骤201中的初至时间得到A个待选择炮点分别激发地震波至海底时第k个海底节点采集到的初至波初至时间；
步骤3032、以炮检距为x轴，以初至时间为y轴，将步骤3031中的A个炮检距和初至时间进行绘制并拟合得到关系曲线；
步骤3033、将偏移关系曲线且偏移的初至时间大于100ms所对应的待选择炮点剔除，得到第k个海底节点的有效炮点。
上述的一种浅海海底节点二次定位方法，其特征在于：步骤五中采用模拟退火算法，得到目标函数Q(k)的实际坐标的最优解，具体过程如下：
步骤501、设定第k个海底节点的迭代的初始位置坐标为第k个海底节点所对应的第k个海底检波点的设计位置并记作P0(Xk,Yk,Zk)，地震波在海水中的传播速度V1的初始值V10为V1c，地震波在海底介质中的传播速度V2的初始值V20为V2c；
步骤502、设定直达波加权因子w1的初始值为w1(0)，折射波加权因子w2的初始值为w2(0)，并将w1(0)、w2(0)、P0(Xk,Yk,Zk)、V1c和V2c代入步骤405，得到初始目标函数值Q(k0)；其中，w1(0)和w2(0)均是(0，1)的随机数，且w1(0)+w2(0)＝1；
步骤503、设定第k个海底节点的迭代上限位置坐标为P(XkH,YkH,ZkH)，第k个海底节点的迭代下限位置坐标为P(Xkx,Ykx,Zkx)，V1的上限为V1s，V2的上限为V2s，V1的下限为V1x，V2的下限为V2x；
并设定模拟退火算法中初始温度T0，退火因子λ；
步骤504、采用模拟退火算法，输入步骤503中的参数，得到第k个海底节点的迭代的第一个新解P1(Xk1,Yk1,Zk1)、地震波在海水中的传播速度的第一次迭代值V11、地震波在海底介质中的传播速度的第一次迭代值V21；其中，Xk1,Yk1,Zk1分别表示第k个海底节点的第一个新解在X轴方向、Y轴方向、Z轴方向的坐标分量；
步骤505、将P1(Xk1,Yk1,Zk1)、V11和V21分别代入步骤402和步骤404，得到直达波第一次迭代的误差Qd(k1)和折射波第一次迭代的误差Qz(k1)；
步骤506、根据公式
Figure BDA0003308593480000061
和/>
Figure BDA0003308593480000062
得到直达波加权因子w1的第一次迭代新解w1(1)，折射波加权因子w2的第一次迭代新解w2(1)；
步骤507、将P1(Xk1,Yk1,Zk1)、V11、V21、w1(1)和w2(1)代入步骤405，得到第一次迭代的目标函数值Q(k1)；
步骤508、当Q(k1)<Q(k0)时，则接受新解；否则，按照Metropol is准则接受新解；
步骤509、根据T1＝T0×λ，得到第一次迭代退火温度T1，完成第一次模拟退火迭代；
步骤50A、将第k个海底节点的迭代的第一个新解P1(Xk1,Yk1,Zk1)、地震波在海水中的传播速度的第一次迭代值V11、地震波在海底介质中的传播速度的第一次迭代值V21、直达波加权因子w1的第一次迭代新解w1(1)，折射波加权因子w2的第一次迭代新解w2(1)输入，并按照步骤504至步骤509所述的方法，进行第二次模拟退火迭代；
步骤50B、按照步骤50A所述的方法，将第k个海底节点的迭代的第e个新解Pe(Xke,Yke,Zke)、地震波在海水中的传播速度的第e次迭代值V1e、地震波在海底介质中的传播速度的第e次迭代值V2e、直达波加权因子w1的第e次迭代新解w1(e)，折射波加权因子w2的第e次迭代新解w2(e)输入，并按照步骤504至步骤509所述的方法，进行第e+1次模拟退火迭代；其中，e为正整数；
步骤50C、多次重复步骤50B，完成第E次模拟退火迭代，当E满足设定的模拟退火迭代次数，且第E次迭代退火温度TE小于设定的温度值时，完成模拟退火迭代，并将进行第E次模拟退火迭代得到的新解记作目标函数Q(k)的实际坐标的最优解。
上述的一种浅海海底节点二次定位方法，其特征在于：步骤503中初始温度T0为0.1℃，退火因子λ的取值为[0.9，1)。
步骤50C中设定的模拟退火迭代次数为500～1000；
步骤50C中设定的温度值为10-4。
本发明与现有技术相比具有以下优点：
1、本发明步骤简单、设计合理且提高了海底节点二次定位的精度。
2、本发明使用效果好，首先是炮点和海底节点的布设及地震波数据的采集，其次是初至波初至时间的拾取及分离，然后是参与定位的有效炮点的选择，接着是利用直达波初至和折射波初至进行海底节点二次定位目标函数的获取，最后采用模拟退火算法，得到目标函数Q(k)的实际坐标的最优解，并将的实际坐标最优解中P*(Xk*,Yk*,Zk*)记作第k个海底节点二次定位位置，提高了海底节点二次定位位置的准确性，便于后续处理。
3、本发明在浅海区域由于海水较浅，可以准确识别的直达波仅限于检波点偏移距较小的几个炮点，随着炮点和检波点偏移距增大，海底产生的折射波成为初至波，这就导致折射波可利用的数据量远远超过直达波，利用直达波初至和折射波初至联合进行二次定位，能够解决浅海仅利用直达波初至进行二次定位时存在的初至数量少、方位分布不均且受炮端干扰的问题。
本发明方法步骤简单，设计合理，利用直达波初至和折射波初至联合进行二次定位，准确且快速获取海底节点的位置，能够解决浅海仅利用直达波初至进行二次定位时存在的初至数量少、方位分布不均且受炮端干扰的问题，实用性强。
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
图1为本发明的方法流程框图。
图2a为本发明地震波数据图。
图2b为利用本发明设计位置对图2a作线性动校正的示意图。
图2c为利用本发明二次定位位置对图2a作线性动校正的示意图。
具体实施方式
如图1所示的一种浅海海底节点二次定位方法，包括以下步骤：
步骤一、炮点和海底节点的布设及地震波数据的采集：
步骤101、在海面上布设多条炮线，每条炮线上设置多个炮点；其中，炮线的数量为L，L条炮线按照震源船放炮的顺序依次标记为第1条炮线，…，第l条炮线，...，第L条炮线，l和L均为正整数，1≤l≤L，每条炮线上炮点的数量为N个，N个炮点位于同一直线上，N个炮点的编号按照震源船的运行方向由前至后依次标记为第1个炮点，...，第n个炮点，...，第N个炮点，n为炮点编号，n和N为正整数，且1≤n≤N；
步骤102、在海底设计多个海底检波点，并在每个海底检波点处投放海底节点；其中，海底检波点和海底节点的数量均为K，K个海底检波点位于同一直线上，K个海底检波点按照震源船的运行方向由前至后依次标记为第1个海底检波点，...，第k个海底检波点，...，第K个海底检波点，相应的K个海底节点依次标记为第1个海底节点，...，第k个海底节点，...，第K个海底节点，k和K均为正整数，且1≤k≤K；
步骤103、建立O-XYZ坐标系；其中，原点为O，X轴和Y轴形成的平面XOY处于海平面上，Z轴垂直平面XOY且向下为正；
步骤104、在O-XYZ坐标系下，设定第k个海底检波点的设计位置P(Xk,Yk,Zk)和第l条炮线中第n个炮点的坐标位置Q(xl,n,yl,n,zl,n)；其中，Xk,Yk,Zk分别表示第k个海底检波点在X轴方向、Y轴方向、Z轴方向的坐标分量；xl,n,yl,n,zl,n分别表示第l条炮线中第n个炮点在X轴方向、Y轴方向、Z轴方向的坐标分量；
步骤105、采用震源船依次在L条炮线上N个炮点激发地震波至海底的过程中，K个海底节点分别采集到L条炮线上N个炮点激发的地震波数据；
步骤二、初至波初至时间的拾取及分离：
步骤201、根据K个海底节点分别采集L条炮线上N个炮点激发的地震波数据，拾取第l个炮线上第n个炮点激发地震波至海底时第k个海底节点采集到的初至波初至时间Tk(l,n)；
步骤202、根据公式
Figure BDA0003308593480000101
得到出现折射波初至的临界距离d；其中，V1c表示地震波在海水中的传播速度测量值，V2c表示地震波在海底介质中的传播速度测量值，h表示海水深度；
步骤203、根据第k个海底节点采集到的初至波初至时间Tk(l,n)以及出现折射波初至的临界距离d，从步骤105中的地震波数据中获取第k个海底节点采集到的初至波为直达波所对应的多个炮点和第k个海底节点采集的初至波为折射波时所对应的多个炮点，并将第k个海底节点采集到的初至波为直达波所对应的多个炮点记作初至波为直达波集合，第k个海底节点采集的初至波为折射波时所对应的多个炮点记作初至波为折射波集合；
步骤三、参与定位的有效炮点的选择：
步骤301、将每一条炮线上相邻两个炮点之间的间距记作Ld；
步骤302、获取第k个海底节点的待选择炮点；
步骤303、剔除异常炮点获取第k个海底节点的有效炮点；
步骤304、根据步骤203中初至波为直达波集合和初至波为折射波集合，从步骤303中第k个海底节点的有效炮点中得到第k个海底节点采集的初至波为直达波的有效炮点和第k个海底节点采集的初至波为折射波的有效炮点；
步骤305、第k个海底节点采集的初至波为直达波的有效炮点按照炮线顺序和炮点顺序排序并依次标记为第1个直达波有效炮点，...，第p个直达波有效炮点，...，第P个直达波有效炮点；其中，p和P均为正整数，且1≤p≤P；
将第k个海底节点采集的初至波为折射波的有效炮点按照炮线顺序和炮点顺序排序并依次标记为第1个折射波有效炮点，...，第p′个折射波有效炮点，...，第P′个折射波有效炮点；其中，p′和P′均为正整数，且1≤p′≤P′；
步骤四、海底节点二次定位目标函数的获取：
步骤401、设定第k个海底节点的实际坐标为P′(Xk′,Yk′,Zk′)，根据公式
Figure BDA0003308593480000111
得到第p个直达波有效炮点的直达波到达第k个海底节点的理论初至时间td(p,k)；其中，xd(p),yd(p),zd(p)分别表示第p个直达波有效炮点在X轴方向、Y轴方向、Z轴方向的坐标分量；
步骤402、根据公式
Figure BDA0003308593480000112
得到P个直达波有效炮点的直达波到达第k个海底检波点的理论和实际初至时间的误差Qd(k)；其中，Td(p,k)表示第k个海底节点采集到第p个直达波有效炮点的初至波的实际初至时间；
步骤403、根据公式
Figure BDA0003308593480000113
得到第p′个折射波有效炮点的折射波到达第k个海底节点的理论初至时间tz(p′,k)；其中，xz(p′),yz(p′),zz(p′)分别表示第p′个折射波有效炮点在X轴方向、Y轴方向、Z轴方向的坐标分量，V1表示地震波在海水中的传播速度，V2表示地震波在海底介质中的传播速度；
步骤404、根据公式
Figure BDA0003308593480000114
得到P′个折射波有效炮点的折射波到达第k个海底检波点的理论和实际初至时间的误差Qz(k)；其中，Tz(p′,k)表示第k个海底节点采集到第p′个折射波有效炮点的初至波的实际初至时间；
步骤405、根据公式Q(k)＝w1×Qd(k)+w2×Qz(k)，得到第k个海底节点采集到的初至波的理论和实际初至时间的目标函数Q(k)；w1表示直达波加权因子，w2表示折射波加权因子；
步骤五、海底节点二次定位位置的获取：
采用模拟退火算法，得到目标函数Q(k)的实际坐标的最优解，并将实际坐标的最优解中P*(Xk*,Yk*,Zk*)记作第k个海底节点二次定位位置。
本实施例中，步骤302中获取第k个海底节点的待选择炮点，具体过程如下：
步骤3021、以第k个海底节点所对应的第k个海底检波点的设计位置在XOY平面的投影为圆心，以半径[d+10Ld,d+30Ld]画圆，得到第k个海底节点待选择的炮点区域，并将第k个海底节点待选择的炮点区域记作第k个待选择炮点区域；
步骤3022、将第k个待选择炮点区域中的各个炮点记作第k个海底节点的待选择炮点；
步骤3023、将第k个海底节点的待选择炮点按照炮线顺序和炮点顺序排序分别记作第k1个待选择炮点，...，第ka个待选择炮点，...，第kA个待选择炮点；其中，a和A均为正整数，且1≤a≤A，A表示第k个海底节点的待选择炮点的总数，k1，ka和kA均为正整数。
本实施例中，步骤303中剔除异常炮点获取第k个海底节点的有效炮点，具体过程如下：
步骤3031、将A个待选择炮点和第k个海底检波点的水平间距分别记作A个炮检距；根据步骤201中的初至时间得到A个待选择炮点分别激发地震波至海底时第k个海底节点采集到的初至波初至时间；
步骤3032、以炮检距为x轴，以初至时间为y轴，将步骤3031中的A个炮检距和初至时间进行绘制并拟合得到关系曲线；
步骤3033、将偏移关系曲线且偏移的初至时间大于100ms所对应的待选择炮点剔除，得到第k个海底节点的有效炮点。
本实施例中，步骤五中采用模拟退火算法，得到目标函数Q(k)的实际坐标的最优解，具体过程如下：
步骤501、设定第k个海底节点的迭代的初始位置坐标为第k个海底节点所对应的第k个海底检波点的设计位置并记作P0(Xk,Yk,Zk)，地震波在海水中的传播速度V1的初始值V10为V1c，地震波在海底介质中的传播速度V2的初始值V20为V2c；
步骤502、设定直达波加权因子w1的初始值为w1(0)，折射波加权因子w2的初始值为w2(0)，并将w1(0)、w2(0)、P0(Xk,Yk,Zk)、V1c和V2c代入步骤405，得到初始目标函数值Q(k0)；其中，w1(0)和w2(0)均是(0，1)的随机数，且w1(0)+w2(0)＝1；
步骤503、设定第k个海底节点的迭代上限位置坐标为P(XkH,YkH,ZkH)，第k个海底节点的迭代下限位置坐标为P(Xkx,Ykx,Zkx)，V1的上限为V1s，V2的上限为V2s，V1的下限为V1x，V2的下限为V2x；
并设定模拟退火算法中初始温度T0，退火因子λ；
步骤504、采用模拟退火算法，输入步骤503中的参数，得到第k个海底节点的迭代的第一个新解P1(Xk1,Yk1,Zk1)、地震波在海水中的传播速度的第一次迭代值V11、地震波在海底介质中的传播速度的第一次迭代值V21；其中，Xk1,Yk1,Zk1分别表示第k个海底节点的第一个新解在X轴方向、Y轴方向、Z轴方向的坐标分量；
步骤505、将P1(Xk1,Yk1,Zk1)、V11和V21分别代入步骤402和步骤404，得到直达波第一次迭代的误差Qd(k1)和折射波第一次迭代的误差Qz(k1)；
步骤506、根据公式
Figure BDA0003308593480000131
和/>
Figure BDA0003308593480000132
得到直达波加权因子w1的第一次迭代新解w1(1)，折射波加权因子w2的第一次迭代新解w2(1)；
步骤507、将P1(Xk1,Yk1,Zk1)、V11、V21、w1(1)和w2(1)代入步骤405，得到第一次迭代的目标函数值Q(k1)；
步骤508、当Q(k1)<Q(k0)时，则接受新解；否则，按照Metropol is准则接受新解；
步骤509、根据T1＝T0×λ，得到第一次迭代退火温度T1，完成第一次模拟退火迭代；
步骤50A、将第k个海底节点的迭代的第一个新解P1(Xk1,Yk1,Zk1)、地震波在海水中的传播速度的第一次迭代值V11、地震波在海底介质中的传播速度的第一次迭代值V21、直达波加权因子w1的第一次迭代新解w1(1)，折射波加权因子w2的第一次迭代新解w2(1)输入，并按照步骤504至步骤509所述的方法，进行第二次模拟退火迭代；
步骤50B、按照步骤50A所述的方法，将第k个海底节点的迭代的第e个新解Pe(Xke,Yke,Zke)、地震波在海水中的传播速度的第e次迭代值V1e、地震波在海底介质中的传播速度的第e次迭代值V2e、直达波加权因子w1的第e次迭代新解w1(e)，折射波加权因子w2的第e次迭代新解w2(e)输入，并按照步骤504至步骤509所述的方法，进行第e+1次模拟退火迭代；其中，e为正整数；
步骤50C、多次重复步骤50B，完成第E次模拟退火迭代，当E满足设定的模拟退火迭代次数，且第E次迭代退火温度TE小于设定的温度值时，完成模拟退火迭代，并将进行第E次模拟退火迭代得到的新解记作目标函数Q(k)的实际坐标的最优解。
本实施例中，步骤503中初始温度T0为0.1℃，退火因子λ的取值为[0.9，1)。
步骤50C中设定的模拟退火迭代次数为500～1000；
步骤50C中设定的温度值为10-4。
本实施例中，需要说明的是，Td(p,k)和Tz(p′,k)均可通过步骤201中的初至时间获取。
本实施例中，XkH,YkH,ZkH分别表示第k个海底节点的迭代上限位置坐标在X轴方向、Y轴方向、Z轴方向的坐标分量，Xkx,Ykx,Zkx分别表示第k个海底节点的迭代下限位置坐标在X轴方向、Y轴方向、Z轴方向的坐标分量。
本实施例中，炮线的总数为L＝21。
本实施例中，每条炮线上炮点的数量N的取值为500～900。
本实施例中，相邻两个炮点之间的间距Ld的取值范围为30米～40米。
本实施例中，K＝100，相邻两个海底检波点之间的间距为100米～110米。
本实施例中，浅海是指水深在200米以内的海域。
本实施例中，针对第100个海底节点进行二次定位，第100个海底节点所对应的第100个海底检波点的设计位置(1205，1469，80)作为第100个海底节点的迭代的初始位置坐标，V10＝1500，V20＝1700；
第k个海底节点的迭代最大位置坐标为(1300，1600，85),V1s＝1550，V2s＝1750；第k个海底节点的迭代下限位置坐标为(1100，1400，80)，V1x＝1450，V2x＝1650。
本实施例中，经过步骤五得到第100个海底节点的二次定位位置为1115.90，1516.13，81.82)。
本实施例中，如图2a、2b、2c分别表示原始地震波数据、设计位置动校正结果、二次定位位置动校正结果，纵横为时间，横轴为任一个炮线上的炮点序号，可以看出，使用二次定位位置结果做线性动校正，初至波的形态趋于直线，则说明海底节点的二次定位位置准确。
本发明方法步骤简单，设计合理，利用直达波初至和折射波初至联合进行二次定位，准确且快速获取海底节点的位置，能够解决浅海仅利用直达波初至进行二次定位时存在的初至数量少、方位分布不均且受炮端干扰的问题，实用性强。
以上所述，仅是本发明的较佳实施例，并非对本发明作任何限制，凡是根据本发明技术实质对以上实施例所作的任何简单修改、变更以及等效结构变化，均仍属于本发明技术方案的保护范围内。

Claims (5)

1.一种浅海海底节点二次定位方法，其特征在于，该方法包括以下步骤：
步骤一、炮点和海底节点的布设及地震波数据的采集：
步骤101、在海面上布设多条炮线，每条炮线上设置多个炮点；其中，炮线的数量为L，L条炮线按照震源船放炮的顺序依次标记为第1条炮线，…，第l条炮线，...，第L条炮线，l和L均为正整数，1≤l≤L，每条炮线上炮点的数量为N个，N个炮点位于同一直线上，N个炮点的编号按照震源船的运行方向由前至后依次标记为第1个炮点，...，第n个炮点，...，第N个炮点，n为炮点编号，n和N为正整数，且1≤n≤N；
步骤102、在海底设计多个海底检波点，并在每个海底检波点处投放海底节点；其中，海底检波点和海底节点的数量均为K，K个海底检波点位于同一直线上，K个海底检波点按照震源船的运行方向由前至后依次标记为第1个海底检波点，...，第k个海底检波点，...，第K个海底检波点，相应的K个海底节点依次标记为第1个海底节点，...，第k个海底节点，...，第K个海底节点，k和K均为正整数，且1≤k≤K；
步骤103、建立O-XYZ坐标系；其中，原点为O，X轴和Y轴形成的平面XOY处于海平面上，Z轴垂直平面XOY且向下为正；
步骤104、在O-XYZ坐标系下，设定第k个海底检波点的设计位置P(Xk,Yk,Zk)和第l条炮线中第n个炮点的坐标位置Q(xl,n,yl,n,zl,n)；其中，Xk,Yk,Zk分别表示第k个海底检波点在X轴方向、Y轴方向、Z轴方向的坐标分量；xl,n,yl,n,zl,n分别表示第l条炮线中第n个炮点在X轴方向、Y轴方向、Z轴方向的坐标分量；
步骤105、采用震源船依次在L条炮线上N个炮点激发地震波至海底的过程中，K个海底节点分别采集到L条炮线上N个炮点激发的地震波数据；
步骤二、初至波初至时间的拾取及分离：
步骤201、根据K个海底节点分别采集L条炮线上N个炮点激发的地震波数据，拾取第l个炮线上第n个炮点激发地震波至海底时第k个海底节点采集到的初至波初至时间Tk(l,n)；
步骤202、根据公式
Figure FDA0004176513660000021
得到出现折射波初至的临界距离d；其中，V1c表示地震波在海水中的传播速度测量值，V2c表示地震波在海底介质中的传播速度测量值，h表示海水深度；
步骤203、根据第k个海底节点采集到的初至波初至时间Tk(l,n)以及出现折射波初至的临界距离d，从步骤105中的地震波数据中获取第k个海底节点采集到的初至波为直达波所对应的多个炮点和第k个海底节点采集的初至波为折射波时所对应的多个炮点，并将第k个海底节点采集到的初至波为直达波所对应的多个炮点记作初至波为直达波集合，第k个海底节点采集的初至波为折射波时所对应的多个炮点记作初至波为折射波集合；
步骤三、参与定位的有效炮点的选择：
步骤301、将每一条炮线上相邻两个炮点之间的间距记作Ld；
步骤302、获取第k个海底节点的待选择炮点；
步骤303、剔除异常炮点获取第k个海底节点的有效炮点；
步骤304、根据步骤203中初至波为直达波集合和初至波为折射波集合，从步骤303中第k个海底节点的有效炮点中得到第k个海底节点采集的初至波为直达波的有效炮点和第k个海底节点采集的初至波为折射波的有效炮点；
步骤305、第k个海底节点采集的初至波为直达波的有效炮点按照炮线顺序和炮点顺序排序并依次标记为第1个直达波有效炮点，...，第p个直达波有效炮点，...，第P个直达波有效炮点；其中，p和P均为正整数，且1≤p≤P；
将第k个海底节点采集的初至波为折射波的有效炮点按照炮线顺序和炮点顺序排序并依次标记为第1个折射波有效炮点，...，第p′个折射波有效炮点，...，第P′个折射波有效炮点；其中，p′和P′均为正整数，且1≤p′≤P′；
步骤四、海底节点二次定位目标函数的获取：
步骤401、设定第k个海底节点的实际坐标为P′(X′k,Y′k,Z′k)，根据公式
Figure FDA0004176513660000031
得到第p个直达波有效炮点的直达波到达第k个海底节点的理论初至时间td(p,k)；其中，xd(p),yd(p),zd(p)分别表示第p个直达波有效炮点在X轴方向、Y轴方向、Z轴方向的坐标分量；
步骤402、根据公式
Figure FDA0004176513660000032
得到P个直达波有效炮点的直达波到达第k个海底检波点的理论和实际初至时间的误差Qd(k)；其中，Td(p,k)表示第k个海底节点采集到第p个直达波有效炮点的初至波的实际初至时间；
步骤403、根据公式
Figure FDA0004176513660000033
得到第p′个折射波有效炮点的折射波到达第k个海底节点的理论初至时间tz(p′,k)；其中，xz(p′),yz(p′),zz(p′)分别表示第p′个折射波有效炮点在X轴方向、Y轴方向、Z轴方向的坐标分量，V1表示地震波在海水中的传播速度，V2表示地震波在海底介质中的传播速度；
步骤404、根据公式
Figure FDA0004176513660000034
得到P′个折射波有效炮点的折射波到达第k个海底检波点的理论和实际初至时间的误差Qz(k)；其中，Tz(p′,k)表示第k个海底节点采集到第p′个折射波有效炮点的初至波的实际初至时间；
步骤405、根据公式Q(k)＝w1×Qd(k)+w2×Qz(k)，得到第k个海底节点采集到的初至波的理论和实际初至时间的目标函数Q(k)；w1表示直达波加权因子，w2表示折射波加权因子；
步骤五、海底节点二次定位位置的获取：
采用模拟退火算法，得到目标函数Q(k)的实际坐标的最优解，并将实际坐标的最优解中P*(Xk*,Yk*,Zk*)记作第k个海底节点二次定位位置。
2.按照权利要求1所述的一种浅海海底节点二次定位方法，其特征在于：步骤302中获取第k个海底节点的待选择炮点，具体过程如下：
步骤3021、以第k个海底节点所对应的第k个海底检波点的设计位置在XOY平面的投影为圆心，以半径[d+10Ld,d+30Ld]画圆，得到第k个海底节点待选择的炮点区域，并将第k个海底节点待选择的炮点区域记作第k个待选择炮点区域；
步骤3022、将第k个待选择炮点区域中的各个炮点记作第k个海底节点的待选择炮点；
步骤3023、将第k个海底节点的待选择炮点按照炮线顺序和炮点顺序排序分别记作第k1个待选择炮点，...，第ka个待选择炮点，...，第kA个待选择炮点；其中，a和A均为正整数，且1≤a≤A，A表示第k个海底节点的待选择炮点的总数，k1，ka和kA均为正整数。
3.按照权利要求2所述的一种浅海海底节点二次定位方法，其特征在于：步骤303中剔除异常炮点获取第k个海底节点的有效炮点，具体过程如下：
步骤3031、将A个待选择炮点和第k个海底检波点的水平间距分别记作A个炮检距；根据步骤201中的初至时间得到A个待选择炮点分别激发地震波至海底时第k个海底节点采集到的初至波初至时间；
步骤3032、以炮检距为x轴，以初至时间为y轴，将步骤3031中的A个炮检距和初至时间进行绘制并拟合得到关系曲线；
步骤3033、将偏移关系曲线且偏移的初至时间大于100ms所对应的待选择炮点剔除，得到第k个海底节点的有效炮点。
4.按照权利要求1所述的一种浅海海底节点二次定位方法，其特征在于：步骤五中采用模拟退火算法，得到目标函数Q(k)的实际坐标的最优解，具体过程如下：
步骤501、设定第k个海底节点的迭代的初始位置坐标为第k个海底节点所对应的第k个海底检波点的设计位置并记作P0(Xk,Yk,Zk)，地震波在海水中的传播速度V1的初始值V10为V1c，地震波在海底介质中的传播速度V2的初始值V20为V2c；
步骤502、设定直达波加权因子w1的初始值为w1(0)，折射波加权因子w2的初始值为w2(0)，并将w1(0)、w2(0)、P0(Xk,Yk,Zk)、V1c和V2c代入步骤405，得到初始目标函数值Q(k0)；其中，w1(0)和w2(0)均是(0，1)的随机数，且w1(0)+w2(0)＝1；
步骤503、设定第k个海底节点的迭代上限位置坐标为P(XkH,YkH,ZkH)，第k个海底节点的迭代下限位置坐标为P(Xkx,Ykx,Zkx)，V1的上限为V1s，V2的上限为V2s，V1的下限为V1x，V2的下限为V2x；
并设定模拟退火算法中初始温度T0，退火因子λ；
步骤504、采用模拟退火算法，输入步骤503中的参数，得到第k个海底节点的迭代的第一个新解P1(Xk1,Yk1,Zk1)、地震波在海水中的传播速度的第一次迭代值V11、地震波在海底介质中的传播速度的第一次迭代值V21；其中，Xk1,Yk1,Zk1分别表示第k个海底节点的第一个新解在X轴方向、Y轴方向、Z轴方向的坐标分量；
步骤505、将P1(Xk1,Yk1,Zk1)、V11和V21分别代入步骤402和步骤404，得到直达波第一次迭代的误差Qd(k1)和折射波第一次迭代的误差Qz(k1)；
步骤506、根据公式
Figure FDA0004176513660000051
和/>
Figure FDA0004176513660000052
得到直达波加权因子w1的第一次迭代新解w1(1)，折射波加权因子w2的第一次迭代新解w2(1)；
步骤507、将P1(Xk1,Yk1,Zk1)、V11、V21、w1(1)和w2(1)代入步骤405，得到第一次迭代的目标函数值Q(k1)；
步骤508、当Q(k1)<Q(k0)时，则接受新解；否则，按照Metropolis准则接受新解；
步骤509、根据T1＝T0×λ，得到第一次迭代退火温度T1，完成第一次模拟退火迭代；
步骤50A、将第k个海底节点的迭代的第一个新解P1(Xk1,Yk1,Zk1)、地震波在海水中的传播速度的第一次迭代值V11、地震波在海底介质中的传播速度的第一次迭代值V21、直达波加权因子w1的第一次迭代新解w1(1)，折射波加权因子w2的第一次迭代新解w2(1)输入，并按照步骤504至步骤509所述的方法，进行第二次模拟退火迭代；
步骤50B、按照步骤50A所述的方法，将第k个海底节点的迭代的第e个新解Pe(Xke,Yke,Zke)、地震波在海水中的传播速度的第e次迭代值V1e、地震波在海底介质中的传播速度的第e次迭代值V2e、直达波加权因子w1的第e次迭代新解w1(e)，折射波加权因子w2的第e次迭代新解w2(e)输入，并按照步骤504至步骤509所述的方法，进行第e+1次模拟退火迭代；其中，e为正整数；
步骤50C、多次重复步骤50B，完成第E次模拟退火迭代，当E满足设定的模拟退火迭代次数，且第E次迭代退火温度TE小于设定的温度值时，完成模拟退火迭代，并将进行第E次模拟退火迭代得到的新解记作目标函数Q(k)的实际坐标的最优解。
5.按照权利要求4所述的一种浅海海底节点二次定位方法，其特征在于：步骤503中初始温度T0为0.1℃，退火因子λ的取值为[0.9，1)；
步骤50C中设定的模拟退火迭代次数为500～1000；
步骤50C中设定的温度值为10-4℃。
CN202111210053.2A 2021-10-18 2021-10-18 一种浅海海底节点二次定位方法 Expired - Fee Related CN113945981B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111210053.2A CN113945981B (zh) 2021-10-18 2021-10-18 一种浅海海底节点二次定位方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111210053.2A CN113945981B (zh) 2021-10-18 2021-10-18 一种浅海海底节点二次定位方法

Publications (2)

Publication Number Publication Date
CN113945981A CN113945981A (zh) 2022-01-18
CN113945981B true CN113945981B (zh) 2023-05-26

Family

ID=79331017

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111210053.2A Expired - Fee Related CN113945981B (zh) 2021-10-18 2021-10-18 一种浅海海底节点二次定位方法

Country Status (1)

Country Link
CN (1) CN113945981B (zh)

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101852867A (zh) * 2009-04-03 2010-10-06 中国石油集团东方地球物理勘探有限责任公司 一种矢量合成检波点二次定位方法

Family Cites Families (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101930080B (zh) * 2009-06-26 2012-08-15 中国石油集团东方地球物理勘探有限责任公司 曲面拟合海底电缆二次定位方法
EP2607920B1 (en) * 2011-12-19 2018-04-11 Sercel Method and device for estimating an inter-node distance between nodes arranged along towed acoustic linear antennas.
WO2017024536A1 (zh) * 2015-08-11 2017-02-16 深圳朝伟达科技有限公司 一种自动切除地震波波至的方法
NO342823B1 (en) * 2016-10-19 2018-08-13 Polarcus Dmcc Positioning of seismic equipment
CN108387867A (zh) * 2018-02-09 2018-08-10 华南理工大学 一种水下源节点定位方法
NO345833B1 (en) * 2019-08-13 2021-08-30 Magseis Fairfield ASA Ocean bottom seismic node deployment and retrieval system
CN112649875B (zh) * 2019-10-11 2024-08-27 中国石油天然气集团有限公司 海底节点检波器定位方法及装置
CN110687606B (zh) * 2019-10-25 2021-04-20 长安大学 一种海底节点地震仪三分量定向校正方法
CN110780350B (zh) * 2019-10-25 2021-10-01 中国海洋大学 一种海底地震检波器动态定位方法、系统及检波器
EP4115212A4 (en) * 2020-03-06 2024-03-20 Services Pétroliers Schlumberger MARINE SEISMIC IMAGING
CN111538089B (zh) * 2020-05-08 2023-04-07 中国石油天然气集团有限公司 海洋节点勘探多船多源震源阵列同步控制方法和装置
CN111580165A (zh) * 2020-05-27 2020-08-25 中国科学院地质与地球物理研究所 一种海底地震仪到时差定位装置及方法
CN112147695B (zh) * 2020-09-30 2022-09-16 长安大学 一种海底节点检波器水下姿态定向方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101852867A (zh) * 2009-04-03 2010-10-06 中国石油集团东方地球物理勘探有限责任公司 一种矢量合成检波点二次定位方法

Also Published As

Publication number Publication date
CN113945981A (zh) 2022-01-18

Similar Documents

Publication Publication Date Title
CN110780350B (zh) 一种海底地震检波器动态定位方法、系统及检波器
CN103477244B (zh) 检测水下结构的结构变化
CN107664758B (zh) 基于长基线或超短基线组网的深海导航定位系统及方法
CN112505750B (zh) 一种深拖多道地震拖缆姿态确定方法及处理终端
CN105676181A (zh) 基于分布式传感器能量比的水下运动目标扩展卡尔曼滤波跟踪方法
CN106443776B (zh) 一种基于时间切片法的海底地震仪重定位方法
CN110146924B (zh) 基于水波初至偏振方位的海底地震仪位置及方位反演方法
CN109738902A (zh) 一种基于同步信标模式的水下高速目标高精度自主声学导航方法
CN108761470B (zh) 一种基于拖缆形态方程解析的目标定位方法
CN111551942B (zh) 一种基于解卷积算法的水下自主航行器回坞方法
CN114460649A (zh) 一种深海近底拖曳式多道地震接收阵列形态重建方法
CN1226318A (zh) 定位水听器的方法
CN109212598B (zh) 基于直达波反演的三维空间二次定位方法
CN114779340A (zh) 一种深拖式高分辨率多道地震数据精细处理方法
CN110441736B (zh) 多关节水下无人潜行器变基线三维空间定位方法
CN116840786A (zh) 一种面向长基线水下定位系统的声线修正方法
CN109116416B (zh) 基于三维方式的倾斜海底检波器二次定位方法
CN113945981B (zh) 一种浅海海底节点二次定位方法
CN115308800B (zh) 一种利用海底反射波走时和地形数据定位海底地震仪的方法及处理终端
CN112147695B (zh) 一种海底节点检波器水下姿态定向方法
CN113703050B (zh) 一种深海地震垂直缆二次定位方法
CN113075732B (zh) 一种消除高分辨率小多道地震地层异常起伏的方法
CN114675331A (zh) 水面走航式探测海底气泡型浅层气的装置及方法
CN110243939B (zh) 一种基于机械波定位的水坝缺陷检测系统及其工作方法
CN117368849A (zh) 一种基于炮点差分的全息二次定位方法及其应用

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
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20230526