WO2024057462A1 - 位相ｏｔｄｒにおける信号処理方法 - Google Patents

位相ｏｔｄｒにおける信号処理方法 Download PDF

Info

Publication number
WO2024057462A1
WO2024057462A1 PCT/JP2022/034477 JP2022034477W WO2024057462A1 WO 2024057462 A1 WO2024057462 A1 WO 2024057462A1 JP 2022034477 W JP2022034477 W JP 2022034477W WO 2024057462 A1 WO2024057462 A1 WO 2024057462A1
Authority
WO
WIPO (PCT)
Prior art keywords
scattered light
optical
vector
signal processing
rotation angle
Prior art date
Application number
PCT/JP2022/034477
Other languages
English (en)
French (fr)
Inventor
佳史 脇坂
央 高橋
大輔 飯田
優介 古敷谷
Original Assignee
日本電信電話株式会社
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 日本電信電話株式会社 filed Critical 日本電信電話株式会社
Priority to PCT/JP2022/034477 priority Critical patent/WO2024057462A1/ja
Publication of WO2024057462A1 publication Critical patent/WO2024057462A1/ja

Links

Images

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

Definitions

  • the present disclosure relates to a phase OTDR that measures the phase of scattered light from each point of an optical fiber to be measured.
  • DAS Distributed Acoustic Sensing
  • a method for measuring physical vibrations applied to an optical fiber distributed in the longitudinal direction of the optical fiber. is known (Non-Patent Document 1).
  • phase OTDR Optical Time Domain Reflectometry
  • DAS-phase DAS-P
  • the phase changes linearly as the optical path length of the optical fiber changes due to vibration, and the rate of change can be considered approximately the same at each point in the longitudinal direction of the optical fiber, making it possible to quantitatively measure vibration.
  • it is possible to faithfully reproduce the vibration waveform applied to the optical fiber to be measured for example, Non-Patent Document 2).
  • Non-Patent Document 3 Patent Document 1
  • This method utilizes the fact that the point at which the intensity of scattered light decreases changes with different optical frequencies.
  • Step S01 Using the measurement data, calculate the difference in phase offset value between the scattered light vector of the optical frequency component selected as the reference among the multiplexed optical frequency components and the scattered light vector of the other optical frequency components. Using the difference in the calculated phase offset values, a rotation angle for correcting the phase offset of each optical frequency is calculated.
  • Step S02 After rotating the scattered light vector at each time by the rotation angle, the rotated vectors of different optical frequencies are averaged to calculate a frequency average vector. An angular change in the frequency average vector is calculated using the frequency average vector at different times.
  • Step S03 Using the angle change, calculate the phase difference between two points separated by the gauge length on the optical fiber, perform phase connection processing, etc., and calculate the vibration waveform generated in the section between the two points. do.
  • step S01 the details of step S01 are as follows. However, assuming that the repetition interval of pulses of the same optical frequency is t, that is, the sampling interval of vibration is t, and that measurement data of M points from measurement time 0 to (M-1)t (M is a natural number) is used, each The scattered light vector of the reference optical frequency (f 1 ) at time mt (m is an integer from 0 to (M-1)) is defined as r 1 (mt,z) (distance from the input end as it depends on the fiber point). z), and the scattered light vector of each other optical frequency f i (i is an integer from 2 to N, and N represents the number of multiplexed optical frequencies) is defined as r i (mt, z).
  • M is set to include all points during the measurement time or to include the first few points.
  • Step S01-1 Calculate the angle ⁇ 1 (mt, z) of the scattered light vector r 1 (mt, z) of the reference optical frequency at each time and each fiber point. If the scattered light vector r 1 (mt, z) is considered to be a complex vector on a complex plane, that is, a complex number, the angle ⁇ 1 (mt, z) can be calculated as arg[r 1 (mt, z)].
  • arg is an operator that gives the argument of a complex number.
  • Step S01-2 Rotate the scattered light vector r i (mt, z) of each optical frequency at each time and each fiber point by an angle ⁇ 1 (mt, z) and set it as r i_rot (mz, t). If the scattered light vector r i (mt, z) is a complex vector on a complex plane, the rotated vector r i_rot (mz, t) is exp[- ⁇ 1 (mt, z)]*r i (mt , z).
  • Step S01-3 The vector r i_rot (mt, z) after rotating the scattered light vector of each optical frequency at each fiber point is calculated for M points where m, which is the measurement data to be used, is from 0 to (M-1). Averaging and calculating a time-averaged vector r i_avet (z). The value of +arg[r i_avet (z)] is the difference in phase offset value. In actual calculations, a vector sum may be used instead of a vector average for M points. No matter which method is used, the same result will be obtained in subsequent processing.
  • Step S01-4 The rotation angle ⁇ i (z) when rotating the scattered light vector of each optical frequency at each fiber point in step S02 is calculated by -arg[r i_avet ( z)]. By rotating by -arg[r i_avet (z)], it is possible to correct the difference in phase offset.
  • step S01-4 Since the rotation angle ⁇ i (z) given in step S01-4 is time-averaged by the steps up to step S01-3, it is a value in which noise is efficiently reduced, and the rotation angle ⁇ i
  • the S/N ratio of the final vibration waveform can be improved and waveform distortion can be suppressed.
  • step S01 of the prior art the rotation angle of each optical frequency is calculated using the measurement data of M points from measurement time 0 to (M-1)t (M is a natural number), and the calculated rotation angle is used to calculate the rotation angle of each optical frequency.
  • step S02 frequency averaging is performed on the measurement data at M points from measurement time 0 to (M-1)t (M is a natural number) and the measurement data at other times.
  • the optimum value of the rotation angle calculated by the method of step S01 of the prior art is determined based on temporal changes in optical characteristics such as the oscillation frequency of the laser, temperature changes in the optical fiber to be measured itself, and optical fibers to be measured with large dynamic strains. It fluctuates over time depending on the voltage applied to it. Therefore, in long-term distributed vibration measurement, it is necessary to continuously update the optimum value of the rotation angle and perform the processing from step S02 onwards.
  • Non-Patent Document 3 and Patent Document 1 there is no method of performing the processes from step S02 onwards while continuously updating the optimum value of the rotation angle. Also, when frequency averaging is performed using the calculated rotation angle for the measurement data used to calculate the rotation angle, streaming processing of the data is not possible, and the measurement data is stored in the computer's memory until the frequency averaging is completed. There is a need to hold it in place.
  • the present disclosure aims to enable streaming processing of data and prevent deterioration of the S/N ratio and vibration measurement accuracy without increasing the memory size requirements of the computer even during long-time measurements.
  • the present disclosure relates to a signal processing method that performs processing from step S02 onwards while continuously updating the optimum value of the rotation angle.
  • the present disclosure enables data to be streamed and eliminates the need to hold measurement data in a computer's memory until frequency averaging is completed.
  • the measurement system of the present disclosure includes: a measurement device that measures scattered light obtained by scattering multiple light pulses with different optical frequencies through an optical fiber; Obtain measurement data of scattered light measured by the measuring device, calculate a scattered light vector r i of scattered light at an arbitrary point of the optical fiber, and use the calculated scattered light vector r i to a signal processing device of the present disclosure that calculates a phase change at an arbitrary point; Equipped with
  • the signal processing device of the present disclosure executes the signal processing method of the present disclosure.
  • the signal processing method of the present disclosure is a signal processing method that performs signal processing of measurement data of scattered light in which N optical pulses with different optical frequencies are scattered by an optical fiber, Calculating the scattered light vector r i of scattered light at each optical frequency in a predetermined time range in which any point on the optical fiber can be measured, Average the scattered light vector r i of each optical frequency over the time range, Using the time average vector r i_avet obtained by the averaging, calculate the rotation angle ⁇ i for correcting the phase of each optical frequency, By rotating the phase of the scattered light vector r i of each optical frequency by a rotation angle ⁇ i , the phase of the scattered light vector r i of each optical frequency is corrected.
  • the time range may be the time range of any one block when the measurement data acquisition time is divided into a plurality of blocks.
  • the signal processing device of the present disclosure corrects the phase of the scattered light vector r i of each optical frequency for each block.
  • the signal processing device of the present disclosure may correct the phase of the scattered light vector r i using the rotation angle ⁇ i calculated in the previous block.
  • the signal processing device of the present disclosure includes: an AD conversion element that individually converts the in-phase components of the scattered light into digital signals; an AD conversion functional element that individually converts the orthogonal components of the scattered light into digital signals; an optical frequency separation unit that separates an in-phase component I digital and a quadrature component Q digital of the digital signal output from the AD conversion functional element into a signal I i and a signal Q i of each optical frequency; A scattered light vector r i is calculated using the signal I i and signal Q i output from the optical frequency separation section, the calculated scattered light vector r i is averaged in the time range, and the time obtained by the average is calculated.
  • a rotation angle calculation unit that calculates the rotation angle ⁇ i using the average vector r i_avet ; a phase correction unit that corrects the phase of the scattered light vector r i calculated by the rotation angle calculation unit using the rotation angle ⁇ i calculated by the rotation angle calculation unit; Equipped with Processing in the AD conversion functional element, the optical frequency separation section, and the rotation angle calculation section may be able to be executed independently and in parallel.
  • the signal processing device of the present disclosure includes: The frequency of the vector R i after rotating the rotation angle ⁇ i is averaged, Calculate the angle ⁇ avef of the frequency average vector R avef after frequency averaging, The vibration waveform may be calculated using the angle ⁇ avef .
  • the present disclosure allows data to be streamed while reducing the memory size requirements of a computer, and is theoretically applicable to real-time vibration waveform monitoring. Therefore, the present disclosure can perform streaming processing of data, and can prevent deterioration of the S/N ratio and vibration measurement accuracy without increasing the memory size requirements of the computer, even during long-time measurements.
  • 1 shows a configuration example of a measuring device according to the present disclosure. 1 shows an example of a light pulse of the present disclosure. Detailed processing of the signal processing unit 17d will be shown for the K-th block.
  • FIG. 1 shows a configuration example of a measurement system according to the present disclosure.
  • the measurement system according to the present disclosure includes a measurement device 31 that performs coherent detection using a 90-degree optical hybrid 7 in a receiving system, and a signal processing device 17 that executes the signal processing method of the present disclosure.
  • the measuring device 31 measures the scattered light from the optical fiber 6 to be measured using a plurality of optical pulses 4 having different optical frequencies.
  • a single wavelength continuous light having an optical frequency f0 is emitted from a CW light source 1, and is split into a reference light and a probe light by a coupler 2.
  • the probe light is shaped by the optical modulator 3 into an optical pulse 4 as shown in FIG.
  • ⁇ f i is a modulation frequency given from the modulator 3, and is a shift amount from the laser frequency f 0 of the CW light source 1.
  • the optical frequency f i is selected so that the intensity of the scattered light at each time and point is sufficiently far apart to the extent that different i's can be considered uncorrelated.
  • the pulse width W corresponds to the spatial resolution.
  • the optical modulator 3 is any device that can generate an optical pulse 4 with a pulse width W at a vibration sampling interval t, and may be composed of a plurality of units.
  • modulation using an LN modulator-based modulator such as an SSB (single side band) modulator may be used, a frequency variable AO modulator, etc. may be used, and the extinction ratio in pulsing may be
  • LN modulator-based modulator such as an SSB (single side band) modulator
  • a frequency variable AO modulator etc.
  • the extinction ratio in pulsing may be
  • further intensity modulation may be performed using an SOA (semiconductor optical amplifier) or the like.
  • the optical pulse 4 is input into the optical fiber 6 to be measured via the circulator 5.
  • the light scattered at each point in the longitudinal direction of the optical fiber 6 returns to the circulator 5 as backscattered light and enters one input of the 90-degree optical hybrid 7.
  • the reference light branched by the coupler 2 is input to the other input of the 90-degree optical hybrid 7.
  • the 90-degree optical hybrid 7 may have any internal configuration as long as it has the functions of a 90-degree optical hybrid. Two outputs among the four outputs of the 90-degree optical hybrid 7 are detected by the balance detector 13, and an electric signal 15 which is an analog in-phase component I analog is output. The remaining two outputs of the 90-degree optical hybrid are detected by a balance detector 14, and an electrical signal 16, which is an analog orthogonal component Q analog, is output. As a result, the light intensity of the in-phase component and the light intensity of the orthogonal component are individually measured for the scattered light obtained by scattering a plurality of light pulses having different optical frequencies by the optical fiber.
  • the electrical signal 15 and the electrical signal 16 are sent to a signal processing device 17 that includes an AD conversion functional element 17a and an AD conversion functional element 17b that can sample the frequency band of the signal without aliasing.
  • the signal processing device 17 performs signal processing on the digitized digital signals of the in-phase component I digital and the orthogonal component Q digital output from the AD conversion functional element 17a and the AD conversion functional element 17b.
  • the signal is separated into a signal I i and a signal Q i .
  • the signal processing unit 17c functions as an "optical frequency separation unit", and the in-phase component signal I i and the orthogonal component signal Q obtained when optical pulses of each optical frequency f 0 +f i component are input independently. i is separated by performing signal processing on I digital , which is a superposition of in-phase components for all i, and Q digital , which is a superposition of orthogonal components for all i.
  • I digital and Q digital may be passed through a digital bandpass filter whose center frequency is f 0 + ⁇ f i and whose passband is 2/W.
  • the filter characteristics of the digital bandpass filter are set to a specification that allows streaming processing of measurement data, taking into consideration the memory size of the computer used, the measurement data will be streamed up to the signal processing unit 17c.
  • the number of filter taps is too large or the number of effective decimal points of the filter coefficients is too large, the memory and calculation speed of the computer used may not be suitable for streaming processing. Therefore, the number of taps of the filter and the number of effective decimal points of the filter coefficients are set according to the memory and calculation speed of the computer used to enable streaming processing.
  • each signal processing section is as follows. Let t be the incident interval of the optical pulse of the i-th optical frequency, that is, the vibration sampling interval, and assume that vibration data is acquired at time kt using an integer k. The interval between t1 and t1' in FIG. 2 is t. Also, let z be the distance from the incident end. By performing phase calculation for each distance z, the phase at each point in the optical fiber 6 to be measured can be calculated.
  • Signal processing section 17d functions as a "rotation angle calculation section". Specifically, using the signal I i (kt, z) and the signal Q i (kt, z), the i-th optical frequency is determined in a predetermined time range in which the point of the optical fiber 6 to be measured can be measured. The scattered light vector r i (kt,z) of the scattered light at is calculated. The calculated scattered light vector r i (kt,z) data is continuously streamed to the signal processing unit 17e.
  • the rotation angle ⁇ i ( z) of each optical frequency fi is calculated from the scattered light vector r i (kt, z).
  • the measurement data is divided into blocks of predetermined time ranges on the time axis, and the rotation angle is calculated and updated for each block of the time range.
  • the rotation angle ⁇ i (z) calculated and updated for each block is delivered to the signal processing unit 17e.
  • one block is divided so that a predetermined number of optical pulses having the same frequency are included.
  • the predetermined time range is the predetermined number M of time points, and the rotation angle is calculated and updated every measurement time Mt.
  • all measurement data is treated as one block, or the rotation angle is calculated using the measurement data of the first few points, and the rotation angle calculated for all the data is used.
  • the rotation angle is updated for each block as needed.
  • K is used as a number symbol for distinguishing blocks.
  • the K-th block is measurement data where k is from M(K-1) to MK-1. Since the rotation angle is calculated for each block, the value changes depending on the block K, and the rotation angle is expressed as ⁇ i (z, K).
  • the rotation angle calculation method for each block is performed in the same manner as the procedure S01 of the prior art.
  • Processing S17d-0 As input to the signal processing unit 17d, the signal I i (kt, z) and the signal Q i (kt, z) are streamed from the signal processing unit 17c in ascending order of k.
  • Processing S17d-6 When the process S17d-5 is completed, the rotation angle ⁇ i (z) is delivered to the signal processing unit 17e.
  • Signal processing section 17e functions as a "phase correction section” and executes the above-mentioned step S02.
  • the signal processing unit 17e corrects the phase of the scattered light vector r i using the rotation angle ⁇ i (z) calculated by the signal processing unit 17d. For example, the signal processing unit 17e rotates the scattered light vector r i of each optical frequency f i included in the K block by a rotation angle ⁇ i (z), averages the rotation angle ⁇ i (z), and calculates the frequency average vector R avef .
  • the rotation angle ⁇ i is updated for each block of time Mt, so in this disclosure, the rotation angle used when calculating the frequency average of the scattered light vector belonging to a certain block K is
  • the processing unit 17d uses the rotation angle calculated using the previous block K-1. That is, the detailed processing of the signal processing unit 17e for the Kth block is as follows.
  • Process S17e-1 Rotate the scattered light vector r i (kt, z) of each optical frequency at time kt by the rotation angle ⁇ i (z, K-1) calculated in the K-1st block, and after rotation Calculate the vector R i (kz,t).
  • R i (kz,t) is expressed by the following formula. (Number 6)
  • R i (kz, t) exp(j ⁇ i (z, K-1)) ⁇ r i (kt, z) (6)
  • Process S17e-2 Calculate a vector obtained by frequency-averaging the rotated vector R i (kz, t). Note that since the final result is the same whether it is vector averaging or vector synthesis (simple vector sum), vector synthesis is performed in the actual calculation procedure to obtain the frequency average vector R avef (kz, t). R avef (kz, t) is expressed by the following formula.
  • Processing S17e-3 Calculate the angle ⁇ avef (kz, t) of the frequency average vector R avef (kz, t) and deliver it to the signal processing unit 17f.
  • ⁇ avef (kz, t) is expressed by the following equation.
  • the first Mt measurement time is treated as a pre-measurement, and no processing is performed after the signal processing section 17e.
  • Signal processing unit 17f executes the above-mentioned step S03.
  • the signal processing unit 17f calculates the phase change at an arbitrary point of the optical fiber 6 to be measured using the angle ⁇ avef (kz,t) of the frequency average vector R ave obtained by the signal processing unit 17e.
  • the signal processing unit 17f calculates a phase difference between two points separated by a gauge length, performs phase connection processing, etc., and calculates a vibration waveform.
  • the specific calculation method is the same as the conventional method.
  • the feature of the present disclosure lies in the signal processing units 17d and 17e.
  • the signal processing unit 17d divides the signals I i and Q i obtained by measurement into blocks of time points M, and calculates and updates the rotation angle ⁇ i for each block of measurement time Mt. Responds to temporal changes in the optical properties such as the oscillation frequency of the provided laser, temperature changes in the optical fiber 6 to be measured, application of large dynamic strain to the optical fiber 6 to be measured, etc. .
  • the measurement time Mt per block is determined based on temporal changes in optical characteristics such as the oscillation frequency of the laser provided in the CW light source 1, temperature changes in the optical fiber 6 itself, and large dynamic strain on the optical fiber 6 under test.
  • the rotation angle can be updated with sufficient detail by keeping it small relative to the time scale of the application.
  • the rotation angle ⁇ i used when the signal processing unit 17e calculates the frequency average of the scattered light vector r i belonging to the block K is the rotation angle ⁇ i of the immediately preceding block K ⁇ calculated by the signal processing unit 17d. There is no problem even if a rotation angle ⁇ i of 1 is used. Further, the signal processing unit 17d outputs the rotation angle ⁇ i of the immediately preceding block K-1 to the signal processing unit 17e, so that the scattered light vector r i_avet (z) is updated by the signal processing unit 17d. There is no need to hold the data of the vector r i (kt,z) in the memory of the signal processing device 17.
  • the present disclosure provides streaming processing in which the processing in the AD conversion functional elements 17a and 17b, the signal processing unit 17c functioning as an optical frequency separation unit, and the signal processing unit 17d functioning as a rotation angle calculation unit is performed independently and in parallel. is possible.
  • the present disclosure uses the rotation angle itself calculated by the signal processing unit 17d in the block K as the rotation angle ⁇ i used when the signal processing unit 17e calculates the frequency average of the scattered light vector r i belonging to the block K.
  • the memory of the computer can be reduced compared to the conventional technology used in which the signal processing unit 17d has to hold all the data of the scattered light vectors belonging to the block K in the memory of the computer.
  • the length of the measurement time Mt per block is also related to the calculation accuracy of the rotation angle ⁇ i .
  • the measurement time Mt per block is set to be long enough to calculate the rotation angle ⁇ i with sufficient accuracy, and the measurement time is set to be long enough to calculate the rotation angle ⁇ i with sufficient accuracy.
  • Non-Patent Document 3 it has been experimentally found that sufficient accuracy of the rotation angle ⁇ i can be obtained by setting M to about 100 (Non-Patent Document 3); however, even if the optical pulse transmission period t is 1 ms, Mt is approximately 100 ms, which is due to temporal changes in optical characteristics such as the oscillation frequency of the laser provided in the CW light source 1, temperature changes in the optical fiber 6 itself, and large dynamic strain applied to the optical fiber 6 under test. If the time scale is 1 s or more, it is considered to be sufficiently fine.
  • the device of the present disclosure can also be realized by a computer and a program, and the program can be recorded on a recording medium or provided through a network.
  • the program of the present disclosure is a program for realizing a computer as each functional unit included in the device of the present disclosure, and is a program for causing the computer to execute each step of a method executed by the device of the present disclosure. .
  • CW light source 2 Coupler 3: Optical modulator 4: Optical pulse 5: Circulator 6: Optical fiber to be measured 7: 90 degree optical hybrid 13, 14: Balance detectors 15, 16: Electrical signal 17: Signal processing device 17a , 17b: AD conversion functional element 17c, 17d: Signal processing section 31: Measuring device

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Optical Transform (AREA)

Abstract

本開示は、データをストリーミング処理することができ、長時間の測定であっても、計算機のメモリサイズ要件を高めることなく、ＳＮ比や振動測定精度の劣化を防ぐことを目的とする。　本開示の信号処理装置及び信号処理方法は、光周波数の異なる複数の光パルスが光ファイバで散乱された散乱光の測定データの信号処理を行う信号処理方法及び信号処理方法であって、前記光ファイバの任意の地点を測定可能な予め定められた時間範囲における、各光周波数での散乱光の散乱光ベクトルｒｉを計算し、各光周波数の散乱光ベクトルｒｉを前記時間範囲で平均し、前記平均によって得られた時間平均ベクトルｒｉ＿ａｖｅｔを用いて、各光周波数の位相を補正するための回転角度αｉを計算し、各光周波数の散乱光ベクトルｒｉの位相を回転角度αｉ回転させることで、各光周波数の散乱光ベクトルｒｉの位相を補正する。

Description

位相ＯＴＤＲにおける信号処理方法
　本開示は、被測定光ファイバの各地点からの散乱光の位相を測定する位相ＯＴＤＲに関する。
　光ファイバに加わった物理的な振動を、光ファイバ長手方向に分布的に計測する手段として、被測定光ファイバにパルス試験光を入射し、レイリー散乱による後方散乱光を検出するＤＡＳ（Ｄｉｓｔｒｉｂｕｔｅｄ　Ａｃｏｕｓｔｉｃ　Ｓｅｎｓｉｎｇ）と呼ばれる手法が知られている（非特許文献１）。
　ＤＡＳの手法の一つとして、被測定光ファイバの各地点からの散乱光の位相を測定し、位相の時間変化を観測する位相ＯＴＤＲ（Ｏｐｔｉｃａｌ　Ｔｉｍｅ　Ｄｏｍａｉｎ　Ｒｅｆｌｅｃｔｏｍｅｔｒｙ）がある。位相ＯＴＤＲはＤＡＳ－Ｐ（ＤＡＳ－ｐｈａｓｅ）とも呼ばれる。ＤＡＳ－Ｐは、振動による光ファイバの光路長変化に対して位相が線形に変化し、その変化率も光ファイバ長手方向の各地点でおよそ同一とみなせるため、振動の定量的な測定が可能となり、被測定光ファイバに加わった振動波形を忠実に再現することができる（例えば非特許文献２）。
　単一光周波数の光パルスを用いたＤＡＳ－Ｐでは、パルス幅内の散乱光同士の干渉効果により、散乱光強度が小さい地点が発生し、感度の劣化が生じる。感度の劣化を防ぐ方法として、光周波数多重を行い、異なる光周波数で信号を平均する方法がある（非特許文献３、特許文献１）。この方法は、異なる光周波数では散乱光強度が小さくなる地点が変化することを利用している。
　非特許文献３及び特許文献１では効率的な平均方法として具体的に以下の手順を行う。
　手順Ｓ０１：測定データを使用して、多重した光周波数成分の内基準として選んだ光周波数成分の散乱光ベクトルと、その他光周波数成分の散乱光ベクトルとの、位相オフセット値の違いを計算する。計算した位相オフセット値の違いを用いて、各光周波数の位相オフセットを補正するための回転角度を計算する。
　手順Ｓ０２：各時刻における散乱光ベクトルを前記回転角度だけ回転させた上で、異なる光周波数の回転後のベクトルを平均し、周波数平均ベクトルを計算する。異なる時刻の周波数平均ベクトルを用いて、周波数平均ベクトルの角度変化を計算する。
　手順Ｓ０３：前記角度変化を用いて、光ファイバ上でゲージ長だけ離れた２地点の位相の差分を計算し、位相接続処理などを行い、前記２地点の間の区間に生じた振動波形を計算する。
　特に手順Ｓ０１の詳細は以下の具体的な手順を行う。ただし、同じ光周波数のパルスの繰り返し間隔をｔ、つまり、振動のサンプリング間隔をｔとして、測定時刻０から（Ｍ－１）ｔ（Ｍは自然数）のＭ点の測定データを使用するとして、各時刻ｍｔ（ｍは０から（Ｍ－１）の整数）における前記基準光周波数（ｆ１とおく）の散乱光ベクトルをｒ１（ｍｔ，ｚ）（ファイバ地点に依存するため入射端からの距離ｚに依存）として、その他各光周波数ｆｉ（ｉは２からＮまでの整数でＮは光周波数の多重数を表す）の散乱光ベクトルをｒｉ（ｍｔ，ｚ）とする。ここで、従来では、Ｍは測定時間全点を含むか、最初の数点を含むように設定する。
　手順Ｓ０１－１：各時刻・各ファイバ地点における基準光周波数の散乱光ベクトルｒ１（ｍｔ，ｚ）の角度θ１（ｍｔ，ｚ）を計算する。散乱光ベクトルｒ１（ｍｔ，ｚ）を、複素平面上の複素ベクトル、つまり複素数だと考えれば、角度θ１（ｍｔ，ｚ）はａｒｇ［ｒ１（ｍｔ，ｚ）］で計算できる。ここでａｒｇは複素数の偏角を与える演算子である。
　手順Ｓ０１－２：各時刻・各ファイバ地点における各光周波数の散乱光ベクトルｒｉ（ｍｔ，ｚ）を角度－θ１（ｍｔ，ｚ）だけ回転させ、ｒｉ＿ｒｏｔ（ｍｚ，ｔ）とする。散乱光ベクトルｒｉ（ｍｔ，ｚ）を、複素平面上の複素ベクトルとすれば、回転後のベクトルｒｉ＿ｒｏｔ（ｍｚ，ｔ）はｅｘｐ［－θ１（ｍｔ，ｚ）］＊ｒｉ（ｍｔ，ｚ）と計算できる。
　手順Ｓ０１－３：各ファイバ地点における各光周波数の散乱光ベクトルを回転した後のベクトルｒｉ＿ｒｏｔ（ｍｔ，ｚ）を、使用する測定データであるｍが０から（Ｍ－１）のＭ点分平均し、時間平均ベクトルｒｉ＿ａｖｅｔ（ｚ）を計算する。＋ａｒｇ［ｒｉ＿ａｖｅｔ（ｚ）］の値が位相オフセット値の違いとなる。実際の計算においてはＭ点分のベクトル平均ではなくてベクトル合計でも良い。どちらを使用しても引き続く処理で同一の結果が得られる。
　手順Ｓ０１－４：各ファイバ地点における各光周波数の散乱光ベクトルを前記手順Ｓ０２で回転させる際の回転角度αｉ（ｚ）を、前記時間平均ベクトルの偏角を用いて－ａｒｇ［ｒｉ＿ａｖｅｔ（ｚ）］と計算する。－ａｒｇ［ｒｉ＿ａｖｅｔ（ｚ）］だけ回転させることで、位相オフセットの違いを補正することが可能である。
　手順Ｓ０１－４で与えられる回転角度αｉ（ｚ）は手順Ｓ０１－３までの手順により時間平均を行っていることから雑音が効率的に低減された値となっており、その回転角度αｉ（ｚ）を用いて手順Ｓ０２で光周波数平均ベクトルを計算することにより、最終的な振動波形のＳＮ比が向上し、波形歪みを抑えることができる。
　前記従来技術の手順Ｓ０１では測定時刻０から（Ｍ－１）ｔ（Ｍは自然数）のＭ点の測定データを使用して各光周波数の回転角度を計算し、計算した回転角度を使用して手順Ｓ０２で測定時刻０から（Ｍ－１）ｔ（Ｍは自然数）のＭ点の測定データおよびその他の時刻の測定データに対して周波数平均を行う。
特開２０２０－１６９９０４号公報
Ａｌｉ．　Ｍａｓｏｕｄｉ，　Ｔ．　Ｐ．　Ｎｅｗｓｏｎ，　"Ｃｏｎｔｒｉｂｕｔｅｄ　Ｒｅｖｉｅｗ：　Ｄｉｓｔｒｉｂｕｔｅｄ　ｏｐｔｉｃａｌ　ｆｉｂｒｅ　ｄｙｎａｍｉｃ　ｓｔｒａｉｎ　ｓｅｎｓｉｎｇ．"　Ｒｅｖｉｅｗ　ｏｆ　Ｓｃｉｅｎｔｉｆｉｃ　Ｉｎｓｔｒｕｍｅｎｔｓ，　ｖｏｌ．８７，　ｐｐ０１１５０１　（２０１６） 西口憲一，李哲賢，グジクアーター，横山光徳，増田欣増「光ファイバによる分布型音波センサの試作とその信号処理」信学技報，１１５（２０２），　ｐｐ２９－３４　（２０１５） Ｙｏｓｈｉｆｕｍｉ　Ｗａｋｉｓａｋａ，　Ｄａｉｓｕｋｅ　Ｉｉｄａ，　Ｈｉｒｏｙｕｋｉ　Ｏｓｈｉｄａ，　ａｎｄ　Ｎａｚｕｋｉ　Ｈｏｎｄａ，　"Ｆａｄｉｎｇ　Ｓｕｐｐｒｅｓｓｉｏｎ　ｏｆ　Φ－ＯＴＤＲ　Ｗｉｔｈ　ｔｈｅ　Ｎｅｗ　Ｓｉｇｎａｌ　Ｐｒｏｃｅｓｓｉｎｇ　Ｍｅｔｈｏｄｏｌｏｇｙ　ｏｆ　Ｃｏｍｐｌｅｘ　Ｖｅｃｔｏｒｓ　Ａｃｒｏｓｓ　Ｔｉｍｅ　ａｎｄ　Ｆｒｅｑｕｅｎｃｙ　Ｄｏｍａｉｎｓ，"　Ｊ．　Ｌｉｇｈｔｗａｖｅ　Ｔｅｃｈｎｏｌ．　３９，　４２７９－４２９３　（２０２１）
　前記従来技術の手順Ｓ０１の方法で計算する回転角度の最適値は、レーザの発振周波数などの光学特性の時間的な変化や被測定光ファイバ自体の温度変化、大きな動的歪みの被測定光ファイバへの印加などにより時間的に変動する。したがって、長時間での分布振動測定においては、回転角度の最適値を連続的に更新し、手順Ｓ０２以降の処理を行う必要がある。
　しかし、非特許文献３、特許文献１に記載のような従来技術においては、回転角度の最適値を連続的に更新しながら手順Ｓ０２以降の処理を行う方法は示されていない。また、回転角度の計算に使用した測定データに対して、計算した回転角度で周波数平均を行う場合には、データをストリーミング処理することができず、周波数平均が完了するまで測定データを計算機のメモリに保持する必要などが生じる。
　本開示は、データをストリーミング処理することができ、長時間の測定であっても、計算機のメモリサイズ要件を高めることなく、ＳＮ比や振動測定精度の劣化を防ぐことを目的とする。
　本開示は、回転角度の最適値を連続的に更新しながら手順Ｓ０２以降の処理を行う信号処理方法に関する。本開示は、データをストリーミング処理することを可能とし、周波数平均が完了するまでの測定データの計算機のメモリ上での保持などを不要とする。
　具体的には、本開示の測定システムは、
　光周波数の異なる複数の光パルスが光ファイバで散乱された散乱光を測定する測定装置と、
　前記測定装置で測定された散乱光の測定データを取得し、前記光ファイバの任意の地点での散乱光の散乱光ベクトルｒｉを計算し、計算した散乱光ベクトルｒｉを用いて前記光ファイバの任意の地点での位相変化を計算する、本開示の信号処理装置と、
　を備える。
　本開示の信号処理装置は、本開示の信号処理方法を実行する。本開示の信号処理方法は、光周波数の異なるＮ個の光パルスが光ファイバで散乱された散乱光の測定データの信号処理を行う信号処理方法であって、
　前記光ファイバの任意の地点を測定可能な予め定められた時間範囲における、各光周波数での散乱光の散乱光ベクトルｒｉを計算し、
　各光周波数の散乱光ベクトルｒｉを前記時間範囲で平均し、
　前記平均によって得られた時間平均ベクトルｒｉ＿ａｖｅｔを用いて、各光周波数の位相を補正するための回転角度αｉを計算し、
　各光周波数の散乱光ベクトルｒｉの位相を回転角度αｉ回転させることで、各光周波数の散乱光ベクトルｒｉの位相を補正する。
　前記時間範囲が、測定データの取得時間を複数のブロックに区切ったときの、任意の１つのブロックの時間範囲であってもよい。この場合、本開示の信号処理装置は、ブロックごとに各光周波数の散乱光ベクトルｒｉの位相を補正する。ここで、本開示の信号処理装置は、１つ前のブロックで計算された前記回転角度αｉを用いて、前記散乱光ベクトルｒｉの位相を補正してもよい。
　本開示の信号処理装置は、
　散乱光の同相成分を個別にデジタル信号に変換するＡＤ変換素子と、
　散乱光の直交成分を個別にデジタル信号に変換するＡＤ変換機能素子と、
　前記ＡＤ変換機能素子から出力されたデジタル信号の同相成分Ｉｄｉｇｉｔａｌと直交成分Ｑｄｉｇｉｔａｌを、各光周波数の信号Ｉｉと信号Ｑｉに分離する光周波数分離部と、
　前記光周波数分離部から出力された信号Ｉｉと信号Ｑｉを用いて散乱光ベクトルｒｉを計算し、計算した散乱光ベクトルｒｉを前記時間範囲で平均し、前記平均によって得られた時間平均ベクトルｒｉ＿ａｖｅｔを用いて前記回転角度αｉを計算する回転角度計算部と、
　前記回転角度計算部で計算された前記回転角度αｉを用いて、前記回転角度計算部で計算された前記散乱光ベクトルｒｉの位相を補正する位相補正部と、
　を備え、
　前記ＡＤ変換機能素子、前記光周波数分離部及び前記回転角度計算部における処理が、独立かつ並列に実行可能あってもよい。
　本開示の信号処理装置は、
　回転角度αｉを回転後のベクトルＲｉを周波数平均し、
　周波数平均後の周波数平均ベクトルＲａｖｅｆの角度θａｖｅｆを計算し、
　角度θａｖｅｆを用いて振動波形を計算してもよい。
　なお、上記各開示は、可能な限り組み合わせることができる。
　本開示を用いることで、長時間の分布振動測定においても、ＳＮ比の劣化や振動測定精度の劣化を防ぎながら、振動波形を計測することが可能となる。また本開示は、計算機のメモリサイズ要件を低減しつつ、データをストリーミング処理することが可能で、リアルタイムな振動波形モニタリングへ適用することが原理的に可能である。したがって、本開示は、データをストリーミング処理することができ、長時間の測定であっても、計算機のメモリサイズ要件を高めることなく、ＳＮ比や振動測定精度の劣化を防ぐことができる。
本開示に係る測定装置の構成例を示す。 本開示の光パルスの一例を示す。 信号処理部１７ｄの詳細処理をＫブロック目について示す。
　以下、本開示の実施形態について、図面を参照しながら詳細に説明する。なお、本開示は、以下に示す実施形態に限定されるものではない。これらの実施の例は例示に過ぎず、本開示は当業者の知識に基づいて種々の変更、改良を施した形態で実施することができる。なお、本明細書及び図面において符号が同じ構成要素は、相互に同一のものを示すものとする。
　図１に、本開示に係る測定システムの構成例を示す。本開示に係る測定システムは、受信系に９０度光ハイブリッド７を用いてコヒーレント検波を行う測定装置３１と、本開示の信号処理方法を実行する信号処理装置１７とを備える。
　測定装置３１は、光周波数の異なる複数の光パルス４を用いて、被測定光ファイバ６からの散乱光を測定する。ＣＷ光源１から光周波数がｆ０の単一波長の連続光が射出され、カプラ２により参照光とプローブ光に分岐される。プローブ光は、光変調器３によって、図２に示すような光パルス４に整形される。光パルス４は、ｉ番目の光周波数ｆｉがｆｉ＝ｆ０＋Δｆｉ（ｉは整数）に設定されかつパルス幅Ｗが光ファイバ長手方向での測定の空間分解能に対応する値に設定された光パルスが、ｉ＝１，２，…，Ｎ（Ｎは整数で光周波数の多重数を表す）だけ並んだ構成となっている。Δｆｉは変調器３から与える変調周波数であり、ＣＷ光源１のレーザ周波数ｆ０からのシフト量である。光周波数ｆｉは、各時刻・各地点における散乱光の強度が、異なるｉ同士で無相関とみなせる程度まで十分に離れているように選択をする。パルス幅Ｗが空間分解能に対応する。
　光変調器３は、パルス幅Ｗの光パルス４を振動のサンプリング間隔ｔで生成できる任意の装置であり、複数台で構成されている場合もある。例えば、ＳＳＢ（ｓｉｎｇｌｅ　ｓｉｄｅ　ｂａｎｄ）変調器などのＬＮ変調器ベースの変調器を利用した変調を用いてもよいし、周波数可変なＡＯ変調器などを用いても良いし、パルス化における消光比を大きくするためにさらにＳＯＡ（半導体光増幅器）などによる強度変調を行っても良い。
　光パルス４は、サーキュレータ５を介して、被測定光ファイバ６に入射される。光ファイバ６の長手方向の各点で散乱された光が、後方散乱光としてサーキュレータ５に戻り、９０度光ハイブリッド７の一方のインプットに入射される。カプラ２により分岐された参照光は、９０度光ハイブリッド７のもう一方のインプットに入射される。
　９０度光ハイブリッド７の内部構成は、９０度光ハイブリッドの機能さえ備えていれば、なんでもよい。９０度光ハイブリッド７の４つのアウトプットの内、２つのアウトプットがバランス検出器１３によって検出され、アナログの同相成分Ｉａｎａｌｏｇである電気信号１５が出力される。９０度光ハイブリッドの残り２つのアウトプットがバランス検出器１４によって検出され、アナログの直交成分Ｑａｎａｌｏｇである電気信号１６が出力される。これにより、光周波数の異なる複数の光パルスが光ファイバで散乱された散乱光について、同相成分の光強度及び直交成分の光強度が個別に測定される。
　電気信号１５と電気信号１６は、信号の周波数帯域をエイリアシングなくサンプリングが可能なＡＤ変換機能素子１７ａとＡＤ変換機能素子１７ｂを備えた信号処理装置１７に送られる。信号処置装置１７では、ＡＤ変換機能素子１７ａとＡＤ変換機能素子１７ｂから出力されたデジタル化された同相成分Ｉｄｉｇｉｔａｌと直交成分Ｑｄｉｇｉｔａｌのデジタル信号に対して信号処理を行う。具体的には、信号処理部１７ｃは、同相成分Ｉｄｉｇｉｔａｌと直交成分Ｑｄｉｇｉｔａｌを、光パルス４を構成する各光周波数ｆｉ（ｉ＝１，２，…，Ｎ）の光パルスによる散乱光の信号Ｉｉと信号Ｑｉに分離する。
　つまり、信号処理部１７ｃは、「光周波数分離部」として機能し、各光周波数ｆ０＋ｆｉ成分の光パルスを単独で入射した場合に得られる同相成分の信号Ｉｉと直交成分の信号Ｑｉを、全てのｉに関する同相成分の重ね合わせとなっているＩｄｉｇｉｔａｌと、全てのｉに関する直交成分の重ね合わせとなっているＱｄｉｇｉｔａｌに対して信号処理を行うことで分離する。具体的な信号処理の方法は、例えばＩｄｉｇｉｔａｌとＱｄｉｇｉｔａｌを中心周波数がｆ０＋Δｆｉであり通過帯域が２／Ｗであるデジタルバンドパスフィルタに通すことなどが考えられる。例えばＦＩＲ（Ｆｉｎｉｔｅ　Ｉｍｐｕｌｓｅ　Ｒｅｓｐｏｎｓｅ）フィルタを用いる場合には、フィルタ特性に応じた遅延も生じるため、その遅延を補償するために、フィルタ後の信号を遅延した時間分だけ逆に前倒しするなどの処理もあわせて実施する。
　使用する計算機のメモリサイズなどを考慮して、デジタルバンドパスフィルタのフィルタ特性を測定データがストリーミング処理可能な仕様に設定しておけば、信号処理部１７ｃまでは測定データはストリーミング処理される。例えば、デジタルバンドパスフィルタとしてＦＩＲフィルタを用いる場合には、フィルタのタップ数が大きすぎたり、フィルタ係数の有効小数点数が多すぎたりすると、使用する計算機のメモリや計算速度がストリーミング処理には適さなくなるため、使用する計算機のメモリや計算速度に応じて、フィルタのタップ数やフィルタ係数の有効小数点数を設定し、ストリーミング処理が可能なようにする。
　信号処理部１７ｃによって取得された信号Ｉｉと信号Ｑｉを元に、信号処理部１７ｄと信号処理部１７ｅと信号処理部１７ｆが引き続き位相計算を行う。それぞれの信号処理部の役割を以下とする。ｉ番目の光周波数の光パルスの入射間隔、つまり振動のサンプリング間隔をｔとし、整数ｋを用いて時刻ｋｔで振動データを取得しているとする。図２のｔ１とｔ１’の間隔がｔである。また、入射端からの距離をｚとする。距離ｚごとに位相計算を行うことで、被測定光ファイバ６における各地点での位相を計算することができる。
　信号処理部１７ｄ：「回転角度計算部」として機能する。
　具体的には、信号Ｉｉ（ｋｔ，ｚ）と信号Ｑｉ（ｋｔ，ｚ）を用いて、被測定光ファイバ６の地点を測定可能な予め定められた時間範囲における、ｉ番目の光周波数での散乱光の散乱光ベクトルｒｉ（ｋｔ，ｚ）を計算する。計算した散乱光ベクトルｒｉ（ｋｔ，ｚ）データを信号処理部１７ｅに連続的にストリーミングする。
　また、散乱光ベクトルｒｉ（ｋｔ，ｚ）から各光周波数ｆｉの回転角度αｉ（ｚ）を計算する。ただし、従来方法と違い、測定データを時間軸上で予め定められた時間範囲のブロックに区切り、前記時間範囲のブロック単位で回転角度を計算・更新する。ブロックごとに計算・更新した回転角度αｉ（ｚ）を信号処理部１７ｅに引き渡す。ここで、本実施形態では、１つのブロックに、同じ周波数の光パルスが予め定められた個数含まれるように区切る。以下、予め定められた時間範囲が予め定められた時間ポイント数Ｍであり、測定時間Ｍｔごとに回転角度を計算・更新する例について説明する。
　尚、従来方法では、全ての測定データを１ブロックとして扱うか、あるいは、最初の数点の測定データを用いて回転角度を計算し、全てのデータに対して計算した回転角度を使用するが、本実施形態例では回転角度をブロックごとに随時更新していく。
　すなわち、仮にｋ＝０を測定データの１点目として、１点目から回転角度の計算を開始するとすれば、例えばｋが０からＭ－１の散乱光ベクトルｒｉ（ｋｔ，ｚ）の測定データが一つ目のブロックとなり、ｋがＭから２Ｍ－１が二つ目のブロックとなる。ブロックを区別する番号の記号としてＫを用いる。例えば上記の例では、ｋが０からＭ－１はＫ＝１のブロック、ｋがＭから２Ｍ－１がＫ＝２のブロックである。Ｋブロック目は、ｋがＭ（Ｋ－１）からＭＫ－１の測定データである。各ブロックで回転角度を計算するため、ブロックＫに応じて値は変化し、回転角度をαｉ（ｚ，Ｋ）と表記する。それぞれのブロックでの回転角度の計算方法は従来技術の手順Ｓ０１と同様に行う。
　信号処理部１７ｄの詳細処理をＫブロック目について図２に示す。
　処理Ｓ１７ｄ－０：信号処理部１７ｄへの入力としては１７ｃからｋが小さい順に信号Ｉｉ（ｋｔ，ｚ）と信号Ｑｉ（ｋｔ，ｚ）がストリーミングされてくる。
　処理Ｓ１７ｄ－１：時刻ｋｔの信号Ｉｉ（ｋｔ，ｚ）と信号Ｑｉ（ｋｔ，ｚ）から散乱光ベクトルｒｉ（ｋｔ，ｚ）を計算していく。例えば、虚数単位をｊとして、以下を計算する。
（数１）
ｒｉ（ｋｔ，ｚ）＝Ｉｉ（ｋｔ，ｚ）＋ｊＱｉ（ｋｔ，ｚ）　　（１）
　計算した散乱光ベクトルのデータを順次信号処理部１７ｅにストリーミングする。
　処理Ｓ１７ｄ－２：
　基準光周波数をｆ１とする。ただし基準光周波数は任意に選べるのでｆ１でなくても良い。時刻ｋｔの各ファイバ地点における基準光周波数の散乱光ベクトルｒ１（ｋｔ，ｚ）の角度θ１（ｋｔ，ｚ）を計算する。θ１（ｋｔ，ｚ）は、例えば、次式で表すことができる。
（数２）
θ１（ｋｔ，ｚ）＝ａｒｇ［ｒ１（ｋｔ，ｚ）］　　（２）
　処理Ｓ１７ｄ－３：
　時刻ｋｔの各ファイバ地点における各光周波数の散乱光ベクトルｒｉ（ｋｔ，ｚ）を角度－θ１（ｋｔ，ｚ）だけ回転させ、ｒｉ＿ｒｏｔ（ｋｚ，ｔ）とする。ｒｉ＿ｒｏｔ（ｋｚ，ｔ）は、例えば、次式で表すことができる。
（数３）
ｒｉ＿ｒｏｔ（ｋｚ，ｔ）＝ｅｘｐ［－θ１（ｋｔ，ｚ）］＊ｒｉ（ｋｔ，ｚ）　　（３）
　処理Ｓ１７ｄ－４：
　時刻ｋｔがＫブロック目の最初の時刻Ｍ（Ｋ－１）の場合には、時間平均ベクトルｒｉ＿ａｖｅｔ（ｚ，Ｋ）を新たに用意し、ｒｉ＿ａｖｅｔ（ｚ）＝ｒｉ＿ｒｏｔ（ｋｚ，ｔ）とする。時刻ｋｔがＫブロック目の２番目以降の場合には、時刻（ｋ－１）ｔまでの時間平均ベクトルに時刻ｋｔでのｒｉ＿ｒｏｔ（ｋｚ，ｔ）を加算して新たな時間平均ベクトルに更新する。時間平均ベクトルｒｉ＿ａｖｅｔ（ｚ，Ｋ）は、次式で表される。
（数４）
ｒｉ＿ａｖｅｔ（ｚ，Ｋ）＝ｒｉ＿ｒｏｔ（ｋｚ，ｔ）＋ｒｉ＿ａｖｅｔ（ｚ，Ｋ）　　（４）
　処理Ｓ１７ｄ－５：
　時刻ｋｔがＫブロック目最後の時刻ＭＫ－１の際の処理Ｓ１７ｄ－４まで終了した段階で得られたｒｉ＿ａｖｅｔ（ｚ，Ｋ）を使用して、Ｋブロック目を使用した各光周波数の回転角度αｉ（ｚ）を計算する。αｉ（ｚ）は次式で表される。
（数５）
αｉ（ｚ，Ｋ）＝－ａｒｇ［ｒｉ＿ａｖｅｔ（ｚ，Ｋ）］　　（５）
　これにより、ブロックに含まれる各光周波数の位相オフセット値の違いを得ることができる。
　処理Ｓ１７ｄ－６：
　処理Ｓ１７ｄ－５が終了した段階で、回転角度αｉ（ｚ）を信号処理部１７ｅに引き渡す。
　信号処理部１７ｅ：「位相補正部」として機能し、前述の手順Ｓ０２を実行する。
　信号処理部１７ｅは、信号処理部１７ｄで計算した回転角度αｉ（ｚ）を用いて、散乱光ベクトルｒｉの位相を補正する。例えば、信号処理部１７ｅは、Ｋブロックに含まれる各光周波数ｆｉの散乱光ベクトルｒｉを回転角度αｉ（ｚ）回転させた上で平均し、周波数平均ベクトルＲａｖｅｆを計算する。
　ただし、従来方法と違い、回転角度αｉは時間Ｍｔのブロックごとに更新されているため、本開示では、あるブロックＫに属する散乱光ベクトルの周波数平均を計算する際に用いる回転角度は、信号処理部１７ｄにおいて一つ前のブロックＫ－１を用いて計算した回転角度を用いる。すなわち、信号処理部１７ｅの詳細処理はＫブロック目について、以下となる。
　処理Ｓ１７ｅ－１：時刻ｋｔの各光周波数の散乱光ベクトルｒｉ（ｋｔ，ｚ）を、Ｋ－１番目のブロックで計算した回転角度αｉ（ｚ，Ｋ－１）だけ回転させ、回転後のベクトルＲｉ（ｋｚ，ｔ）を計算する。Ｒｉ（ｋｚ，ｔ）は、次式で表される。
（数６）
Ｒｉ（ｋｚ，ｔ）＝ｅｘｐ（ｊ・αｉ（ｚ，Ｋ－１））・ｒｉ（ｋｔ，ｚ）　　（６）
ベクトルＲｉ（ｋｚ，ｔ）を計算することで、各光周波数の散乱光ベクトルｒｉの位相オフセットの違いを補正することができる。
　処理Ｓ１７ｅ－２：回転後のベクトルＲｉ（ｋｚ，ｔ）を周波数平均したベクトルを計算する。尚、ベクトル平均でもベクトル合成（単純なベクトル和）でも最終的な結果は変わらないため、実際の計算手順ではベクトル合成して周波数平均ベクトルＲａｖｅｆ（ｋｚ，ｔ）とする。Ｒａｖｅｆ（ｋｚ，ｔ）は次式で表される。
Figure JPOXMLDOC01-appb-M000001
　処理Ｓ１７ｅ－３：周波数平均ベクトルＲａｖｅｆ（ｋｚ，ｔ）の角度θａｖｅｆ（ｋｚ，ｔ）を計算し、信号処理部１７ｆに引き渡す。θａｖｅｆ（ｋｚ，ｔ）は次式で表される。
Figure JPOXMLDOC01-appb-M000002
　尚、Ｋ＝１の最初のブロックにおいては、それ以前のブロックのデータがないため回転角度が未取得の状態であり、信号処理部１７ｅの手順を行うことができない。したがって、最初のＭｔの測定時間分はプレ測定として、信号処理部１７ｅ以降の処理は行わない。あるいは、最初のＭｔの測定時間分の散乱光ベクトルｒｉ（ｋｔ，ｚ）だけ別途計算機内で保存し、従来手法により振動波形まで計算することも可能である。
　信号処理部１７ｆ：前述の手順Ｓ０３を実行する。
　信号処理部１７ｆは、信号処理部１７ｅで得られた周波数平均ベクトルＲａｖｅｆの角度θａｖｅｆ（ｋｚ，ｔ）を用いて、被測定光ファイバ６の任意の地点での位相変化を計算する。例えば、信号処理部１７ｆは、ゲージ長だけ離れた２点の位相の差分を計算し、位相接続処理などを行い、振動波形を計算する。具体的な計算方法は従来方法と同様である。
　本開示の特徴は信号処理部１７ｄと１７ｅにある。信号処理部１７ｄでは、測定で得られた信号Ｉｉ及びＱｉを時間ポイント数Ｍのブロックに区切り、測定時間Ｍｔのブロック単位で回転角度αｉを計算・更新することで、ＣＷ光源１に備わるレーザの発振周波数などの光学特性の時間的な変化や被測定光ファイバ６自体の温度変化、大きな動的歪みの被測定光ファイバ６への印加などによる最適値の時間的な変化に対応する。一ブロックあたりの測定時間Ｍｔを、ＣＷ光源１に備わるレーザの発振周波数などの光学特性の時間的な変化や被測定光ファイバ６自体の温度変化、大きな動的歪みの被測定光ファイバ６への印加の時間スケールに対して小さくとることで、十分な細かさでの回転角度の更新が可能なようにする。
　このような設定にしておけば、信号処理部１７ｅでブロックＫに属する散乱光ベクトルｒｉの周波数平均を計算する際に用いる回転角度αｉとして、信号処理部１７ｄで計算した直前のブロックＫ－１の回転角度αｉを用いても問題なくなる。また、信号処理部１７ｄが直前のブロックＫ－１の回転角度αｉを信号処理部１７ｅに出力することで、信号処理部１７ｄが散乱光ベクトルｒｉ＿ａｖｅｔ（ｚ）を更新した段階で、散乱光ベクトルｒｉ（ｋｔ，ｚ）のデータを信号処理装置１７のメモリに保持する必要がなくなる。
　このように、本開示は、ＡＤ変換機能素子１７ａ及び１７ｂ、光周波数分離部として機能する信号処理部１７ｃ、及び回転角度計算部として機能する信号処理部１７ｄにおける処理が、独立かつ並列にストリーミング処理を実行可能である。これにより、本開示は、信号処理部１７ｅでブロックＫに属する散乱光ベクトルｒｉの周波数平均を計算する際に用いる回転角度αｉとして、ブロックＫで信号処理部１７ｄで計算した回転角度自体を用いる従来技術では信号処理部１７ｄでブロックＫに属する散乱光ベクトルのデータを全て計算機のメモリに保持する必要があるのに比べて、計算機のメモリを削減することができる。
　ただし、一ブロックあたりの測定時間Ｍｔの長さは回転角度αｉの計算精度にも関係する。一ブロックあたりの測定時間Ｍｔを、十分な精度で回転角度αｉが計算できるだけ長くとり、ＣＷ光源１に備わるレーザの発振周波数などの光学特性の時間的な変化や被測定光ファイバ６自体の温度変化、大きな動的歪みの被測定光ファイバ６への印加の時間スケールに対しては小さくとることで、回転角度αｉの計算精度も確保できる。このような設定は、多くの状況において可能である。例えば、Ｍを１００程度にとれば十分な回転角度αｉの精度が得られることが実験的な知見として得られているが（非特許文献３）、光パルスの送出周期ｔを１ｍｓとしてもＭｔは１００ｍｓ程度となり、これはＣＷ光源１に備わるレーザの発振周波数などの光学特性の時間的な変化や被測定光ファイバ６自体の温度変化、大きな動的歪みの被測定光ファイバ６への印加の時間スケールが１ｓ以上の場合には十分に細かいと考えられる。
　本開示の装置はコンピュータとプログラムによっても実現でき、プログラムを記録媒体に記録することも、ネットワークを通して提供することも可能である。本開示のプログラムは、本開示に係る装置に備わる各機能部としてコンピュータを実現させるためのプログラムであり、本開示に係る装置が実行する方法に備わる各ステップをコンピュータに実行させるためのプログラムである。
１：ＣＷ光源
２：カプラ
３：光変調器
４：光パルス
５：サーキュレータ
６：被測定光ファイバ
７：９０度光ハイブリッド
１３、１４：バランス検出器
１５、１６：電気信号
１７：信号処置装置
１７ａ、１７ｂ：ＡＤ変換機能素子
１７ｃ、１７ｄ：信号処理部
３１：測定装置

Claims (7)

  1. 　光周波数の異なる複数の光パルスが光ファイバで散乱された散乱光の測定データを取得し、
    　前記光ファイバの任意の地点を測定可能な予め定められた時間範囲における、各光周波数での散乱光の散乱光ベクトルｒｉを計算し、
    　各光周波数の散乱光ベクトルｒｉを前記時間範囲で平均し、
    　前記平均によって得られた時間平均ベクトルｒｉ＿ａｖｅｔを用いて、各光周波数の位相を補正するための回転角度αｉを計算し、
    　各光周波数の散乱光ベクトルｒｉの位相を回転角度αｉ回転させることで、各光周波数の散乱光ベクトルｒｉの位相を補正する、
    　信号処理装置。
  2. 　前記時間範囲が、測定データの取得時間を複数のブロックに区切ったときの、任意の１つのブロックの時間範囲であり、
    　ブロックごとに各光周波数の散乱光ベクトルｒｉの位相を補正する、
    　請求項１に記載の信号処理装置。
  3. 　１つ前のブロックで計算された前記回転角度αｉを用いて、前記散乱光ベクトルｒｉの位相を補正する、
    　請求項２に記載の信号処理装置。
  4. 　散乱光の同相成分と直交成分を個別にデジタル信号に変換するＡＤ変換機能素子と、
    　前記ＡＤ変換機能素子から出力されたデジタル信号の同相成分Ｉｄｉｇｉｔａｌと直交成分Ｑｄｉｇｉｔａｌを、各光周波数の信号Ｉｉと信号Ｑｉに分離する光周波数分離部と、
    　前記光周波数分離部から出力された信号Ｉｉと信号Ｑｉを用いて散乱光ベクトルｒｉを計算し、計算した散乱光ベクトルｒｉを前記時間範囲で平均し、前記平均によって得られた時間平均ベクトルｒｉ＿ａｖｅｔを用いて前記回転角度αｉを計算する回転角度計算部と、
    　前記回転角度計算部で計算された前記回転角度αｉを用いて、前記回転角度計算部で計算された前記散乱光ベクトルｒｉの位相を補正する位相補正部と、
    　を備え、
    　前記ＡＤ変換機能素子、前記光周波数分離部及び前記回転角度計算部における処理が、独立かつ並列に実行可能である、
    　請求項１に記載の信号処理装置。
  5. 　回転角度αｉを回転後のベクトルＲｉを周波数平均し、
    　周波数平均後の周波数平均ベクトルＲａｖｅｆの角度θａｖｅｆを計算し、
    　角度θａｖｅｆを用いて振動波形を計算する、
    　請求項１に記載の信号処理装置。
  6. 　光周波数の異なる複数の光パルスが光ファイバで散乱された散乱光を測定する測定装置と、
    　前記測定装置で測定された散乱光の測定データを取得し、前記光ファイバの任意の地点での散乱光の散乱光ベクトルｒｉを計算し、計算した散乱光ベクトルｒｉの位相変化を計算する、請求項１から５のいずれかに記載の信号処理装置と、
    　を備える測定システム。
  7. 　光周波数の異なる複数の光パルスが光ファイバで散乱された散乱光の測定データの信号処理を行う信号処理方法であって、
    　前記光ファイバの任意の地点を測定可能な予め定められた時間範囲における、各光周波数での散乱光の散乱光ベクトルｒｉを計算し、
    　各光周波数の散乱光ベクトルｒｉを前記時間範囲で平均し、
    　前記平均によって得られた時間平均ベクトルｒｉ＿ａｖｅｔを用いて、各光周波数の位相を補正するための回転角度αｉを計算し、
    　各光周波数の散乱光ベクトルｒｉの位相を回転角度αｉ回転させることで、各光周波数の散乱光ベクトルｒｉの位相を補正する、
    　信号処理方法。
PCT/JP2022/034477 2022-09-14 2022-09-14 位相ｏｔｄｒにおける信号処理方法 WO2024057462A1 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
PCT/JP2022/034477 WO2024057462A1 (ja) 2022-09-14 2022-09-14 位相ｏｔｄｒにおける信号処理方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2022/034477 WO2024057462A1 (ja) 2022-09-14 2022-09-14 位相ｏｔｄｒにおける信号処理方法

Publications (1)

Publication Number Publication Date
WO2024057462A1 true WO2024057462A1 (ja) 2024-03-21

Family

ID=90274612

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2022/034477 WO2024057462A1 (ja) 2022-09-14 2022-09-14 位相ｏｔｄｒにおける信号処理方法

Country Status (1)

Country Link
WO (1) WO2024057462A1 (ja)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2020169904A (ja) * 2019-04-03 2020-10-15 日本電信電話株式会社 位相測定方法及び信号処理装置
US20200370949A1 (en) * 2019-05-22 2020-11-26 Nec Laboratories America, Inc Performance of rayleigh-based phase-otdr with correlation-based diversity combining and bias removal
US20220120608A1 (en) * 2020-10-21 2022-04-21 Nec Laboratories America, Inc Polarization diversity combining method in coherent das maintaining phase continuity

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2020169904A (ja) * 2019-04-03 2020-10-15 日本電信電話株式会社 位相測定方法及び信号処理装置
US20200370949A1 (en) * 2019-05-22 2020-11-26 Nec Laboratories America, Inc Performance of rayleigh-based phase-otdr with correlation-based diversity combining and bias removal
US20220120608A1 (en) * 2020-10-21 2022-04-21 Nec Laboratories America, Inc Polarization diversity combining method in coherent das maintaining phase continuity

Similar Documents

Publication Publication Date Title
WO2020203332A1 (ja) 位相測定方法及び信号処理装置
JP4097758B2 (ja) 光信号の伝送方法、光伝送システム、および光伝送システムの監視装置
JP6277147B2 (ja) 光ファイバ振動測定方法及びシステム
US7426021B2 (en) Interferometric optical analyzer and method for measuring the linear response of an optical component
JP7265057B2 (ja) Ｏｆｄｒ呼掛け監視及び最適化のための方法及び装置
EP2972409A1 (en) System and method for increasing coherence length in lidar systems
JP7435160B2 (ja) 光ファイバ振動検知装置及び振動検知方法
WO2020008886A1 (ja) 分布光ファイバ振動計測装置および分布光ファイバ振動計測方法
JP2020518146A (ja) スペクトル相関を使用した変調信号のスペクトル形状の無雑音測定
US6407805B1 (en) Time delay measurement using optical signals modulated at multiple frequencies
WO2024057462A1 (ja) 位相ｏｔｄｒにおける信号処理方法
CN113316710B (zh) 振动检测方法、信号处理装置和存储介质
CN114448500A (zh) 相频响应测量方法和装置
JP7156132B2 (ja) 光コヒーレントセンサ及び光コヒーレントセンシング方法
US20220006521A1 (en) Phase measurement method, signal processing device, and program
WO2022259437A1 (ja) 振動測定器及び振動測定方法
EP2972555B1 (en) System and method for increasing coherence length in lidar systems
NO330324B1 (no) Fremgangsmate for a beregne et mal for lysforplantningstidsdifferansen for to lysforplantende veier gjennom et lysforplantende medium
WO2020234989A1 (ja) 位相測定方法及び信号処理装置
JP2011043344A (ja) 光信号波形計測装置
WO2022259436A1 (ja) 信号処理装置、振動検出システム及び信号処理方法
WO2022059050A1 (ja) 光周波数領域反射計測装置及び方法
RU2664692C1 (ru) Измеритель фазовых шумов узкополосных лазеров, основанный на состоящем из рм-волокна интерферометре маха-цендера
WO2023053263A1 (ja) 光パルス試験方法及び光パルス試験装置
WO2024024094A1 (ja) モード間群遅延差測定装置及び方法

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 22958788

Country of ref document: EP

Kind code of ref document: A1