CN106441386B - 基于分布式光纤传感系统的数据处理方法及装置 - Google Patents
基于分布式光纤传感系统的数据处理方法及装置 Download PDFInfo
- Publication number
- CN106441386B CN106441386B CN201610861568.1A CN201610861568A CN106441386B CN 106441386 B CN106441386 B CN 106441386B CN 201610861568 A CN201610861568 A CN 201610861568A CN 106441386 B CN106441386 B CN 106441386B
- Authority
- CN
- China
- Prior art keywords
- scheme
- data
- bilateral filtering
- filtering algorithm
- standard deviation
- 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
Links
- 239000013307 optical fiber Substances 0.000 title claims abstract description 44
- 238000003672 processing method Methods 0.000 title claims abstract description 26
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 50
- 238000001914 filtration Methods 0.000 claims abstract description 48
- 230000002146 bilateral effect Effects 0.000 claims abstract description 43
- 238000012545 processing Methods 0.000 claims abstract description 29
- 238000000034 method Methods 0.000 claims description 43
- 238000006243 chemical reaction Methods 0.000 claims description 9
- 238000013139 quantization Methods 0.000 claims description 7
- 238000004364 calculation method Methods 0.000 abstract description 4
- 238000010586 diagram Methods 0.000 description 16
- 239000000835 fiber Substances 0.000 description 14
- 230000006870 function Effects 0.000 description 10
- 230000000694 effects Effects 0.000 description 9
- 230000003287 optical effect Effects 0.000 description 9
- 238000005516 engineering process Methods 0.000 description 7
- 230000008569 process Effects 0.000 description 7
- 230000002093 peripheral effect Effects 0.000 description 4
- 238000005259 measurement Methods 0.000 description 3
- 238000012986 modification Methods 0.000 description 3
- 230000004048 modification Effects 0.000 description 3
- 101100398248 Arabidopsis thaliana KIN10 gene Proteins 0.000 description 2
- 101000835860 Homo sapiens SWI/SNF-related matrix-associated actin-dependent regulator of chromatin subfamily B member 1 Proteins 0.000 description 2
- 102100025746 SWI/SNF-related matrix-associated actin-dependent regulator of chromatin subfamily B member 1 Human genes 0.000 description 2
- 230000003321 amplification Effects 0.000 description 2
- 230000001174 ascending effect Effects 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 230000003993 interaction Effects 0.000 description 2
- 230000004807 localization Effects 0.000 description 2
- 239000000203 mixture Substances 0.000 description 2
- 238000003199 nucleic acid amplification method Methods 0.000 description 2
- 230000002441 reversible effect Effects 0.000 description 2
- 230000003595 spectral effect Effects 0.000 description 2
- 230000002269 spontaneous effect Effects 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 230000009466 transformation Effects 0.000 description 2
- 230000002411 adverse Effects 0.000 description 1
- 210000001367 artery Anatomy 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 238000004590 computer program Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000013016 damping Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 230000005611 electricity Effects 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 230000008447 perception Effects 0.000 description 1
- 230000010287 polarization Effects 0.000 description 1
- 238000011002 quantification Methods 0.000 description 1
- 239000000523 sample Substances 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 230000001568 sexual effect Effects 0.000 description 1
- 210000003462 vein Anatomy 0.000 description 1
- 238000005406 washing Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01D—MEASURING NOT SPECIALLY ADAPTED FOR A SPECIFIC VARIABLE; ARRANGEMENTS FOR MEASURING TWO OR MORE VARIABLES NOT COVERED IN A SINGLE OTHER SUBCLASS; TARIFF METERING APPARATUS; MEASURING OR TESTING NOT OTHERWISE PROVIDED FOR
- G01D5/00—Mechanical means for transferring the output of a sensing member; Means for converting the output of a sensing member to another variable where the form or nature of the sensing member does not constrain the means for converting; Transducers not specially adapted for a specific variable
- G01D5/26—Mechanical means for transferring the output of a sensing member; Means for converting the output of a sensing member to another variable where the form or nature of the sensing member does not constrain the means for converting; Transducers not specially adapted for a specific variable characterised by optical transfer means, i.e. using infrared, visible, or ultraviolet light
- G01D5/32—Mechanical means for transferring the output of a sensing member; Means for converting the output of a sensing member to another variable where the form or nature of the sensing member does not constrain the means for converting; Transducers not specially adapted for a specific variable characterised by optical transfer means, i.e. using infrared, visible, or ultraviolet light with attenuation or whole or partial obturation of beams of light
- G01D5/34—Mechanical means for transferring the output of a sensing member; Means for converting the output of a sensing member to another variable where the form or nature of the sensing member does not constrain the means for converting; Transducers not specially adapted for a specific variable characterised by optical transfer means, i.e. using infrared, visible, or ultraviolet light with attenuation or whole or partial obturation of beams of light the beams of light being detected by photocells
- G01D5/353—Mechanical means for transferring the output of a sensing member; Means for converting the output of a sensing member to another variable where the form or nature of the sensing member does not constrain the means for converting; Transducers not specially adapted for a specific variable characterised by optical transfer means, i.e. using infrared, visible, or ultraviolet light with attenuation or whole or partial obturation of beams of light the beams of light being detected by photocells influencing the transmission properties of an optical fibre
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01D—MEASURING NOT SPECIALLY ADAPTED FOR A SPECIFIC VARIABLE; ARRANGEMENTS FOR MEASURING TWO OR MORE VARIABLES NOT COVERED IN A SINGLE OTHER SUBCLASS; TARIFF METERING APPARATUS; MEASURING OR TESTING NOT OTHERWISE PROVIDED FOR
- G01D3/00—Indicating or recording apparatus with provision for the special purposes referred to in the subgroups
- G01D3/028—Indicating or recording apparatus with provision for the special purposes referred to in the subgroups mitigating undesired influences, e.g. temperature, pressure
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Optical Communication System (AREA)
Abstract
本发明提供了一种基于分布式光纤传感系统的数据处理方法及装置，属于数据处理技术领域。该数据处理方法包括：获取原始数据，根据预设规则将所述原始数据转换为二维图；获取所述二维图的噪声标准差；根据所述噪声标准差得到所述二维图的灰度标准差；根据二维双边滤波算法对所述二维图进行滤波，获得去噪后的二维图，其中，所述二维双边滤波算法包括所述灰度标准差和预设的距离权重。相比于现有技术，本发明提供的数据处理方法及装置有效地提高了原始数据的信噪比，且降低了计算量，更加简单实用。
Description
技术领域
本发明涉及分布式光纤传感技术领域，具体而言，涉及一种基于分布式光纤传感系统的数据处理方法及装置。
背景技术
近年来，随着油气管道、高铁、大型建筑等的高速发展，其安全越来越受到各界关注。而分布式光纤传感技术由于其本身大量的优点，成为长距离、恶劣环境下进行外界信息感知的关键技术。然而，较低的信噪比是分布式光纤传感技术中一直以来都存在的急需解决的问题。为了解决这一难题，很多研究人员做了大量工作，主要包括基于硬件的改进方法和基于去噪算法的改进方法。其中，基于硬件改进方法虽然能有效地提高系统信噪比，但是其技术难度大并且会大大增加系统成本，性价比较低。而相比于基于硬件改进方法，将图像去噪算法应用于分布式光纤传感技术中能在不改进硬件的前提下获得比较理想的去噪效果，在获得令人满意的去噪效果条件下拥有更好的性价比。然而，目前应用到分布式光纤传感技术中且能够提高分布式光纤传感系统输出的解调信号的信噪比的去噪算法参数设置复杂，计算难度较大。
发明内容
鉴于此，本发明的目的在于提供一种基于分布式光纤传感系统的数据处理方法及装置，以有效地改善上述问题。
为了实现上述目的，本发明采用的技术方案如下：
第一方面，本发明实施例提供了一种基于分布式光纤传感系统的数据处理方法，所述方法包括：获取原始数据，根据预设规则将所述原始数据转换为二维图；获取所述二维图的噪声标准差；根据所述噪声标准差得到所述二维图的灰度标准差；根据二维双边滤波算法对所述二维图进行滤波，获得去噪后的二维图，其中，所述二维双边滤波算法包括所述灰度标准差和预设的距离权重。
第二方面，本发明实施例还提供了一种基于分布式光纤传感系统的数据处理装置，所述装置包括：转换模块、第一获取模块、第二获取模块及滤波模块。其中，转换模块用于获取原始数据，根据预设规则将所述原始数据转换为二维图；第一获取模块用于获取所述二维图的噪声标准差；第二获取模块用于根据所述噪声标准差得到所述二维图的灰度标准差；滤波模块用于根据二维双边滤波算法对所述二维图进行滤波，获得去噪后的二维图，其中，所述二维双边滤波算法包括所述灰度标准差和预设的距离权重。
相比于现有技术，本发明实施例提供的基于分布式光纤传感系统的数据处理方法及装置将获取到的原始数据转换为二维图后，获取该二维图的噪声标准差，并根据所获取的噪声标准差自适应地得到上述二维图的灰度标准差，采用二维双边滤波算法对所获得的二维图进行滤波。将该数据处理方法应用于分布式光纤传感系统输出的解调信号的处理时，能够有效地提高解调信号的信噪比。此外，本发明实施例提供的二维双边滤波算法的权系数由距离权重和灰度权重组成，其中，灰度权重根据所获取二维图的灰度标准差调节，距离权重直接预先设置，有效地降低了计算量，更加简单实用。
为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。
附图说明
为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。
图1为本发明较佳实施例提供的计算机的方框示意图；
图2为本发明较佳实施例提供的一种基于分布式光纤传感系统的数据处理方法的方法流程图；
图3为本发明较佳实施例提供的一种基于分布式光纤传感系统的数据处理方法的步骤S210的流程图；
图4为本发明较佳实施例提供的一种基于分布式光纤传感系统的数据处理方法的步骤S220的流程图；
图5为本实施例提供的一种基于相位敏感光时域反射技术的分布式光纤传感系统的结构示意图；
图6为在图5所示的系统中的27.6km的传感光纤上加载单个振动时去噪前、后的解调结果示意图；
图7为在图5所示的系统中的27.6km的传感光纤上加载两个振动时去噪前、后的解调结果示意图；
图8为本发明较佳实施例提供的一种基于分布式光纤传感系统的数据处理装置的功能模块图；
图9为本发明较佳实施例提供的另一种基于分布式光纤传感系统的数据处理装置的功能模块图。
具体实施方式
下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。同时，在本发明的描述中，术语“第一”、“第二”等仅用于区分描述，而不能理解为指示或暗示相对重要性。
如图1所示，是本发明较佳的实施例提供的计算机的方框示意图。所述计算机100包括基于分布式光纤传感系统的数据处理装置110、存储器120、存储控制器130、处理器140、外设接口150、输入输出装置160及显示装置170。
所述存储器120、存储控制器130、处理器140、外设接口150、输入输出装置160、显示装置170各元件相互之间直接或间接地电性连接，以实现数据的传输或交互。例如，这些元件相互之间可通过一条或多条通讯总线或信号线实现电性连接。所述数据处理装置110包括至少一个可以软件或固件(firmware)的形式存储于所述存储器120中的软件功能模块。所述处理器140用于执行存储器120中存储的可执行模块，例如所述数据处理装置110包括的软件功能模块或计算机100程序。
其中，存储器120可以是，但不限于，随机存取存储器(Random Access Memory，RAM)，只读存储器(Read Only Memory，ROM)，可编程只读存储器(Programmable Read-OnlyMemory，PROM)，可擦除只读存储器(Erasable Programmable Read-Only Memory，EPROM)，电可擦除只读存储器(Electric Erasable Programmable Read-Only Memory，EEPROM)等。其中，存储器120用于存储程序，所述处理器140在接收到执行指令后，执行所述程序，前述本发明实施例任一实施例揭示的流过程定义的服务器所执行的方法可以应用于处理器140中，或者由处理器140实现。
处理器140可能是一种集成电路芯片，具有信号的处理能力。上述的处理器140可以是通用处理器，包括中央处理器(Central Processing Unit，简称CPU)、网络处理器(Network Processor，简称NP)等；还可以是数字信号处理器(DSP)、专用集成电路(ASIC)、现成可编程门阵列(FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件。可以实现或者执行本发明实施例中的公开的各方法、步骤及逻辑框图。通用处理器可以是微处理器或者该处理器140也可以是任何常规的处理器等。
所述外设接口150将各种输入/输出装置耦合至处理器140以及存储器120。在一些实施例中，外设接口150、处理器140以及存储控制器130可以在单个芯片中实现。在其他一些实例中，他们可以分别由独立的芯片实现。
输入输出装置160用于提供给用户输入数据实现用户与所述计算机100的交互。所述输入输出装置160可以是，但不限于，鼠标和键盘等。
显示装置170在计算机100与用户之间提供一个交互界面(例如用户操作界面)或用于显示图像数据给用户参考。在本实施例中，所述显示装置170可以是液晶显示器或触控显示器。若为触控显示器，其可为支持单点和多点触控操作的电容式触控屏或电阻式触控屏等。支持单点和多点触控操作是指触控显示器能感应到来自该触控显示器上一个或多个位置处同时产生的触控操作，并将该感应到的触控操作交由处理器140进行计算和处理。
相比于基于硬件改进方法，将图像去噪算法应用于分布式光纤传感技术中能在不改进硬件的前提下获得比较理想的去噪效果，在获得令人满意的去噪效果条件下拥有更好的性价比。然而，目前应用到分布式光纤传感技术中且能够提高分布式光纤传感系统输出的解调信号的信噪比的去噪算法参数设置复杂，计算难度较大。
例如，目前应用到分布式光纤传感技术中的去噪算法主要有以下几种：(1)基于小波变换类。该方法虽然能够提高测量数据的信噪比，但是这种方法参数设置比较复杂，实用性较差；(2)基于二维边界探测类。该方法能够实现的信噪比提高较小，不足以满足实际需求；(3)基于非局部平均算法。该方法虽然能够实现较高的信噪比提升，但是复杂度较高，需要的处理时间很长，对于动态量测量系统则无法满足实时性要求，可行性较差。
鉴于此，本发明实施例提供了一种基于分布式光纤传感系统的数据处理方法，应用于分布式光纤传感系统输出的解调信号的处理。如图2所示，所述方法至少包括以下步骤S210至步骤S240。
步骤S210，获取原始数据，根据预设规则将所述原始数据转换为二维图。
原始数据包括多组子数据，每组子数据包括按预设顺序排列的多个数据。将获取到的原始数据转换为二维图，即可以将对原始数据的去噪问题转换为对图像的去噪问题。其中，原始数据为分布式光纤传感系统输出的多组解调数据，其中，所述多组解调数据是按照预设时间间隔采样得到的，每一组解调数据即为上述子数据。此时，每一组子数据均包括多个对应于不同光纤距离处的光强值，且每一组子数据所包括的多个光强值按照光纤距离大小排列。
如图3所示，作为一种实施方式，根据预设规则将所述原始数据转换为二维图的方法包括步骤S211至步骤S213。
步骤S211，获取所述多组子数据中位于同一排列位置处的数据的均值获得均值序列。
每组子数据所包括的数据的排列方式及采样间隔均相同，多组子数据中位于同一排列位置处的数据即为对应于同一个采样点的数据。例如，原始数据包括三组子数据，分别用D1、D2及D3表示，D1为{A1，A2，A3，A4，A5}，D2为{B1，B2，B3，B4，B5}，D3为{C1，C2，C3，C4，C5}。其中，A1、B1及C1均对应于采样点P1，A2、B2及C2均对应于采样点P2，A3、B3及C3均对应于采样点P3，A4、B4及C4均对应于采样点P4，A5、B5及C5均对应于采样点P5。此时，对多组子数据中位于同一排列位置处的数据求均值，即为分别获取A1、B1及C1的均值AVE1，A2、B2及C2的均值AVE2，A3、B3及C3的均值AVE3，A4、B4及C4的均值AVE4以及A5、B5及C5的均值AVE5，从而得到均值序列{AVE1，AVE2，AVE3，AVE4，AVE5}。
步骤S212，获取每组子数据中的每个数据与所述均值序列中与该数据排列位置相同的均值的差值。
获取到均值序列后，多组子数据中对于同一排列位置处的数据的均值即为均值序列中位于相同排列位置的数据。获取每一组子数据中的每个数据与所述均值序列中与该数据排列位置相同的均值的差值的具体实施方式可以为：根据每个数据的排列位置获取该数据对应的均值，获取该数据与所获取到的均值的差值。例如，获取到的上述子数据D1{A1，A2，A3，A4，A5}中的每个数据与该数据对应的均值的差值分别为：A1-AVE1，A2-AVE2，A3-AVE3，A4-AVE4，A5-AVE5。
步骤S213，将所述差值量化到[0～255]区间内得到二维图。
其中，所述二维图为灰度图。可以理解的是，灰度图可以理解为一个a×b的像素点阵列，a表示行，b表示列，阵列的元素值用灰度值表示，即每个像素点均具有一个灰度值，且灰度值的取值范围为[0～255]，0表示黑色，255表示白色。每一个像素点在灰度图中的位置可以用(i，j)表示，其中，i表示第i行，j表示第j列。因此，为了完成原始数据向二维图的转换，需要将所述差值量化到[0～255]区间内。
本实施例中，将所述差值量化到[0～255]区间内的实施方式优选为：将步骤S212中获取到的差值取绝对值后，获取所有差值的绝对值的最大值，将所有差值的绝对值均除以该最大值后再乘以255，即可以将步骤S212中获取到的差值量化到[0～255]区间内，得到二维图。这种量化方式能够将每组子数据所包括的数据之间的差异性有效地体现在二维图中，所得到的二维图更加准确。
当然，将所述差值量化到[0～255]区间内的实施方式也可以为：将步骤S212中获取到的每一组子数据与对应的均值的差值取绝对值后进行排序，可以按照由小到大的顺序排列，也可以按照由大到小的顺序排列。根据用户的需要将排列好的数据划分为依次排列的256个区间，256个区间与[0～255]区间包括的256个整数值一一对应，将属于相同区间内的数据均转换为该区间对应的整数值，从而得到二维图。需要说明的是，相邻的区间对应的整数值也为相邻的。例如，当每一组子数据与对应的均值的差值取绝对值后按照由小到大的顺序排列时，所划分的区间分别为{K1，K2，K3，…，K256}，若K1对应于0，则K2对应于1，…，K256对应于255，若K1对应于255，则K2对应于254，…，K256对应于0。
步骤S220，获取所述二维图的噪声标准差。
由于图像处理技术中，二维双边滤波算法的最优灰度标准差与图像的噪声标准差之间呈正比例关系。因此，为了获取最优灰度标准差，需要先获取上述步骤S210所获得的二维图的噪声标准差。
如图4所示，作为一种实施方式，步骤S220可以包括步骤S221至步骤S223。
步骤S221，将所述二维图划分为多个区域，每个所述区域包括多个像素点。
其中，具体的划分规则可以根据二维图的像素数量即原始数据的数据量、用户所需要的信噪比及对处理时间的要求设置。例如，当原始数据包括100组子数据，且每组子数据包括30000个数据时，所获取的二维图可以包括100(行)×30000(列)个像素点，此时，考虑到所需要的信噪比及对处理时间，可以将20×20个像素点划分为一个区域。
步骤S222，获取每个所述区域所包括的多个像素点的灰度值的标准差；
获取步骤S221中获得的每个区域所包括的多个像素点的灰度值后，得到每个区域所包括的多个像素点的灰度值的标准差。
步骤S223，获取所有所述区域的标准差的均值作为所述二维图的噪声标准差。
步骤S230，根据所述噪声标准差得到所述二维图的灰度标准差。
图像处理技术中，二维双边滤波算法的最优灰度标准差与图像的噪声标准差之间呈正比例关系。因此，获取到二维图的噪声标准差后，根据所述噪声标准差及预设的正比例系数即可以获得该二维图的最优灰度标准差。
其中，预设的正比例系数可以根据多次试验获得。具体的，可以获取多组原始数据，每组原始数据对应于一个二维图，获取每一个二维图的噪声标准差。根据经验设置一个最优灰度标准差的最小值σmin。将该最小值作为最优灰度标准差，将当前二维图、预设的距离权重以及最优灰度标准差带入二维双边滤波算法即可得到去噪后的二维图，并计算去噪后的二位图的信噪比SNR1。以该最小值σmin为起点按照预设步长Δσ递增更新当前二维图的最优灰度标准差为σmin+Δσ，重复上述过程得到信噪比SNR2；再次按照预设步长Δσ更新当前二维图的最优灰度标准差为σmin+2Δσ，重复上述过程得到信噪比SNR3，依此类推，可以得到一系列的最优灰度标准差及其对应的去噪后的二维图的信噪比。以上述一系列的最优灰度标准差为横坐标，去噪后的二维图的信噪比为纵坐标，可以得到关于最优灰度标准差与去噪后的二维图的信噪比的曲线，根据具体试验可知，该曲线呈高斯分布，存在一个最优灰度标准差σp对应的信噪比为最大值，此时可以得到该二维图的噪声标准差σn对应的最优灰度标准差σp。
以获取到的每一个二维图分别作为当前二维图，重复上述过程即可以得到每一个二维图的噪声标准差对应的最优灰度标准差，此时，可以得到多个相互对应的噪声标准差及最优灰度标准差：(σn1，σp1)、(σn2，σp2)、(σn3，σp3)、…、(σnk，σpk)。将上述多个相互对应的噪声标准差及最优灰度标准差进行曲线拟合，即可以得到正比例系数ξ，同时，也验证了二维双边滤波算法的最优灰度标准差与由原始数据转换成的二维图的噪声标准差之间的正比例关系。
可以理解的是，在本实施例的另一种实施方式中，也可以根据二维图的噪声方差获取所述二维图的灰度标准差。二维图的噪声方差的获取方式与二维图的噪声标准差的获取方式类似，不同之处在于，上述步骤S222为获取每个所述区域所包括的多个像素点的灰度值的方差；上述步骤S223为获取所有所述区域的方差的均值作为所述二维图的噪声方差。
步骤S240，根据二维双边滤波算法对所述二维图进行滤波，获得去噪后的二维图。
其中，所述二维双边滤波算法包括所述灰度标准差和预设的距离权重。
数字图像处理中，滤波的基本原理是将图像中的每个点与其相邻的像素点的灰度值作加权平均卷积运算。相比于其它图像滤波算法，双边滤波算法能够在平滑图像的同时保持图像的边缘细节。双边滤波器在二维邻域即滤波窗口的加权系数是由距离权重和灰度权重两部分组成。具体来讲，距离权重和灰度权重的乘积即为双边滤波器在滤波窗口的加权系数。其中，距离权重是由像素的空间距离决定的，灰度权重是由像素间的亮度差值决定的。
可以理解的是，对于含有噪声的图像X中的像素点(i，j)处的灰度值，假设滤波窗口的像素范围为N×N，所述滤波窗口包括一个中心像素点和位于中心像素点的N×N邻域内的多个参考像素点。此时，采用现有的双边滤波算法对像素点(i，j)进行双边滤波后的输出值可以表示为：
式(1)中，Y表示去噪之后的图像；(i，j)表示滤波窗口的中心像素点O，m表示滤波窗口中每个参考像素点到中心像素点之间的横向距离，n表示滤波窗口中每个参考像素点到中心像素点之间的纵向距离。需要说明的是，上述横向距离和纵向距离均为欧式距离，由于滤波窗口为方形窗口，m和n的取值范围均为[-P，P]。W[i,j；m,n]为加权系数，W[i,j；m,n]＝WG[i,j；m,n]·WR[i,j；m,n]，其中，WG[i,j；m,n]表示距离权重，WR[i,j；m,n]表示灰度权重。
现有的二维双边滤波算法中，上述距离权重WG[i,j；m,n]和灰度权重WR[i,j；m,n]可以分别表示为：
其中，σG为距离标准差，σR为灰度标准差。由式(2)和式(3)可以看出，为了调节距离权重和灰度权重以调节滤波窗口的加权系数得到较好的滤波效果，需要同时调节距离标准差σG和灰度标准差σR，参数设置较为复杂，实际适用性较差。因此，本实施例中，通过预先设置距离权重，自适应地调节灰度标准差σR，以保证得到较好的滤波效果的同时简化参数设置。
考虑到距离权重WG[i,j；m,n]为高斯函数，WG的值随着参考像素点到中心像素点的距离的增加而呈现指数衰减，并且其衰减快慢由σG决定。本实施例中，忽略σG对WG的影响，将距离权重WG设置为1/(m2+n2+1)。
需要说明的是，当本实施例应用于分布式光纤传感系统输出的解调信号的处理，即上述原始数据为分布式光纤传感系统输出的解调数据时，假设该原始数据包括的每组子数据对应于二维图的一行像素点，携带外界扰动信息的数据通常主要存在于二维图的某一列像素点，而每组子数据中携带外界扰动信息的数据较少，即对应于二维图的每一行像素点中携带外界扰动信息的像素点较少。因此，本发明实施例提供的数据处理方法中，将二维双边滤波算法的滤波窗口为矩形窗口，即此时滤波窗口的像素范围为M×N，且M＜N。相比于现有的方形滤波窗口，能够在保证滤波效果的同时减少处理时间，使得本数据处理方法具有更好的实时性。
因此，本实施例中采用的二维双边滤波算法的表达式可以表示为：
式(4)中，[-P，P]表示滤波窗口M×N中横向距离m的取值范围，[-Q，Q]表示滤波窗口M×N中纵向距离n的取值范围，当原始数据包括的每组子数据对应于二维图的一行像素点时，P<Q。此时，将步骤S210获取到的二维图及步骤S230获取到的最优灰度标准差带入式(4)，根据式(4)所示的二维双边滤波算法对所述二维图进行滤波，即可以获得去噪后的二维图。
获得去噪后的二维图后，可以根据步骤S210将原始数据转换为二维图的规则将去噪后的二维图逆向转换为去噪后的原始数据。当然，对于分布式光纤传感系统来讲，可以直接根据去噪后的二维图定位扰动点。
可以理解的是，为了进一步提高去噪后的二维图的信噪比，可以对获取到的去噪后的二维图重复上述步骤S220至步骤S240。也就是说，再次获取去噪后的二维图的噪声标准差，根据所获取的噪声标准差得到所述去噪后的二维图的灰度标准差，根据式(4)所示的二维双边滤波算法对去噪后的二维图再次进行滤波，获得再次去噪后的二维图。同时考虑到去噪性能和处理时间，优选采用两次迭代以实现好的去噪性能和实时性同时兼顾。当然，在处理时间允许的条件下，也可以采用两次以上的迭代以达到更好的去噪性能。
为了更好的对本实施例的技术方案及效果进行说明，将本发明实施例提供的基于分布式光纤传感系统的数据处理方法应用于一种基于相位敏感光时域反射技术的分布式光纤传感系统中，用于实现该系统的解调信号的滤波，提高该解调信号的信噪比，以便于实现动态应变的测量。
如图5所示，上述基于相位敏感光时域反射技术的分布式光纤传感系统由光路和电路两部分组成。窄线宽激光器501输出的连续激光经第一光耦合器502后分为两路光信号。一路光信号由经脉冲发生器503与声光调制器驱动504驱动的声光调制器505调制为脉冲光之后经第一掺铒光纤放大器506放大后经环形器507注入到光纤光栅滤波器508。经光纤光栅滤波器508滤除自发辐射噪声后注入传感光纤509，传感光纤509中散射回来的后向瑞利散射光经第二掺铒光纤放大器510放大后进入光滤波器511，由光滤波器511滤除所述后向瑞利散射光中包含的自发辐射噪声。另一路光信号通过可调衰减器512和偏振控制器513调节后通过第二光耦合器514与上述后向散射信号耦合，然后经平衡光电探测器515进行光电转换，获得的电信号经过电低噪放大器516放大后由带通滤波器517滤除带外噪声，滤除带外噪声后的电信号接着采用3dB电功分器518分为两路，其中一路经可调延时线519调节之后与另一路信号在混频器520上进行混频解调，解调结果经过电低通滤波器521后可以通过数据采集卡522进行数据采集，将采集到的数据发送到计算机100后，采用本发明实施例提供的数据处理方法对采集到的信号进行数据处理，以提高采集到的数据的信噪比。
实际使用时，平衡光电探测器515电域带宽大于声光调制器505引入的频移值。带通滤波器517中心频率与声光调制器505引入的频移相同，带宽需大于窄线宽激光器501频率漂移值和探测脉冲脉宽对应的频率值中的较大值。电功分器518的功率比为50:50。
当在27.6km的传感光纤509的某一位置处加载单个振动时，将采集到的原始的解调信号转换为的二维图如图6(a)所示。通过本实施例提供的数据处理方法对图6(a)所示的二维图进行滤波处理后得到的去噪后的二维图如图6(b)所示。将去噪后的二维图逆向转换为去噪后的解调信号后，与原始的解调信号进行对比，信噪比显著提高。如图6(c)所示，原始的解调信号OR的信噪比为6.43dB，而去噪后的解调信号AF的信噪比为14.31dB。此外，图6(d)示出了采用本发明实施例提供的数据处理方法处理去噪前后上述基于相位敏感光时域反射技术的分布式光纤传感系统所获得的空间分辨率的对比图。以下将原始的解调信号简称为原始信号，采用本实施例提供的数据处理方法去噪后的解调信号简称为去噪信号。由图6(d)示出的去噪前的原始信号可得上述系统的空间分辨率约为7.2m，由图6(d)示出的去噪信号可得上述系统的空间分辨率约为6m。说明本实施例提供的数据处理方法不仅能够有效地提高上述分布式光纤传感系统输出的解调信号的信噪比，还能有效地提高系统的空间分辨率。
上述图6(a)和图6(b)中，横坐标为光纤距离，单位为千米(km)，纵坐标为振动轨迹，即可以每振动一次采集一组解调数据，即可以用振动轨迹表示所采集的解调数据的组数。上述图6(c)和图6(d)中，横坐标为光纤距离，单位为千米(km)，纵坐标为归一化强度。可以理解的是，随着光在光纤中传播距离的增加，由于光纤的色散性质及存在的其它损耗，光的强度随着传播距离的增加而衰减。为了避免光强的衰减对处理结果造成影响，需要对解调数据的光强进行归一化处理。
当在27.6km的传感光纤509上加载两个振动时，测得的解调数据采用本实施例提供的数据处理方法去噪前、后的结果如图7所示。其中，图7(a)和图7(b)示出了去噪前的原始信号转换的二维图；图7(c)示出了去噪前的原始信号的定位曲线和对应的信噪比；图7(d)和图7(e)示出了去噪后的二维图；图7(f)示出了去噪信号的定位曲线和对应的信噪比。
分别对比图7(a)和图7(d)、图7(b)和图7(e)可以明显看出，图7(d)和图7(e)中的噪声明显减小，提高了携带扰动信号的像素点的对比度。对比图7(c)和图7(f)可以明显看出，去噪前的原始信号对应的信噪比：SNR1为5.32dB，SNR2为5.14dB，去噪信号对应的信噪比：SNR1为12.44dB，SNR2为11.57dB，结果表明本实施例提供的数据处理方法能够有效地提高上述分布式光纤传感系统输出的解调信号的信噪比。
综上所述，本发明实施例提供的基于分布式光纤传感系统的数据处理方法将获取到的原始数据转换为二维图后，采用改进的二维双边滤波算法对所获得的二维图进行滤波。相比于现有技术，有效地提高了原始数据的信噪比，且有效地降低了计算量，更加简单实用。应用于分布式光纤传感系统输出的解调信号的处理时，不仅能够显著提高解调信号的信噪比，还能够提高系统的空间分辨率，从而有利于降低系统的误报率。
另外，本发明实施例还提供了一种基于分布式光纤传感系统的数据处理装置，如图8所示，所述数据处理装置110包括转换模块810、第一获取模块820、第二获取模块830和滤波模块840。
其中，转换模块810用于获取原始数据，根据预设规则将所述原始数据转换为二维图。其中，所述原始数据为分布式光纤传感系统输出的解调数据。第一获取模块820用于获取所述二维图的噪声标准差。第二获取模块830用于根据所述噪声标准差得到所述二维图的灰度标准差。滤波模块840用于根据二维双边滤波算法对所述二维图进行滤波，获得去噪后的二维图，其中，所述二维双边滤波算法包括所述灰度标准差和预设的距离权重。
其中，所述二维双边滤波算法的滤波窗口为矩形窗口。所述二维双边滤波算法的滤波窗口包括一个中心像素点和多个参考像素点，所述距离权重根据公式得到。其中，WG表示预设的距离权重，m表示所述二维双边滤波算法的滤波窗口中参考像素点到中心像素点的横向距离，n表示所述二维双边滤波算法的滤波窗口中参考像素点到中心像素点的纵向距离。
具体的，如图9所示，所述第一获取模块820包括划分单元821、第一获取单元822及第二获取单元823。其中，划分单元821用于将所述二维图划分为多个区域，每个所述区域包括多个像素点。第一获取单元822用于获取每个所述区域所包括的多个像素点的灰度值的标准差。第二获取单元823用于获取所有所述区域的标准差的均值作为所述二维图的噪声标准差。
所述原始数据包括多组子数据，每组子数据包括按预设顺序排列的多个数据。具体的，如图9所示，上述转换模块810包括均值序列获取单元811、差值获取单元812及量化单元813。
其中，均值序列获取单元811用于获取所述多组子数据中位于同一排列位置处的数据的均值获得均值序列。差值获取单元812用于获取每一组子数据中的每个数据与所述均值序列中与该数据排列位置相同的均值的差值。量化单元813用于将所述差值量化到[0～255]区间内得到二维图。
在本申请所提供的几个实施例中，应该理解到，所揭露的装置和方法，也可以通过其它的方式实现。以上所描述的装置实施例仅仅是示意性的，例如，附图中的流程图和框图显示了根据本发明的多个实施例的装置、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或代码的一部分，所述模块、程序段或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现方式中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。
另外，在本发明各个实施例中的各功能模块可以集成在一起形成一个独立的部分，也可以是各个模块单独存在，也可以两个或两个以上模块集成形成一个独立的部分。
所述功能如果以软件功能模块的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。
Claims (8)
1.一种基于分布式光纤传感系统的数据处理方法，其特征在于，所述方法包括：
获取原始数据，根据预设规则将所述原始数据转换为二维图；
获取所述二维图的噪声标准差；
根据所述噪声标准差得到所述二维图的灰度标准差；
根据二维双边滤波算法对所述二维图进行滤波，获得去噪后的二维图，其中，所述二维双边滤波算法包括所述灰度标准差和预设的距离权重，所述二维双边滤波算法的滤波窗口包括一个中心像素点和多个参考像素点，所述距离权重根据公式得到，其中，WG表示距离权重，m表示所述二维双边滤波算法的滤波窗口中所述参考像素点到所述中心像素点的横向距离，n表示所述二维双边滤波算法的滤波窗口中所述参考像素点到所述中心像素点的纵向距离。
2.根据权利要求1所述的方法，其特征在于，所述二维双边滤波算法的滤波窗口为矩形窗口。
3.根据权利要求1所述的方法，其特征在于，所述原始数据包括多组子数据，每组所述子数据包括按预设顺序排列的多个数据，所述的根据预设规则将所述原始数据转换为二维图的步骤，包括：
获取所述多组子数据中位于同一排列位置处的数据的均值获得均值序列；
获取每组所述子数据中的每个数据与所述均值序列中与该数据排列位置相同的均值的差值；
将所述差值量化到[0～255]区间内得到二维图。
4.根据权利要求1所述的方法，其特征在于，所述的获取所述二维图的噪声标准差的步骤，包括：
将所述二维图划分为多个区域，每个所述区域包括多个像素点；
获取每个所述区域所包括的多个像素点的灰度值的标准差；
获取所有所述区域的标准差的均值作为所述二维图的噪声标准差。
5.根据权利要求1所述的方法，其特征在于，所述的根据二维双边滤波算法对所述二维图进行滤波，获得去噪后的二维图的步骤之后，还包括：
获取去噪后的二维图的噪声标准差；
根据所获取的噪声标准差得到所述去噪后的二维图的灰度标准差；
根据所述二维双边滤波算法对所述去噪后的二维图再次进行滤波，获得再次去噪后的二维图。
6.一种基于分布式光纤传感系统的数据处理装置，其特征在于，包括：
转换模块，用于获取原始数据，根据预设规则将所述原始数据转换为二维图；
第一获取模块，用于获取所述二维图的噪声标准差；
第二获取模块，用于根据所述噪声标准差得到所述二维图的灰度标准差；
滤波模块，用于根据二维双边滤波算法对所述二维图进行滤波，获得去噪后的二维图，其中，所述二维双边滤波算法包括所述灰度标准差和预设的距离权重，所述二维双边滤波算法的滤波窗口包括一个中心像素点和多个参考像素点，所述距离权重根据公式得到，其中，WG表示距离权重，m表示所述二维双边滤波算法的滤波窗口中所述参考像素点到所述中心像素点的横向距离，n表示所述二维双边滤波算法的滤波窗口中所述参考像素点到所述中心像素点的纵向距离。
7.根据权利要求6所述的装置，其特征在于，所述二维双边滤波算法的滤波窗口为矩形窗口。
8.根据权利要求6所述的装置，其特征在于，所述第一获取模块包括：
划分单元，用于将所述二维图划分为多个区域，每个所述区域包括多个像素点；
第一获取单元，用于获取每个所述区域所包括的多个像素点的灰度值的标准差；
第二获取单元，用于获取所有所述区域的标准差的均值作为所述二维图的噪声标准差。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610861568.1A CN106441386B (zh) | 2016-09-29 | 2016-09-29 | 基于分布式光纤传感系统的数据处理方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610861568.1A CN106441386B (zh) | 2016-09-29 | 2016-09-29 | 基于分布式光纤传感系统的数据处理方法及装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106441386A CN106441386A (zh) | 2017-02-22 |
| CN106441386B true CN106441386B (zh) | 2019-02-26 |
Family
ID=58171057
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610861568.1A Active CN106441386B (zh) | 2016-09-29 | 2016-09-29 | 基于分布式光纤传感系统的数据处理方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106441386B (zh) |
Families Citing this family (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106908169B (zh) * | 2017-02-28 | 2020-03-17 | 中国科学院合肥物质科学研究院 | 一种分布式光纤温度传感器数据处理方法 |
| CN107270952B (zh) * | 2017-07-27 | 2020-03-31 | 天津求实飞博科技有限公司 | 基于光频域反射长距离光纤分布式扰动传感信号处理方法 |
| CN108665427A (zh) * | 2018-04-17 | 2018-10-16 | 浙江华睿科技有限公司 | 一种图像降噪方法及装置 |
| CN113237431B (zh) * | 2021-05-06 | 2022-03-18 | 山东大学 | 一种提升ofdr系统分布式空间分辨率的测量方法 |
| CN113188461B (zh) * | 2021-05-06 | 2022-05-17 | 山东大学 | 一种高空间分辨率下的ofdr大应变测量方法 |
| CN115001571B (zh) * | 2022-05-11 | 2023-09-12 | 青岛诺克通信技术有限公司 | 一种光时域反射计的自适应数字滤波方法及光时域反射计 |
| CN116399379B (zh) * | 2023-06-07 | 2023-11-03 | 山东省科学院激光研究所 | 分布式光纤声波传感系统及其测量方法 |
| CN116380140B (zh) * | 2023-06-07 | 2023-11-03 | 山东省科学院激光研究所 | 基于均值滤波技术的分布式声波传感系统及其测量方法 |
| CN116907677B (zh) * | 2023-09-15 | 2023-11-21 | 山东省科学院激光研究所 | 用于混凝土结构的分布式光纤温度传感系统及其测量方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101425655A (zh) * | 2007-10-31 | 2009-05-06 | 中国科学院半导体研究所 | 分布反馈光纤激光器调频、降噪的装置 |
| CN105913396A (zh) * | 2016-04-11 | 2016-08-31 | 湖南源信光电科技有限公司 | 一种噪声估计的图像边缘保持混合去噪方法 |
| CN105931402A (zh) * | 2016-06-27 | 2016-09-07 | 上海波汇科技股份有限公司 | 基于图像识别的光纤周界入侵监测方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| ES2392527B1 (es) * | 2011-05-13 | 2013-11-11 | Universidad Pública de Navarra | Dispositivo y procedimiento para la medida de la distribución de magnitudes físicas en una fibra óptica |
-
2016
- 2016-09-29 CN CN201610861568.1A patent/CN106441386B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101425655A (zh) * | 2007-10-31 | 2009-05-06 | 中国科学院半导体研究所 | 分布反馈光纤激光器调频、降噪的装置 |
| CN105913396A (zh) * | 2016-04-11 | 2016-08-31 | 湖南源信光电科技有限公司 | 一种噪声估计的图像边缘保持混合去噪方法 |
| CN105931402A (zh) * | 2016-06-27 | 2016-09-07 | 上海波汇科技股份有限公司 | 基于图像识别的光纤周界入侵监测方法 |
Non-Patent Citations (1)
| Title |
|---|
| 一种新的自适应双边滤波算法;余博等;《应用科学学报》;20120930;第30卷(第3期);正文的第518页左栏第2段至第522页右栏第2段 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN106441386A (zh) | 2017-02-22 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106441386B (zh) | 基于分布式光纤传感系统的数据处理方法及装置 | |
| JP6247362B2 (ja) | 容量型タッチセンサの信号対雑音を局所的に改善するための技法 | |
| Fang et al. | Noise reduction in lidar signal based on discrete wavelet transform | |
| Sanders et al. | Detecting edges in the X-ray surface brightness of galaxy clusters | |
| JP5819425B2 (ja) | 受信器と相互関係を有する容量型タッチセンサ | |
| CN103154870A (zh) | 具有码分和时分发送波形的电容式触摸传感器 | |
| CN107810466B (zh) | 触摸传感器中的多个相关 | |
| DE102016106511A1 (de) | Parametrische Online-Kalibrierung und Kompensation bei der TOF-Bilderzeugung | |
| Zhao et al. | Parameter optimal determination for canny edge detection | |
| CN108415073A (zh) | 角度域逆散射偏移成像方法及装置 | |
| CN108196242A (zh) | 基于边沿检测的激光雷达计时方法及数据处理单元 | |
| CN113419222A (zh) | 基于雷达信号提取桥梁振动频率的方法及系统 | |
| CN109300139A (zh) | 车道线检测方法和装置 | |
| Larsen et al. | Suppression of very low frequency radio noise in transient electromagnetic data with semi-tapered gates | |
| Shang et al. | Chaotic SVD method for minimizing the effect of exponential trends in detrended fluctuation analysis | |
| Cao et al. | Correction algorithm of the frequency-modulated continuous-wave LIDAR ranging system | |
| Wu et al. | Noise suppression of distributed fiber-optical acoustic sensing seismic data by attention-guided multiscale generative adversarial network | |
| CN107066364A (zh) | 性能检测方法及装置 | |
| Fröhlich et al. | Computation of decaying turbulence in an adaptive wavelet basis | |
| Naglič et al. | OpenCL framework for fast estimation of optical properties from spatial frequency domain images | |
| CN106291680A (zh) | 一种数据低频延拓方法 | |
| Raynaud et al. | Heterogeneous filtering of ensemble‐based background‐error variances | |
| Song et al. | InSAR phase filtering using spatially adapted total generalized variation | |
| Wang et al. | Multiple prediction from incomplete data with the focused curvelet transform | |
| CN112305590B (zh) | 粘声介质叠前时间偏移计算方法及装置 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |