WO2023058160A1 - レイリー強度パターン計測装置およびレイリー強度パターン計測方法 - Google Patents

レイリー強度パターン計測装置およびレイリー強度パターン計測方法 Download PDF

Info

Publication number
WO2023058160A1
WO2023058160A1 PCT/JP2021/037005 JP2021037005W WO2023058160A1 WO 2023058160 A1 WO2023058160 A1 WO 2023058160A1 JP 2021037005 W JP2021037005 W JP 2021037005W WO 2023058160 A1 WO2023058160 A1 WO 2023058160A1
Authority
WO
WIPO (PCT)
Prior art keywords
rayleigh
intensity pattern
data
correlation coefficient
frequency
Prior art date
Application number
PCT/JP2021/037005
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/JP2021/037005 priority Critical patent/WO2023058160A1/ja
Priority to JP2023552604A priority patent/JPWO2023058160A1/ja
Publication of WO2023058160A1 publication Critical patent/WO2023058160A1/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

  • This application relates to a Rayleigh intensity pattern measuring device and a Rayleigh intensity pattern measuring method.
  • DFOS Distributed Fiber Optic Sensing
  • DFOS Distributed Fiber Optic Sensing
  • Typical methods of this distributed optical fiber sensing technology include Raman, Brillouin, and Rayleigh methods.
  • DAS in civil engineering or other companies' prior art can measure dynamic changes with high accuracy, but they are based on electrical signals, not on optical fiber measurements. Therefore, reproducibility is poor and it is not suitable for long-term management.
  • the method using optical fiber has reproducibility and can be said to be a more suitable method for long-term management.
  • representative methods of the DFOS technology include the above-described three methods, and among them, the Brillouin method and the Rayleigh method are particularly promising.
  • the Brillouin method can measure the absolute frequency, the accuracy is often insufficient for the required accuracy, and the Rayleigh method is currently the most effective. is assumed.
  • the Rayleigh method has attracted attention due to its many features such as high stability, signal strength, availability in both multimode and single mode, high precision, and high resolution of 1 mm or less. That is, compared to the other two methods, the Rayleigh method is (1) superior in measurement stability and signal strength, and can be used in either single mode or multimode; (3) high resolution of 1 mm or less;
  • FIG. 1 shows the representative performance of each of the six prior arts described in the vertical direction, divided into six indexes described in the horizontal direction, as a list. (See Patent Document 1 and Non-Patent Documents 1 to 5).
  • the presence or absence of the absolute frequency meaning the frequency of the LD whose frequency can be controlled and the frequency error is clear; the same applies hereinafter
  • the measurement range is sufficient
  • the effectiveness of non-uniform strain distribution even if there is non-uniform strain distribution, there is no problem in measurement, in other words, poor correlation in a relatively small strain range
  • FIG. It is a case-by-case method of dependence, and there is no standard to continue to the next model.
  • the effectiveness of the non-uniform strain distribution can be obtained up to several tens of ⁇ m as a resolution. Gender is "Yes”.
  • the TGD-DAS technology shown in the sixth row is also effective because it can be directly converted into distortion if the phase change is targeted.
  • the measurement distance is 10 km or more
  • the performance of velocity and spatial resolution is above the standard
  • the absolute frequency, sufficient measurement range, and effectiveness of non-uniform strain distribution are all satisfied in this table. It can be seen that the six prior arts do not.
  • FIG. 2 is a diagram showing a comparison of the above six prior arts related to DFOS by classifying the characteristic techniques used in these six prior arts into five types (Patent Document 1 and See Non-Patent Documents 1-5).
  • a frequency control means for example, Etalon, Gas Cell, etc.
  • RBS Rayleigh Backscattering Spectrum
  • optical fibers as material coordinates (coordinates embedded in a moving object), to be recognizable as a position indicator (to indicate a place), and to be detectable without being an optical signal.
  • material coordinates coordinates embedded in a moving object
  • position indicator to indicate a place
  • optical signal There are three features: that there is practically no limit to distance resolution and measurement accuracy.
  • the Rayleigh backscattered light spectrum (RBS) method applied to the measurement of the managed object uses the position and frequency as variables, the strength of the electric field (hereinafter also referred to as the electric field), its absolute value, or the absolute value Since it is essentially the same as a Rayleigh Intensity Pattern (hereinafter also referred to as RIP) whose function value is the square, it will be referred to as a Rayleigh intensity pattern measuring device in the present application and will be described below. .
  • the electric field the strength of the electric field
  • RIP Rayleigh Intensity Pattern
  • the spatial position is matched, the cross-correlation between the RIP value at time t1 and the RIP value at time t2 is obtained, and from this correlation peak position, the Rayleigh frequency shift amount Ask for That is, the Rayleigh frequency shift is measured by power correlation.
  • FIG. 3 is a diagram showing an example of RIP distribution obtained by a Rayleigh backscattered light spectrum method obtained by a full-band tunable wavelength distributed feedback LD (LD: Laser Diode).
  • LD Laser Diode
  • the Rayleigh frequency shift amount ⁇ R represented by the formula (F S ) is obtained (see, for example, Patent Document 2).
  • ⁇ R C 21 ⁇ +C 22 ⁇ T+C 23 ⁇ P (F S )
  • C 21 , C 22 , and C 23 are sensitivity coefficients given to Rayleigh frequency shift
  • ⁇ , ⁇ T, and ⁇ P are the strain, temperature change, and pressure of the optical fiber at the measurement position, respectively. Show change.
  • the frequency shifts derived by strain, temperature and pressure changes mentioned above are obtained by advanced correlation analysis (hence the frequency shifts are relative values).
  • the signal used for correlation analysis is the power spectrum over a specific range of distances. and is phase insensitive. Therefore, the above frequency shift can be considered as a distinct feature comparison by frequency of the TW-COTDR system in two states (reference time and observation time) at the same location.
  • the TW-COTDR method and the PPP-BOTDA (Pulse Pre-Pump Brillouin Optical Time Domain Analysis) method which is a high-precision Brillouin technique, are decoded to obtain two methods, the Rayleigh backscattered light method and the Brillouin backscattered light method.
  • the Rayleigh backscattered light method in the time domain has the following advantages over the Brillouin backscattered light method. small changes such as creep monitoring over long distances can be accurately and easily detected; and the signal from the optical fiber is phase independent, resulting in stability and stability. It must have high reproducibility.
  • the prior art including the TW-COTDR system has the following problems.
  • the "absolute frequency (absolute reference frequency)" which is the first index shown in FIG. 1, has not been established. None of the prior arts shown in Prior Art No. 2 to Prior Art No. 6 can clarify the error by controlling the LD frequency or wavelength value as shown in FIG. Also, in the TW-COTDR system shown in Prior Art No. 1, although the wavelength (optical frequency) of the LD is controlled, the absolute value is not managed. Therefore, only relative frequency shift values can be obtained in measurements between different devices (of the same manufacturer) or between devices of different manufacturers.
  • the second indicator is not sufficient (a sufficient measurement range cannot be secured for practical use).
  • the measurement range is too small. Specifically, for example, a temperature change of 100 ° C. is a common case, but in each method using chirped pulse technology (each prior art other than the TW-COTDR method of prior art number 1), the required numerical value It can only be measured in the range of 4 digits or less.
  • FIG. 4A shows measurement results when the frequency shift (hereinafter also referred to as RFS) of Rayleigh backscattered light is measured using a bar with the following measurement parameters. Measurement parameters; sweep start frequency: 192010 GHz, sweep period: 400 GHz, sweep step: 0.2 GHz, spatial resolution: 2 cm, sampling interval: 1 cm.
  • RFS frequency shift
  • the strain is not correctly detected in the area enclosed by the dotted line and surrounded by the oval (within the range of 512.3 m to 513.1 m in the distance position shown on the horizontal axis).
  • the correlation coefficient between the two data that is, the reference data and the measured data is 0.2 or less, as shown in FIG. indicates that there is no
  • data are accumulated in the time direction to finally obtain the deformation. This is because the measurement loses the correlation with the initial value (initial strain change amount).
  • FIG. 5 shows data showing the measurement results when the strain generated in the bar based on the Brillouin scattered light was measured using the same bar as described above with the following measurement parameters. Measurement parameters; sweep start frequency: 10 GHz, sweep period: 1.45 GHz, sweep step: 5 MHz, spatial resolution: 5 cm, sampling interval: 2.5 cm.
  • Fig. 5 shows the distribution of strain when distance (unit: m) is taken on the horizontal axis with respect to strain data of bars measured using the Brillouin backscattered light method. From FIG. 5, it can be seen that the strain value, which was 100 ⁇ or more near the distance of 512.2 m, changed to ⁇ 400 ⁇ near the distance of 512.5 m, and the strain change of 500 ⁇ or more was measured.
  • the unit step of frequency (change) during data processing is 0.04 GHz.
  • TW-COTDR (Prior Art No. 1) uses a full-band LD, so the measurement range is sufficiently wide, but there is a problem that the line width of the LD is too wide.
  • the effective coherence length is only several tens of meters, and in TGD (Time Gated Digital) phase analysis, the phase signal contains large phase noise, so it cannot be used as it is. That is, in the TW-COTDR system, as shown in Fig. 2, the homodyne reception is adopted as the reception system, so the phase information may be lost from the beginning, and the phase information necessary for advanced signal analysis is used. Can not.
  • the reception frequency range (reception bandwidth) of the RIP must be 100 GHz or more, but it is almost impossible to cover the reception band for electrical signals. This is the current situation.
  • the problems described above can be summarized as follows with specific numerical examples for each representative prior art.
  • the TW-COTDR method of prior art number 1 is excellent in static measurement of strain, temperature, etc., and is characterized by high spatial resolution (2 cm) and high accuracy measurement method.
  • the disadvantage is that the time is as slow as several tens of seconds.
  • the TGD-DAS method of Prior Art No. 6 is excellent in dynamic measurement such as sound wave measurement, and is characterized by a high measurement speed of 2000 times or more per second, but the spatial resolution is about 1 m, which is inferior.
  • it is limited to minute change measurement, for example, the temperature change is up to 0.2° C., and it is not suitable for large change measurement.
  • the present application was made to solve the above problems, and three parameters that have a trade-off relationship in realizing high-speed and high-precision measurement: spatial resolution, measurement speed, and sufficient It is an object of the present invention to provide a Rayleigh intensity pattern measuring device capable of realizing a wide range of measurement at the same time and managing data acquired over a long period of time.
  • the Rayleigh intensity pattern measurement device disclosed in the present application is a light source unit having a variable wavelength LD and a controller capable of changing the frequency of laser light emitted from the variable wavelength LD; an optical coupler for inputting a laser beam oscillated from the light source into an optical fiber and outputting Rayleigh scattered light from the optical fiber along a path different from the incident path of the laser beam; a receiver that coherently receives the local light from the light source and the Rayleigh scattered light from the optical coupler; an AD converter for AD-converting the signal output from the receiving unit; a first calculation unit for calculating the signal converted by the AD converter by a predetermined method; A memory database to which a signal is input and stored, and a second computing unit that performs a predetermined computation based on the data stored in the memory database, and is obtained from the electric field signal of the Rayleigh scattered light , from two different Rayleigh intensity patterns at a predetermined measurement position, the cross-correlation coefficient is obtained after performing a predetermined correction on the measurement position in the second calculation unit, and
  • the Rayleigh intensity pattern measurement device disclosed in the present application, three parameters that have a trade-off relationship in realizing high-speed and high-precision measurement, spatial resolution, measurement speed, and a sufficient measurement range are simultaneously realized. It is an object of the present invention to provide a Rayleigh intensity pattern measuring device capable of managing data acquired over a long period of time.
  • FIG. 2 is a diagram for explaining technical indicators to be a problem of the Rayleigh intensity pattern measuring device of Embodiment 1;
  • FIG. 4 is a diagram for explaining a technology characteristic of DFOS;
  • FIG. 4 is a diagram showing an example of RIP distribution obtained by the Rayleigh backscattering method;
  • FIG. 5 shows an example of bar strain and correlation function measured by the Rayleigh backscatter method.
  • FIG. 4 shows bar strain data measured by the Brillouin backscatter method.
  • 1 is a diagram for explaining the overall configuration of a RIP measuring device according to Embodiment 1;
  • FIG. FIG. 4 is a diagram for explaining the operation of the RIP measuring device according to the first embodiment;
  • FIG. 5 is a diagram for explaining peak values of cross-correlation coefficients of Rayleigh frequency shifts in the non-uniform distortion detection and optimization analysis method of the RIP measuring apparatus according to the first embodiment;
  • FIG. 4 is a diagram for explaining Fresnel integrals used in a non-uniform strain detection optimization analysis method;
  • FIG. 4 is a diagram for explaining a finite chirp waveform and its Fourier transform when the distortion change of non-uniform distortion is 0.1 GHz/m;
  • FIG. 4 is a diagram for explaining a finite chirp waveform and its Fourier transform when the distortion change of non-uniform distortion is 5 GHz/m;
  • FIG. 4 is a diagram for explaining a finite chirp waveform and its Fourier transform when the distortion change of non-uniform distortion is 20 GHz/m;
  • FIG. 10 is a diagram showing a strain model and a strain gradient model in the non-uniform strain detection optimization analysis method;
  • FIG. 10 is a diagram showing a Rayleigh frequency shift model and its gradient model in an analysis method for proper detection of non-uniform distortion;
  • FIG. 10 is a diagram for explaining simulation results of a model that does not consider the inclination of Rayleigh frequency shift in detecting non-uniform distortion;
  • FIG. 10 is a diagram for explaining a simulation result when a Rayleigh frequency shift distribution is obtained by applying an analysis method for proper detection of non-uniform distortion;
  • FIG. 10 is a diagram showing a finite chirp waveform and its Fourier transform when the distortion change of non-uniform distortion is 20 GHz/m;
  • FIG. 10 is a diagram showing a strain model and a strain gradient model in
  • FIG. 10 is a diagram for explaining a simulation result when a distribution of slopes of Rayleigh frequency shifts is obtained by applying an analysis method for proper detection of non-uniform distortion
  • FIG. 10 is a diagram for explaining a simulation result when a strain distribution is obtained by applying an analysis method for proper detection of non-uniform strain
  • FIG. 5 is a diagram for explaining strain distribution data for forming a Rayleigh intensity pattern used for explaining an example in the RIP measuring apparatus according to the first embodiment
  • FIG. 5 is a diagram for explaining parameters used when explaining a scanning method of a plurality of chirped pulses simulating laser light, which are used to drive a receiver in a Rayleigh intensity pattern measuring device;
  • FIG. 4 is a model diagram for explaining a method of stepwise scanning with a plurality of chirped pulses;
  • FIG. 4 is a diagram for explaining the form of a matched filter for extracting subbands from each chirped pulse;
  • FIG. 4 is a diagram for explaining a series of methods from each chirped pulse to extracting a trace showing a specific frequency along the distance of the object to be measured;
  • FIG. 10 shows an example of RIP, which is the intensity of the Rayleigh backscatter signal at the corresponding laser frequency divided by the laser linewidth.
  • FIG. 10 is a diagram for explaining a data preparation step in another analysis method for proper detection of non-uniform strain;
  • FIG. 10 is a diagram for explaining a data preparation step in another analysis method for proper detection of non-uniform strain;
  • FIG. 10 is a diagram for explaining a data decomposition step in another analysis method for proper detection of non-uniform distortion;
  • FIG. 10 is a diagram for explaining a correlation coefficient calculation and determination process in another analysis method for proper non-uniform strain detection;
  • FIG. 10 is a diagram showing a data processing flowchart in another analysis method for proper detection of non-uniform distortion.
  • FIG. 10 is a diagram showing results of processing data in which non-uniform strain is detected using another analysis method for non-uniform strain detection optimization;
  • FIG. 30 is a supplementary explanatory diagram of FIG. 29;
  • FIG. 10 is a diagram for explaining a fundamental frequency component used in a method similar to another analysis method for proper detection of non-uniform distortion;
  • FIG. 10 illustrates a data processing flow used in an alternative analysis method-like approach for non-uniform strain detection and optimization
  • FIG. 10 is a diagram showing a correlation coefficient, which is an example of the result of applying a technique similar to another analysis method for proper detection of non-uniform distortion.
  • FIG. 5 is a diagram showing a strain distribution that is an example of the result of applying a method similar to another analysis method for proper detection of non-uniform strain.
  • FIG. 10 is a diagram for explaining the effect of a method similar to another analysis method for proper detection of non-uniform distortion;
  • Embodiment 1 The overall configuration of the Rayleigh intensity pattern measuring device (hereinafter also referred to as RIP measuring device for simplification) according to the first embodiment will be described below with reference to FIGS. 6 and 7.
  • FIG. 6 The overall configuration of the Rayleigh intensity pattern measuring device (hereinafter also referred to as RIP measuring device for simplification) according to the first embodiment will be described below with reference to FIGS. 6 and 7.
  • FIG. 6 The overall configuration of the Rayleigh intensity pattern measuring device (hereinafter also referred to as RIP measuring device for simplification) according to the first embodiment will be described below with reference to FIGS. 6 and 7.
  • FIG. 1 The overall configuration of the Rayleigh intensity pattern measuring device (hereinafter also referred to as RIP measuring device for simplification) according to the first embodiment will be described below with reference to FIGS. 6 and 7.
  • FIG. 6 is a diagram for explaining the overall configuration of the RIP measuring device 100 according to the first embodiment.
  • the RIP measurement apparatus 100 according to the first embodiment is a TW (variable wavelength)-TGD-RFS (Rayleigh frequency shift) type measurement apparatus that uses a chirped pulse by dividing it.
  • the RIP measuring apparatus 100 according to the first embodiment roughly includes four main components. They are a light source unit 11, a transmission/reception optical unit 12, a reception unit 13, and a RIP determination digital processing unit 14 (hereinafter also referred to as a Rayleigh intensity pattern digital processing unit 14 or RIP digital processing unit 14). These details will be described in order below.
  • the light source unit 11 has an LD1 that is a broadband wavelength tunable laser (TW-LD) and a controller that controls laser light emitted from the LD1, and performs automatic temperature control and automatic frequency control, while controlling the built-in optical Absolute frequency control of LD1 by a local oscillator (OLO) (here, absolute frequency control is based on the absorption wavelength of molecules, for example, like Gas Cell.
  • OLO local oscillator
  • absolute frequency control is based on the absorption wavelength of molecules, for example, like Gas Cell.
  • the information is publicly known, and the temperature dependence is also It is low and can be easily realized within a range that does not affect the measurement accuracy), and frequency step control is performed to oscillate divided chirp pulses as an output.
  • the wavelength tunable laser also includes an external modulation type.
  • the light source unit 11 uses a chirped light pulse to control the frequency of the wavelength-variable LD 1, and sweeps while changing the frequency in a plurality of steps (for example, changing in 12 steps at 200 GHz). do.
  • the TW-LD is used to generate a high-precision optical pulse using a TGD (Time Gated Digital) method or the like.
  • the sweep width in the TGD method is, for example, 4 GHz (sample rate is 8 GS/s).
  • the transmission/reception optical unit 12 includes an optical switch 2 that converts a laser signal into pulses, an erbium-doped fiber amplifier 3 (abbreviated as EDFA 3 hereinafter) that amplifies the output signal of the optical switch 2, and an EDFA 3.
  • the signal from is transmitted to the object to be measured (hereinafter also referred to as the object to be measured (specifically, optical fiber)) through the reference optical fiber 5, and the scattered light from the object to be measured of the transmitted signal After receiving the backscattered light through the reference optical fiber 5, it has an optical coupler 4 for transmitting it to an external receiver, which will be described later.
  • the receiving section 13 has a receiver that receives the backscattered light emitted from the optical coupler 4 of the transmitting/receiving optical section 12 and receives the light source light from the light source section 11 .
  • this receiver consists of a heterodyne receiver that receives signals using the polarization diversity system, or a homodyne receiver that receives signals using the phase diversity system and the polarization diversity system. It is
  • the RIP digital processing unit 14 is a part that digitally analyzes the intensity or phase of a signal.
  • the output signal from the receiving unit 13 is converted to digital at high speed by a 2ch AD converter 6 and received, and FPGA (Field Programmable Gate Array), ASIC (Application Specific Integrated Circuit), etc. are used to convert 4 to be described later.
  • FPGA Field Programmable Gate Array
  • ASIC Application Specific Integrated Circuit
  • the RIP digital processing unit 14 has a processor for arithmetic processing of digital signals, and a memory database 8 (hereinafter also referred to as a memory DB 8) for storing and saving the signals processed by this processor.
  • a first computing unit 7 that computes a signal received from (for example, performs a computation of synthesizing a plurality of chirp signals) and outputs it to a memory database 8 as a signal different from the received signal, and this memory DB 8 Select data at two times (for example , time t i and time t j ) for correlation analysis from the information stored in . , resampling such as interpolation processing on the frequency axis, distance correction, correlation (correlation judgment), and other arithmetic processing, and output as a determined Rayleigh frequency shift ⁇ RD (this second arithmetic unit 9).
  • the 2ch AD converter 6 can convert, for example, the amplitude component and the phase component of the signal to each channel using a processor such as an FPGA chip.
  • the RIP digital processing unit 14 can perform signal strength and phase analysis using a digital method.
  • the output of the determined Rayleigh frequency shift ⁇ RD is included inside the RIP digital processing unit 14 , but it may be outside the RIP digital processing unit 14 .
  • the light source unit 11, the receiving unit 13, and the RIP digital processing unit 14, which are surrounded by a solid line frame, are particularly the RIP measurement apparatus according to the first embodiment. It is a component that characterizes the on the other hand, the transmitting/receiving optical unit 12 surrounded by a dotted line frame can be applied with existing OTDR technology, and is not particularly new (see, for example, Patent Document 3).
  • the chirp signal which is the signal heterodyned received by the receiving unit 13, is not subjected to any special processing in the 2ch AD converter 6 (the signal is passed through), and is processed as it is in the memory of the RIP digital processing unit 14.
  • Store in DB 8 (step S1). Since data are stored in the memory/DB 8 in chronological order tk , data at two times, for example, time t i and time t j are selected from this data (step S4). Based on the data at the two selected times, the second calculation unit 9 executes the following series of operations.
  • an initial Rayleigh frequency shift ⁇ RI is calculated (step S5).
  • the obtained initial Rayleigh frequency shift ⁇ RI is subjected to frequency resampling processing including interpolation processing for distance correction (step S6).
  • a cross-correlation coefficient after distance correction is obtained based on the frequency-resampled data (step S7).
  • it is determined whether or not the obtained cross-correlation coefficient value is greater than a predetermined threshold value (step S8). This threshold can be empirically predetermined as, for example, 0.3. If the result of the determination is "Yes", the Rayleigh frequency shift ⁇ R is obtained based on the value of the obtained cross-correlation coefficient, and the value of the obtained Rayleigh frequency shift ⁇ R is determined . D (step S10).
  • step S9 the data used as the basis for recalculating the cross-correlation coefficient after distance correction is obtained by analysis (step S9), and the cross-correlation coefficient after distance correction is obtained again (step S7).
  • step S8 it is determined whether or not the obtained cross-correlation coefficient value is greater than a predetermined threshold value (step S8). Then, the operations of step S9 ⁇ step S7 ⁇ step S8 are repeated until the value of the cross-correlation coefficient becomes larger than the predetermined threshold value, and the value of the cross-correlation coefficient becomes larger than the predetermined threshold value.
  • the Rayleigh frequency shift ⁇ R is determined based on the value of the cross-correlation coefficient in the case, and the obtained value of the Rayleigh frequency shift ⁇ R is defined as the fixed Rayleigh frequency shift ⁇ RD (step S10).
  • Method 1 since there is no processing by the 2ch AD converter 6, high-speed processing of processing data is possible.
  • the chirp signal which is the signal heterodyned received by the receiving unit 13, is AD-converted by the 2ch AD converter 6, and the signal after comb filtering processing by FPGA or the like is stored in the memory/DB 8 of the RIP digital processing unit 14.
  • Store step S2 Since subsequent processing is the same as method 1, detailed description is omitted here.
  • Method 2 by performing comb filtering processing, periodic noise contained in the signal can be removed, and the accuracy of the measurement data can be improved.
  • method 3 After processing in method 2, the signal after digital processing obtained by squaring the processed signal is stored in the memory DB 8 of the RIP digital processing unit 14 (step S3). Since subsequent processing is the same as method 1, detailed description is omitted here. Since the squared signal is used in method 3, the influence of the phase among the signal components is reduced, and the precision of the measurement data can be improved accordingly.
  • Method 4 will be described, focusing on the differences from Method 1, Method 2, and Method 3.
  • a chirp signal that is a signal heterodyned-received by the receiver 13 is used.
  • signal is used.
  • a phase diversity scheme is employed to prevent loss of phase information.
  • a signal obtained by squaring the homodyne-received signal is stored in the memory DB 8 of the RIP digital processing unit 14 (step S3). Since the operation after the above is the same as method 1, detailed description is omitted here.
  • Method 4 requires a narrow linewidth LD and an optical PLL with less phase noise than methods 1 to 3, which use heterodyne received signals. It has the advantage of being easier to process data at high speed.
  • the pulse power can be effectively improved, and the scattered waveforms at a plurality of frequencies can be acquired simultaneously.
  • the number of pulses can be greatly reduced, making it possible to shorten the measurement time.
  • the RIP digital processing unit 14 performs a large amount of analysis.
  • the methods 1 to 4 described above are devised so that the non-uniform strain can be properly detected using an analysis method for detection optimization as shown in FIG. Therefore, an analysis method capable of properly detecting the non-uniform strain, which is commonly used in methods 1 to 4, will be described in detail below.
  • This analysis method is roughly divided into two types, an analysis method N1 and an analysis method N2.
  • the analysis method N1 will be sequentially described below.
  • Analysis method N1 is a reference spectrum obtained from the electric field of the reference data of Rayleigh scattered light, various reference spectra obtained by modifying it, and the observation spectrum obtained from the electric field of the observation data, even if there is non-uniform distortion. It is an analysis method that can properly detect the strain distribution of the object to be measured from the cross-correlation of .
  • the deformation of the reference spectrum is to deal with the strain slope (described in detail later), and the deformation varies according to the magnitude of the slope.
  • the processing method (algorithm) used in this analysis method N1 will be described in detail below.
  • parameters are determined based on the underlying model.
  • a laser beam having a frequency ⁇ and a pulse width D is incident from the input end of an optical fiber, and the scattered and returned light is observed, particularly the Rayleigh scattered light P is observed.
  • the observed values are the intensity of the light obtained by direct detection or the intensity of the electric field obtained by heterodyne reception.
  • the frequency ⁇ For example, spectral data of Rayleigh scattered light can be obtained by obtaining observation values while scanning the frequency ⁇ of the laser light source.
  • the distance interval [z ⁇ l p /2, z+l p /2] (hereafter, the distance interval is simply referred to as the interval)
  • the Rayleigh frequency shift ⁇ R (z) is a constant value in this interval.
  • Y 0 (z, u) and Y 1 (z, u) are given by equations (6) and (7), respectively. Note that heterodyne detection must be performed in order to use the electric field data.
  • the strain or temperature changes linearly in the interval corresponding to the pulse width
  • the Rayleigh frequency shift is expressed in the interval [z ⁇ l p /2, z+l p /2] by the following equation (8)
  • a is the Rayleigh frequency shift at the center of the interval
  • b is the slope of the Rayleigh frequency shift within the interval.
  • Equation (10) is the Fourier transform of the product of two functions ⁇ (z+ ⁇ ), ⁇ l p /2 ⁇ l p /2 and exp( ⁇ i ⁇ 2 ), ⁇ , with respect to ⁇ . ing.
  • ⁇ (z) is the Rayleigh backscatter coefficient at distance z
  • the Fourier transform of the product of two functions is equal to the convolution of the Fourier transforms of their respective functions.
  • the Fourier transform of ⁇ (z+ ⁇ ), ⁇ l p /2 ⁇ l p /2 is equal to Y 0 (z, u) from equation (6). Therefore, the following equation (11) holds except for constant coefficients.
  • * on the right side represents convolution with respect to u.
  • q(u; ⁇ ) is the Fourier transform of the chirp exp( ⁇ i ⁇ 2 ) on the spatial axis and is given by equation (12). This becomes a chirp on the frequency axis.
  • Equation (16) is established using Dirac's ⁇ function.
  • E[ ⁇ ] represents an expected value.
  • Q is a constant.
  • Equation (17) the peak value of the cross-correlation coefficient between P 0 and P 1 is expressed by Equation (17). Substituting each specific expected value into this equation (17), the peak value of the cross-correlation coefficient is given by equation (18).
  • FIGS. 8A and 8B are plots of r peak values obtained using equation (19).
  • the peak value r peak must be sufficiently large. Therefore, it is necessary to suppress the decrease in this peak value, and it is found that it is effective to reduce the pulse width as shown in FIG. 8B.
  • the observation spectrum the spectrum obtained from the observation data (hereinafter also referred to as the observation spectrum) is compared with the above reference spectrum, and by taking the correlation, the shift amount a of the Rayleigh frequency and Obtain an estimate of the slope b.
  • a reference spectrum is created by processing electric field data Y 0 (z, ⁇ ) before change.
  • the electric field data before change is represented by the following equation (21).
  • the electric field of the reference spectrum is expressed by Equation (22) using Equation (21).
  • Equation (24) is the Fourier transform of an infinite chirp, which is expressed as an infinite chirp on the frequency axis as shown in Equation (12).
  • this infinite-length chirp has an infinitely large vibration component, it is difficult to discretize and handle it in signal processing. In other words, aliasing occurs no matter how high the sampling rate is.
  • Equation (28) is the Fourier transform of the infinite-length chirp represented by Equation (28), which is represented by Equation (29) using the Fresnel integral function.
  • ⁇ represents a complex conjugate
  • Z( ⁇ ) is a Fresnel integral function (see formulas (32) to (34) below)
  • arguments X 1 and X 2 are represented by formulas ( 30), given by equation (31).
  • the Fresnel integral function is defined by the following equations (32) to (34) and changes as shown in FIG.
  • the spatial frequency bandwidth of s(z; ⁇ ) is l p ⁇ . If the maximum assumed absolute value of ⁇ is ⁇ max , then the maximum bandwidth is l p ⁇ max .
  • Detection method 1 the reference spectrum is obtained in the frequency domain.
  • a) Obtain the Fourier transform q lp (u; ⁇ ) of the chirp from equations (29)-(31) above. Note that Z( ⁇ ) is the Fresnel integral function defined by Equation (32).
  • b) Take the frequency domain convolution of the Fourier transform of the chirp with the reference data (see equation (35)).
  • c) Obtain the reference spectrum by the following equation (36).
  • Equation (5) the Rayleigh spectrum P 1 (z, u) is calculated using Equation (5).
  • Equation (3) Find the cross-correlation coefficient when u is shifted between P 0 (z, u) and P 1 (z, u) for each distance z, 0 ⁇ z ⁇ L f (see equation (37)) , the maximum value of the correlation (see equation (38)) is sufficiently large, the deviation there is taken as the estimated value of ⁇ (see equation (39)).
  • Equation (40) the estimated value of the Rayleigh frequency shift u R in the interval [z ⁇ l p /2, z+l p /2] of the spatial coordinates is obtained by equation (40). Note that if the spatial sampling points are sufficiently dense, the estimated value may be evaluated using equation (41). (5 ) For z where the maximum value of the correlation coefficient r max (z) is not sufficiently large, P 0 (z, u; A cross-correlation coefficient r(z, ⁇ ; ⁇ j ) when u is shifted between ⁇ j ) and P 1 (z, u ) is obtained by equation (42).
  • Equation (44) the estimated value of ⁇ is determined by the following equation (45).
  • Equation (46) the estimated value of the Rayleigh frequency shift in the interval [z ⁇ l p /2, z+l p /2] of the spatial coordinates is obtained by equation (46). Note that if the spatial sampling points are sufficiently dense, the estimate may be evaluated using equation (47).
  • Detection method 2 In the detection method 1 described above, the reference spectrum was obtained by convolution in the frequency domain, but in the detection method 2, the product is calculated in the space domain using the Fourier transform, and then the inverse Fourier transform is used to restore the frequency. Use This detection method 2 can shorten the calculation time.
  • u k k ⁇ u is the spatial frequency and ⁇ u is its step size.
  • FIGS. 13A and 13B The assumed strain distribution is shown in FIGS. 13A and 13B.
  • FIG. 13A is a diagram showing the strain magnitude distribution in the length direction of the optical fiber.
  • FIG. 13B is a diagram showing the distribution of changes in strain magnitude in the length direction of the optical fiber.
  • the horizontal axis indicates the position (unit: m) in the length direction of the optical fiber. From these figures, it can be seen that there is strain in two sections in the optical fiber, a uniform strain of 10 ⁇ in the section from 20.7 m to 25.9 m, and a constant slope of 30 ⁇ /m in the section from 31.0 m to 34.1 m.
  • Equation (54) the relationship between distortion and Rayleigh frequency shift is given by Equation (54) below.
  • C R ⁇ is the conversion coefficient from distortion to Rayleigh frequency shift and takes a value of the order of ⁇ 151 MHz/ ⁇ .
  • 14A and 14B show the distribution of the Rayleigh frequency shift in the length direction of the optical fiber when the Rayleigh frequency shift is obtained from the strain change using this value.
  • the horizontal axis indicates the position (unit: m) in the length direction of the optical fiber.
  • Pulse width 80 ns
  • pulse waveform rectangular
  • frequency sweep width 4 GHz
  • pulse power 100 mW
  • light source line width 0.3 MHz
  • number of transmitted pulses 16.
  • FIG. 15 shows the result of applying the conventional method of detecting the spectral shift without considering the tilt.
  • the section from 31.0 m to 34.1 m it is about -4.5 GHz/m, but it is not detected at all in this section. Since the spatial resolution is 0.5 m, the tilt detection limit expected from equation (20) is 1.28 GHz/m, exceeding this limit in the section from 31.0 m to 34.1 m.
  • FIG. 16A shows the detection result of the Rayleigh frequency shift. Detection is performed at all points except the four points (20.7 m, 25.9 m, 31.0 m, and 34.0 m) where the Rayleigh frequency shift changes stepwise.
  • FIG. 16B shows the estimation error of the Rayleigh frequency shift. All but two points (25.2m, 33.1m) are within the quantization error range. It should be noted that the quantization error is caused by the fact that the interval of the frequency of the Fourier spectrum is 25 MHz and the lag at the time of obtaining the cross-correlation becomes an integer multiple of 25 MHz.
  • FIG. 17A shows the estimation result of the slope of the Rayleigh frequency shift. It can be seen that the magnitude of the inclination can be estimated well except for one point (33.1 m).
  • FIG. 17B shows peak values and thresholds of cross-correlation coefficients.
  • FIG. 18A shows distortion estimates converted from Rayleigh frequency shifts to distortion according to equation (54).
  • FIG. 18B also shows the estimation error of the distortion estimate converted from the Rayleigh frequency shift to distortion according to equation (54).
  • the strain gradient is as large as 30 ⁇ /m (Rayleigh frequency shift gradient is ⁇ 4.5 GHz/m) by using this detection method. It can be seen that the Rayleigh frequency shift can be detected with a high degree of accuracy even when it has a value.
  • the introduction of the analysis method N1 makes it possible to detect the Rayleigh frequency shift even when strain or temperature changes locally. That is, by using the analysis method N1, the reference spectrum is deformed by chirps of various chirplates, and the shift amount including the slope can be measured by taking the cross-correlation between them and the observed spectrum.
  • the analysis method N1 it was clarified that not only the spectrum shifts but also deformation due to chirp occurs when strain or temperature change has a slope.
  • the analysis method N1 can realize a Rayleigh intensity pattern measuring apparatus that is effective for non-uniform strain distribution.
  • the method of scanning a large number of chirped light pulses in a stepwise manner solves the second problem, that is, a sufficient measurement range. is realized.
  • an LD with a narrow linewidth is generally required.
  • the prior art No. 6 TGD-DAS technology described above requires an LD with a linewidth of 100 Hz.
  • this line width condition is relaxed, and a line width of 300 kHz, which is compatible with commercially available LDs, can also be handled. So, below, these contents are demonstrated below.
  • a plurality of chirped pulses are used to drive the receiving system, and heterodyne Detection is used to obtain the Rayleigh backscatter signal.
  • the Rayleigh intensity pattern for cross-correlation is obtained by using the matched filter method to generate simulated data from multiple chirped pulses, from which a set of subbands is extracted and combined. (RIP). Therefore, the details up to the formation of the Rayleigh intensity pattern in the first embodiment will be specifically described below based on examples.
  • FIG. 19 shows strain distribution data used in this simulation. Basically, it is similar to the strain distribution described in FIG. 14 above.
  • a strain distribution model having a constant strain of 10 ⁇ in a section of 20 m to 25 m and a non-uniform strain with a gradient of 30 ⁇ /m in the spatial resolution space in a section of 30 m to 33 m is used. Note that there is no distortion in sections other than the above (distortion is zero).
  • FIG. 20 shows main parameters used in this simulation.
  • FIG. 21 specifically shows a part of the multi-chirped pulse used in this simulation.
  • the chirped optical pulses are set so that the frequencies overlap (the overlapping frequency band between adjacent pulses is 0.175 GHz), and the frequency step is set between multiple pulses.
  • the laser line width in this embodiment is 300 kHz and the pulse width is 80 ns. Note that the laser linewidth is one of the most important parameters because it induces phase noise in the Rayleigh backscattered optical signal.
  • a matched filter for subband extraction is defined in the form shown in FIG.
  • a single-pulse chirp signal sweeps a frequency band of 4 GHz and is represented by equation (55) below.
  • f0 is the starting frequency
  • ⁇ is the chirp rate
  • tc is the pulse width
  • rect( ⁇ ) is the rectangular function.
  • a trace showing a particular frequency along the distance can be calculated by performing a cross-correlation of the acquired data with a matched filter. This procedure is called subband extraction, and as a result, a total of 2448 traces can be extracted from all simulation data (see Figure 23). The center frequency of the subband is then used to represent the frequency of the corresponding trace and is calculated by equation (56) below.
  • the Rayleigh backscatter spectrum consists of 2448 traces and covers frequencies from 0.1 GHz to 61.275 GHz as calculated by equation (56).
  • RIP is defined as the intensity of the Rayleigh backscatter signal at the corresponding laser frequency, but since the laser frequency before chirp is 193548.487 GHz, the absolute frequency of RIP is expressed as 193548.587 GHz to 193609.762 GHz. .
  • the RIP can be stored in absolute frequency. As long as the condition of the optical fiber does not change, the RIP can be obtained by sweeping the same absolute frequency even if the measuring instrument changes. Examples of RIPs at a distance of 10 m are shown in Figures 24A, 24B and 24C. Here, FIG.
  • FIG. 24A shows the RIP when the laser line width is 2 kHz
  • FIG. 24B shows the RIP when the laser line width is 100 kHz
  • FIG. 24C shows the RIP when the laser line width is 300 kHz. From these figures, it was found that the influence of the laser linewidth can be almost ignored in the linewidth range of 2 kHz to 300 kHz.
  • a sufficient measuring range can be realized by the method of stepwise scanning with multiple chirped light pulses. It was also shown that even with a laser line width of 300 kHz, a RIP that can be used for measurement can be obtained.
  • high-speed and high-precision strain or temperature measurement simultaneously realizes the three features of high spatial resolution, high-speed measurement, and a sufficient measurement range. It turned out that it is possible to provide a measuring device for In addition, it is possible to manage RIP information for a long period of several tens of years or more, and highly accurate RIP information can be continuously obtained during the measurement period. Moreover, even if a non-uniform strain distribution occurs in a section below the spatial resolution, the frequency shift can be obtained using the detection optimization analysis method.
  • Embodiment 2 the Rayleigh intensity pattern measuring apparatus according to the second embodiment will be described below with reference to the drawings, focusing on the differences from the first embodiment.
  • the Rayleigh intensity pattern measuring apparatus according to the second embodiment is different from the above description only in that the analysis method N2, which is a different analysis method from the analysis method N1 for eliminating non-uniform distortion, is used. Therefore, the analysis method N2 will be described below.
  • ACD approximate component decoding
  • the ACD method is roughly composed of the following four steps.
  • (d) Correlation coefficient calculation and determination step (a step of obtaining a cross-correlation coefficient based on the reconstructed data and calculating a correlation coefficient when the obtained correlation coefficient is equal to or greater than a threshold value)
  • the analysis method N2 is called the ACD method because the original data is decomposed into an approximation data portion and a detailed data portion using wavelet transform, which is used in the step combining the above steps (b) and (c). After that, the detailed data part is set to zero, and the data in the approximate data part is inversely wavelet-transformed to reconstruct the data.
  • an appropriate RFS can be calculated even when a correlation coefficient is obtained based on signal data of Rayleigh scattered light, which is detected as uneven distortion in the conventional method. Then, based on this calculated RFS, the target strain can be detected. Therefore, next, the detailed content of each of the above steps will be described in order with reference to the drawings.
  • the frequency spectrum data can be regarded as a waveform containing many frequency components, but a non-uniform strain distribution produces high frequency components that distort the waveform, making it impossible to obtain an appropriate correlation. Therefore, in order to remove the high-frequency components that are the factors that distort the waveform, the above signal data is decomposed using the discrete wavelet transform (hereinbelow, discrete wavelet transform is also referred to as DWT. Here, DWT: Discrete wavelet transform). do.
  • DWT discrete wavelet transform
  • FIG. 28 is a data processing flowchart in the ACD method showing the above processing as a whole.
  • two Rayleigh spectrum data are obtained (step S22).
  • the data at the decomposition level n of each data is decomposed by DWT, the detailed data portion of each data is set to zero, and stored in the memory database (step S23).
  • the Rayleigh spectrum data of the previously stored data are reconstructed by inverse wavelet transform (step S24).
  • the cross-correlation coefficient of the two reconstructed data is calculated (step S25).
  • step S26 The value of the calculated cross-correlation coefficient is compared with a threshold value (for example, 0.4) (step S26), and if the value is greater than the threshold value, the data processing ends. On the other hand, if the value of the calculated cross-correlation coefficient is smaller than the threshold, the value of n is increased by 1 (step S27), the process returns to step S23 of the above process flow, and the processes from step S23 to step S27 are calculated. Continue until the value of the cross-correlation coefficient is greater than the threshold.
  • a threshold value for example, 0.4
  • FIG. 29A is a graph showing the strain distribution produced in the optical fiber corresponding to the optical fiber position (see the distance shown on the horizontal axis) obtained using this ACD method.
  • a threshold value of 0.4 is used here to obtain the strain distribution (see FIG. 28).
  • FIG. 29B shows the distribution of correlation coefficients corresponding to optical fiber positions, obtained using this ACD method.
  • FIG. 29C shows a strain distribution graph obtained by the conventional method shown in FIG. 4A
  • FIG. 29D shows a corresponding correlation coefficient distribution.
  • FIG. 29A to 29D From these figures (FIGS. 29A to 29D), it can be seen that the ACD method is effective even when there is non-uniform strain.
  • FIG. 29B and FIG. 29D a value of 0.3 empirically used as a threshold value is indicated by a dotted line.
  • FIG. 29A In order to investigate the validity of the ACD method from another point of view, next, the strain distribution shown in FIG. 29A was compared with the data obtained using the Brillouin scattering method. The resulting data are shown in FIG. 30A. Similarly, FIG. 30B shows a graph obtained by smoothing the two types of data shown in FIG. 30A by the moving average method.
  • the strain data processed by the above ACD method based on Rayleigh scattered light measures the strain distribution equally compared to the data obtained using the Brillouin scattering method. judged to be done.
  • this ACD method is applied to a place with non-uniform distortion, the number of points that were successfully treated and the number of all treated points (total of points that were successfully treated and points that were not treated well) number) is defined as a success rate (unit: %). From the above, it can be seen that the ACD method is effective.
  • FCD method As the analysis method N2, a Fundamental Components Decoding (FCD) method, which is an analysis method similar to the ACD method, will be described below. That is, instead of the ACD method, the FCD method may be used as the analysis method N2. The contents of the FCD method will be described below, focusing on the differences from the ACD method.
  • FCD Fundamental Components Decoding
  • the ACD method uses discrete wavelet transform (DWT) as a basic analysis method, but the FCD method uses FFT (Fast Fourier Transform).
  • DWT discrete wavelet transform
  • FCD method uses FFT (Fast Fourier Transform).
  • LPF low-pass filter
  • the ACD method requires selection of the wavelet function to be used, but the FCD method does not.
  • the ACD method does not require an initial value
  • the FCD method requires an initial value. Therefore, the FCD method will be described in detail below in order from the above principle.
  • the object to be measured has an uneven strain distribution, its Rayleigh frequency spectrum will have high-frequency components, which will lead to distortion of the distorted waveform.
  • the two Rayleigh frequency data as shown in an example in FIG. (see the frequency components indicated by the arrows in ).
  • a method of obtaining the correlation is based on these common basic components, especially the basic components in the low frequency region.
  • the LPF is used here to cut the high frequency components of the Rayleigh frequency signal (in this case, it is necessary to select the optimum cutoff frequency).
  • the lower the cutoff frequency the higher the correlation coefficient between the two Rayleigh frequency spectra, and the higher the cutoff frequency, the lower the correlation coefficient between the two Rayleigh frequency spectra.
  • the optimum cutoff frequency find the appropriate value
  • FIG. 32 shows the data processing flow of the processing contents described above. The data processing flow in the FCD method will be described below with reference to FIG.
  • FIG. 32 shows a processing flow for calculating a series of Rayleigh frequency shifts with correlation coefficients higher than a given threshold A. Then, those Rayleigh frequency shifts are compared with the Rayleigh frequency shifts of the previous location (position number: n ⁇ 1), and the one with the smallest absolute value of difference is selected, and the current location (position number: n) is selected. be the Rayleigh frequency shift of .
  • two Rayleigh spectrum data (data 1, data 2) are passed through an LPF (low-pass filter) (step S31), and data 1 after passing the filter and data 1 after passing the filter Data 2 are obtained respectively, and the correlation between these two data is obtained (step S32).
  • the correlation coefficient obtained by taking the correlation is compared with a predetermined threshold value A (step S33).
  • the cutoff frequency of the LPF is changed (step S34), and the two Rayleigh spectrum data (data 1 and data 2) are again filtered by an LPF (low-pass filter). , and the correlation between the data 1 and the data 2 after passing through the filter is taken again and compared with the threshold value A.
  • FIG. 32 two Rayleigh spectrum data (data 1, data 2) are passed through an LPF (low-pass filter) (step S31), and data 1 after passing the filter and data 1 after passing the filter Data 2 are obtained respectively, and the correlation between these two data is obtained (step S32).
  • the correlation coefficient obtained by taking the correlation is compared with a predetermined threshold value A (step S33).
  • the value of the Rayleigh frequency shift is obtained from the correlation coefficient at this time (step S35) and stored in a memory or the like as Rayleigh frequency shift data. Then, the absolute value of the difference between the Rayleigh frequency shift data of the previous location obtained in advance and the Rayleigh frequency shift data obtained by the processing of steps S31 to S35 is obtained, and the value with the minimum value is obtained. It is selected (step S36) as the Rayleigh frequency shift of the current location (position number: n).
  • the first location is the location of the known Rayleigh frequency shift.
  • FIG. 33A shows the correlation coefficient obtained through the processing by this FCD method.
  • FIG. 33B shows data obtained without processing by the FCD method. It can be seen that the correlation coefficient is clearly improved in the data processed by the FCD method.
  • FIGS. 34A and 34B strain distribution data at distance positions including the distances showing non-uniform strain was obtained and shown in FIGS. 34A and 34B.
  • FIG. 34A the strain distribution data obtained by the processing by the FCD method are plotted as they are.
  • FIG. 34B shows data obtained by moving average and smoothing the data of FIG. 34A. It can be seen that both data are significantly improved by using the FCD method compared to the strain distribution data previously shown in FIG. 4A. This is similar to the ACD method described above.
  • the analysis method N2 has the same effects as the analysis method N1 described in the first embodiment. Moreover, in the analysis method N2, it was explained that both the ACD method and the FCD method are effective. Which method is selected can be determined according to the purpose of use, the availability of proprietary tools (for example, FFT in the FCD method) used for each method, and the like. Even when this analysis method N2 is used, the same effect as in the first embodiment can be obtained.

Abstract

レイリー強度パターン計測装置（１００）は、広帯域型の波長可変ＬＤ（１）、発振されたＬＤ光を光ファイバに入射すると共に光ファイバからのレイリー散乱光をＬＤ光の入射経路とは別の経路で出射する光カップラ（４）、波長可変ＬＤ（１）からのコヒーレント光とレイリー散乱光とを受信する受信部（１３）、受信部（１３）からの出力信号をＡＤ変換器（６）で受信し、位相情報を含めたデータを基に、レイリー散乱光から得た異なる２つのレイリー強度パターン信号から相互相関係数を演算し、与えられた閾値と比較した結果から得た相互相関係数を保存するＲＩＰデジタル処理部（１４）、を備え、比較した相互相関係数が閾値より小さい場合には、閾値以上となるまで相互相関係数を計算して、閾値以上となった相互相関係数からレイリー周波数シフトを求めることにより、被測定体の歪分布、あるいは温度分布を計測する。

Description

レイリー強度パターン計測装置およびレイリー強度パターン計測方法
　本願は、レイリー強度パターン計測装置およびレイリー強度パターン計測方法に関わる。
　分布型光ファイバセンシング（以下では、ＤＦＯＳ：Distributed Fiber Optic Sensingとも呼ぶ）技術の特徴として、その実用性においては、敷設した光ファイバがどのような状況においても、有効な計測値が提供できること、すなわち、信憑性が高いことにある。この分布型光ファイバセンシング技術の代表的な方式として、ラマン（Raman）、ブリルアン（Brillouin）、レイリー（Rayleigh）の３方式が挙げられる。
　特に、土木工事においては、このＤＦＯＳ技術は、施工品質の監視、および工事完成後の維持管理に至るまで、一貫してその適用が期待されている。しかしながら、未だ、実際に１００年程度の寿命を検証され、製造されているものはない。
　例えば、土木工事におけるＤＡＳ、あるいは他社の先行技術では、動的な変化が、高精度で計測できているが、光ファイバによる計測を基にしたものではなく、電気信号を基にしたものであるため、再現性に乏しく長期管理の対象に用いるには不向きである。
　一方、光ファイバを用いる方式は再現性もあり、長期管理の対象としてよりふさわしい方式と言える。このうち、上記ＤＦＯＳ技術の代表的な方式としては、上述の３方式が挙げられるが、これらのうち、特に、ブリルアン方式とレイリー方式が有力である。ただし、ブリルアン方式によるものは、絶対周波数での計測はできているが、その精度は、多くの場合、要求される精度に対して不十分なものとなっており、現状、レイリー方式が最も有効と想定される。
　このレイリー方式は、最近では、高い安定性、信号強度、マルチモード、シングルモードの双方での使用可能、高精度、１ｍｍ以下の高分解能という多くの特長を持つため、注目されている。すなわち、レイリー方式は、上記他の２つの方式と比較すると、（１）計測安定性、信号強度に優れており、シングルモード、マルチモードのいずれでも使用可能であること、（２）他の２つの方式に比べて、計測精度が高いこと、（３）１ｍｍ以下の高分解能を持つこと、の３つの優れた特長を持っている。
　以上、一般的な説明を行ってきたが、次に、上記ＤＦＯＳに関する先行技術を詳細に比較するため、図１に示した具体例を基に以下説明する。
ここで、図１は、縦方向に記載した６つの先行技術のそれぞれについての代表的な性能を、横方向に展開して記載した６つの指標に区分して、一覧表として示したものである（特許文献１、および非特許文献１－５参照）。
　図１において、横方向に示した、上記６つの指標である、絶対周波数（周波数の制御が可能で周波数の誤差が明確であるＬＤの周波数を意味する。以下同様）の有無、十分な計測範囲（計測範囲が十分である）の有無、不均一ひずみ分布の有効性（不均一ひずみ分布があっても、計測上、問題を生じないこと、言い換えれば、比較的小さな歪範囲においては相関不良の問題はないと判断されること）の有無、速度（すなわち、計測のスピード）、空間分解能、測定距離（測定可能な最大距離）を左から順に示している。以上において、計測範囲とは、論文、あるいは製品の値として公表された歪の計測範囲をいう。
　具体的数値例としては、図１において、先行技術番号１に示したＴＷ－ＣＯＴＤＲ（Tunable Wavelength Coherent Optical Time Domain Reflectometry）技術では、絶対周波数を２．３ＭＨｚの計測精度で制御しているが、部品依存のケースバイケース方式であり、次機種に継続する基準はない。また、４番目の行に示したＯＢＲ／Ｌｕｎａ技術では、不均一ひずみ分布の有効性については、分解能として数１０μｍまで可能となったため、この値まで歪の不均一性は問題にならないため、有効性は「有り」となっている。また、６番目の行に示したＴＧＤ－ＤＡＳ技術でも、位相変化を対象とすれば、そのまま歪に換算できるため、有効性は「有り」である。
　この図１より、測定距離が１０ｋｍ以上で速度と空間分解能の性能が水準以上あり、かつ、絶対周波数、十分な計測範囲、不均一ひずみ分布の有効性を全て満たすものは、この表に挙げた６つの先行技術には無いことが分かる。
　なお、図２は、ＤＦＯＳに関する上記６つの先行技術について、これら６つの先行技術において用いられている特徴的な技術を５種類に分類し、比較して示した図である（特許文献１、および非特許文献１－５参照）。この図２において、上記５種類の特徴となる技術のうち、周波数制御については、周波数制御手段（例えば、Etalon、Gas Cellなど）を有しているか否かを示している。
　ところで、本ＤＦＯＳ技術の実際のニーズとして、監視対象である橋、あるいは油井などにおいては、１００年の寿命が要求されている。また、地震などの計測においては、１００年以上の変化量を確実に捉える必要がある。
　そこで、この厳しい要求に対して、上述のレイリー方式に該当する、レイリー後方散乱光スペクトル（Rayleigh Backscattering Spectrum、以下、略してＲＢＳとも呼ぶ）方式が適用できないか、以下、検討する。レイリー後方散乱光スペクトル（ＲＢＳ）方式を管理対象の計測に適用する場合、この方式が以下に述べる３つの特長を持っているため、上述の要求を満たすことができると期待されるためである。
　すなわち、物質座標（運動する物体中に埋め込まれた座標）として光ファイバを利用することができ、位置を示す（場所を示す）ものとして認知可能であること、光信号でなくても検知可能であること、距離分解能、および計測精度に実用上の制限がなく使用できること、の３つの特長である。
　なお、管理対象の計測に適用されるレイリー後方散乱光スペクトル（ＲＢＳ）方式は、位置および周波数を変数とし、電場の強度（以下、電場は電界ともいう）、その絶対値、あるいは絶対値の２乗を関数値とするレイリー強度パターン（Rayleigh Intensity Pattern、以下、ＲＩＰと呼ぶこともある）と本質的には同じものであるので、本願では、レイリー強度パターン計測装置と称して、以下説明を行う。
　上記ＲＩＰを管理対象の計測手法として適用する場合には、空間位置を合わせて、時刻ｔ１でのＲＩＰ値と時刻ｔ２でのＲＩＰ値との相互相関を求め、この相関ピーク位置からレイリー周波数シフト量を求める。すなわち、パワーの相関でレイリー周波数シフトを計測するものである。
　ここで、ＲＩＰの計測技術について、図３を用いてさらに詳しく説明する。図３は、フルバンド可変波長分布帰還型ＬＤ（ＬＤ：Laser Diode）で取得したレイリー後方散乱光スペクトル方式によって得たＲＩＰ分布の一例を示した図である。
ここで、ＲＩＰの計測量としては、式（ＦＳ）で表されるレイリー周波数シフト量ΔνＲを求める（例えば、特許文献２参照）。
　　ΔνＲ＝Ｃ２１Δε＋Ｃ２２ΔＴ＋Ｃ２３ΔＰ・・・（ＦＳ）
なお、式（ＦＳ）において、Ｃ２１、Ｃ２２、Ｃ２３は、レイリー周波数シフトに与える感度係数であり、Δε、ΔＴ、ΔＰは、それぞれ、測定位置における光ファイバの歪、温度変化、圧力変化を示す。
　上述の歪変化、温度変化、および圧力変化によって導出された周波数シフトは、高度相関分析によって得られる（よって、周波数シフトは相対値である）。相関分析に使用される信号は、特定の距離範囲のパワースペクトルである。そして、位相に影響されない。
　従って、上記周波数シフトは、同じ場所における２つの状態（参照時と観測時）でのＴＷ－ＣＯＴＤＲ方式の周波数による明確な特徴の比較と考えられる。ここで、ＴＷ－ＣＯＴＤＲ方式と、高精度のブリルアン技術であるＰＰＰ－ＢＯＴＤＡ（Pulse Pre-Pump Brillouin Optical Time Domain Analysis）方式を復号して、レイリー後方散乱光方式およびブリルアン後方散乱光方式の２方式を使うハイブリッド計測法の要求を調和させることができる。
　一方、時間領域におけるレイリー後方散乱光方式は、ブリルアン後方散乱光方式に比較して、以下の利点がある。すなわち、シングルエンド方式のみ要求されること、長距離におけるクリープ監視などの小さな変化を正確に、かつ容易に検出可能であること、位相に対して独立した光ファイバからの信号であるため安定性と高い再現性があること、である。
米国特許第７４４００８７号明細書 国際公開２０１４／０８３９８９号 特許第６５５２９８３号公報
Y.Koyamada,"Novel Technique to Improve Spatial Resolution in Brillouin Optical Time-Domain Reflectometry",IEEE PHOTONICS TECHNOLOGY LETTERS, VOL.19, NO.23, DECEMBER 1,2007 Y.Wang et al.,"Distributed Fiber-Optic Dynamic-Strain Sensor With Sub-Meter Spatial Resolution and Single-Shot Measurement",IEEE Photonic Journal, Vol.11, No.6, DECEMBER, 2019 J.Pastor-Graells, et al.,"Single-shot distributed temperature and strain tracking using direct detection phase-sensitive OTDR with chirped pulses", OPTICS EXPRESS 13121, Vol.24, No.12, June 13, 2016 S.Liehr, et al.,"Wavelength-Scanning Distributed Acoustic Sensing for Structual Monitoring and Seismic Application", www.mdpi.com/journal/proceedings, Proceedings 2019,15,30 D.Chen, et al.,"Fiber-optic distributed acoustic sensor based on a chirped pulse and a non-matched filter", OPTICS EXPRESS 29415, Vol.27, No.20, September 30, 2019
　しかしながら、ＴＷ－ＣＯＴＤＲ方式を含む先行技術には、以下に挙げる課題がある。
　まず、図１に示した最初の指標である「絶対周波数（絶対基準の周波数）」が未確立であること。先行技術番号２から先行技術番号６に示した、いずれの先行技術も、図１に示したように、ＬＤ周波数、あるいは波長の値を制御して誤差を明確にできていない。また、先行技術番号１に示したＴＷ－ＣＯＴＤＲ方式では、ＬＤの波長（光周波数）を制御してはいるが、絶対値の管理を行ってはいない。従って、（同一メーカーの）異なる機器間、あるいは異なるメーカーの機器間における計測においては、相対的な周波数シフトの値しか得ることができない。
　次に、２番目の指標である「計測範囲」が十分でないこと（実用化に十分な計測範囲が確保できていないこと）が挙げられる。先行技術番号１以外の各技術では、計測範囲が小さ過ぎる。具体的には、例えば、１００℃の温度変化はよくあるケースであるが、チャープパルス技術を用いる各手法（先行技術番号１のＴＷ－ＣＯＴＤＲ方式以外の各先行技術）においては、要求される数値の４桁以下の範囲でしか計測できていない。
　次に、３番目の指標である「不均一ひずみ分布に対する有効性」が欠如していることが挙げられる。先行技術番号１のＴＷ－ＣＯＴＤＲ方式では、空間分解能の範囲内において、言い換えると、パルス幅の半分以内で不均一なひずみ分布が発生すると、適切なひずみの検出はできなくなる。空間分解能の長さ範囲内の不均一なひずみ分布によって引き起こされるレイリー強度パターン（ＲＩＰ）の歪のために、参照データとの相関が失われるためである。
　また、長期に亘る監視中に、一度、相関関係がとれなくなると、初期の値との関連性が失われるため、長期の監視ができなくなる。これは、レイリー後方散乱光方式が、原理的に、時間軸方向でデータを蓄積することにより、最終的に変形を計測しているためである。
　相関がとれなかった場合の実測例について、図４を用いて説明する。図４Ａは、棒材を用いてレイリー後方散乱光の周波数シフト（以下、ＲＦＳとも呼ぶ）を以下に示す測定パラメータで測定した場合の測定結果を示している。
　測定パラメータ；掃引開始周波：１９２０１０ＧＨｚ、掃引期間：４００ＧＨｚ、掃引ステップ：０．２ＧＨｚ、空間分解能：２ｃｍ、サンプリング間隔：１ｃｍ。
　図４Ａにおいて、点線で示した長丸で囲まれた領域（横軸に示した距離位置で５１２．３ｍから５１３．１ｍの範囲内）では、ひずみが正しく検出されていない。具体的には、例えば、距離位置５１２．５ｍでは、２つのデータ、すなわち、参照データと実測データとの相関係数は、図４Ｂに示すように、０．２以下となっており、相関がないことを示している。
　従来、レイリー後方散乱光方式では、時間方向にデータを蓄積して変形を最終的に求めるようにしているため、データ蓄積中に１個のデータの蓄積が旨く行かなかった場合には、それ以降の計測が、初期値（初期の歪変化量）との相関を失ってしまうからである。
　上記のＲＦＳデータについて検討するため、２つのデータの相関関係を利用しない測定法である、ブリルアン散乱光を基にして測定した棒材の測定データとの比較を行う。図５は、上記と同じ棒材を用いてブリルアン散乱光を基にして棒材に発生したひずみを以下に示す測定パラメータで測定した場合の測定結果を示したデータである。
　測定パラメータ；掃引開始周波：１０ＧＨｚ、掃引期間：１．４５ＧＨｚ、掃引ステップ：５ＭＨｚ、空間分解能：５ｃｍ、サンプリング間隔：２．５ｃｍ。
　図５は、ブリルアン後方散乱光方式を用いて測定した棒材の歪データに関し、横軸に距離（単位：ｍ）を採った場合のひずみの分布を示している。図５から、距離５１２．２ｍ付近で１００με以上であった歪の値が、距離５１２．５ｍ付近で－４００μεに変化し、差し引き５００με以上の歪変化が測定されていることが判る。
　この図５を用いて、距離５１２．２ｍ付近での歪変化を調べると、図中に記載したように、距離変化４．１ｃｍに対して１５９．８４μεの歪変化が発生していることが判る。言い換えると、距離５１２．５ｍ付近での４．１ｃｍ内（４ｃｍは、レイリー後方散乱光方式で用いられるパルス長に相当）の領域において、１ｍ当たりの歪変化量は、３８９８．５４μεであった。
　以上に説明したように、空間分解能の範囲内で急激なひずみ変化がある場合（不均一ひずみがある場合）には、同じ棒材の歪変化を、レイリー後方散乱光方式で測定した場合においては、正確な歪の変化、あるいは相互相関が測定できなくなることが分かった。
　なお、通常、データを補間することにより相関は改善されるが、今回のケースでは、相関は改善されなかった。なお、データ処理時の周波数（変化）の単位ステップは０．０４ＧＨｚである。
　次に、４番目の課題として、ホモダイン受信を受信方式としている先行技術（表の１番目と４番目の技術。図２参照）においては、高度化信号解析に必要な位相情報が、現行のホモダイン受信では、最初から失われるという問題があることである。
　さらに、５番目の課題として、表で指標とした速度が遅い場合には、高速現象に追随できないだけでなく、計測中に歪変化があれば、レイリー後方散乱光の周波数シフトの結果中に、異なる周波数シフト（ΔνＲ）が含まれているために、先に説明したように、相関が取れない。このため、やがて、全体として、監視の有効性が失われることに繋がることとなる。
　次に、図１に示した各先行技術を実現する際の、技術的な課題について説明する。ＴＷ－ＣＯＴＤＲ（先行技術番号１）については、全帯域のＬＤを使用するため、計測範囲としては十分広いが、ＬＤの線幅が広すぎるという問題点がある。例えば、４ＭＨｚの場合においては、有効なコヒーレント長は、数１０ｍしかなく、位相解析のＴＧＤ（Time Gated Digital）において、位相信号に大きな位相ノイズを含むため、そのままでは使えない。すなわち、ＴＷ－ＣＯＴＤＲ方式では、図２に示したように受信方式にはホモダイン受信を採用しているため、位相情報が最初から失われることがあり、高度化信号解析に必要な位相情報が利用できない。
　また、ＴＧＤ－ＤＡＳ方式（図１の先行技術番号６に示した方式）を用いた計測では、線幅が１ＫＨｚ以下の光源を使用して旨く測定できていたが、この場合においては、周波数の変更がほとんどできない（市販品がほぼ無いため）という問題がある。
　さらに、広い計測範囲を確保するためには、ＲＩＰの受信周波数の範囲（受信帯域幅）は、１００ＧＨｚ以上が必要となるが、電気信号による受信帯域としてカバーするのは不可能に近いというのが現状である。
　以上に述べた課題を、代表的な先行技術ごとに具体的な数値例を挙げてまとめると以下のようになる。先行技術番号１のＴＷ-ＣＯＴＤＲ方式は、ひずみ、温度などの静的計測に優れ、高空間分解能（２ｃｍ）、高精度な計測手法であることに特徴があるが、計測の速度については、計測時間が数１０秒と遅いことが欠点である。一方、先行技術番号６のＴＧＤ－ＤＡＳ方式は、音波計測など動的計測に優れ、計測の速度が毎秒２０００回以上と高速であることに特徴があるが、空間分解能は１ｍ程度であり、劣る。また、図１には記載していないが、微小な変化計測に限定され、例えば、温度変化は０．２℃までであり、大変化する計測には向かない。
　本願は、上記のような課題を解決するためになされたものであって、高速・高精度な計測を実現する上でトレードオフの関係にある３つのパラメータ、空間分解能、計測の速度、および十分な計測範囲を同時に実現でき、かつ、長期間にわたり取得したデータの管理が可能となるレイリー強度パターン計測装置を提供することを目的とする。
　本願に開示されるレイリー強度パターン計測装置は、
波長可変ＬＤと、前記波長可変ＬＤから発振されるレーザ光の周波数を変更可能な制御器と、を有する光源部、
前記光源部から発振されたレーザ光を光ファイバに入射するとともに、前記光ファイバからのレイリー散乱光を前記レーザ光の入射経路とは別の経路で出射する光カップラ、
前記光源部からのローカル光と、前記光カップラからの前記レイリー散乱光と、をコヒーレント受信する受信部、
前記受信部から出力された信号をＡＤ変換するＡＤ変換器と、前記ＡＤ変換器で変換された信号を所定の方法で演算する第１の演算部と、前記第１の演算部で演算された信号が入力され保存されるメモリ・データベースと、前記メモリ・データベースに保存されたデータを基に所定の演算を行う第２の演算部と、を有し、前記レイリー散乱光の電界信号から得られる、予め定められた測定位置での異なる２つのレイリー強度パターンから、前記第２の演算部で前記測定位置について所定の補正をした後、相互相関係数を求め、求めた相互相関係数が予め定められた閾値以上となった場合の相互相関係数を前記メモリ・データベースに保存するレイリー強度パターンデジタル処理部、
を備え、
予め定められた閾値以上となった相互相関係数を基にレイリー周波数シフトを求めることにより、被測定体の歪分布、あるいは温度分布を定めることを特徴とするものである。
　本願に開示されるレイリー強度パターン計測装置によれば、高速・高精度な計測を実現する上でトレードオフの関係にある３つのパラメータ、空間分解能、計測の速度、および十分な計測範囲を同時に実現でき、かつ、長期間にわたり取得したデータの管理が可能となるレイリー強度パターン計測装置を提供することを目的とする。
実施の形態１のレイリー強度パターン計測装置の課題となる技術指標を説明するための図である。 ＤＦＯＳの特徴的な技術を説明するための図である。 レイリー後方散乱光方式によって得たＲＩＰ分布の一例を示す図である。 レイリー後方散乱光方式によって測定した、棒材の歪と相関関数の例を示す図である。 ブリルアン後方散乱光方式によって測定した、棒材の歪データを示す図である。 本実施の形態１に係るＲＩＰ計測装置の全体構成を説明するための図である。 本実施の形態１に係るＲＩＰ計測装置の動作を説明するための図である。 本実施の形態１に係るＲＩＰ計測装置の不均一歪の検出適正化解析法において、レイリー周波数シフトの相互相関係数のピーク値を説明するための図である。 不均一歪の検出適正化解析法において用いられるフレネル積分を説明するための図である。 不均一歪の歪変化が０．１ＧＨｚ／ｍの場合の有限チャープ波形とそのフーリエ変換を説明するための図である。 不均一歪の歪変化が５ＧＨｚ／ｍの場合の有限チャープ波形とそのフーリエ変換を説明するための図である。 不均一歪の歪変化が２０ＧＨｚ／ｍの場合の有限チャープ波形とそのフーリエ変換を説明するための図である。 不均一歪の検出適正化解析法における、歪モデルと歪の傾きモデルを示す図である。 不均一歪の検出適正化の解析法における、レイリー周波数シフトモデルとその傾きモデルを示す図である。 不均一歪の検出において、レイリー周波数シフトの傾きを考慮しないモデルでのシミュレーション結果を説明するための図である。 不均一歪の検出適正化の解析法を適用してレイリー周波数シフトの分布を求めた場合のシミュレーション結果を説明するための図である。 不均一歪の検出適正化の解析法を適用してレイリー周波数シフトの傾きの分布を求めた場合のシミュレーション結果を説明するための図である。 不均一歪の検出適正化の解析法を適用して歪の分布を求めた場合のシミュレーション結果を説明するための図である。 本実施の形態１におけるに係るＲＩＰ計測装置において、実施例の説明に用いるレイリー強度パターンを形成するための歪分布データを説明するための図である。 レイリー強度パターン計測装置において、受信部を駆動するために使用される、レーザ光を模擬する複数のチャープパルスの走査方法を説明する際に用いるパラメータを説明するための図である。 複数のチャープパルスをステップ状に走査する方法について説明するためのモデル図である。 各チャープパルスからサブバンドを抽出するための整合フィルタの形態を説明するための図である。 各チャープパルスから被測定体の距離に沿った特定の周波数を示すトレースを抽出するまでの一連の方法を説明するための図である。 レーザ線幅に分けて、対応するレーザ周波数でのレイリー後方散乱信号の強度であるＲＩＰの例を示した図である。 不均一歪の検出適正化の別の解析法におけるデータ準備工程を説明するための図である。 不均一歪の検出適正化の別の解析法におけるデータ分解工程を説明するための図である。 不均一歪の検出適正化の別の解析法における相関係数計算および決定工程を説明するための図である。 不均一歪の検出適正化の別の解析法におけるデータ処理フローチャートを示した図である。 不均一歪の検出適正化の別の解析法を用いて、不均一な歪が検知されたデータを処理した結果を示した図である。 図２９の補足説明図である。 不均一歪の検出適正化の別の解析法に類似した手法に用いる、基本周波数成分を説明するための図である。 不均一歪の検出適正化の別の解析法に類似した手法に用いられる、データ処理フローを示す図である。 不均一歪の検出適正化の別の解析法に類似した手法を適用した結果の一例である相関係数を示した図である。 不均一歪の検出適正化の別の解析法に類似した手法を適用した結果の一例である歪分布を示した図である。 不均一歪の検出適正化の別の解析法に類似した手法の効果を説明するための図である。
実施の形態１．
　実施の形態１に係るレイリー強度パターン計測装置（以下、簡略化してＲＩＰ計測装置とも呼ぶ）の全体構成について、以下、図６、図７を用いて説明する。
　図６は、本実施の形態１に係るＲＩＰ計測装置１００の全体構成を説明するための図である。実施の形態１に係るＲＩＰ計測装置１００は、チャープパルスを分割して用いるＴＷ（可変波長）－ＴＧＤ－ＲＦＳ（レイリー周波数シフト）方式の計測装置である。
　本実施の形態１に係るＲＩＰ計測装置１００は、大別して４つの主要な構成要素を備えている。それらは、光源部１１、送受信光学部１２、受信部１３、ＲＩＰ確定デジタル処理部１４（以下、レイリー強度パターンデジタル処理部１４、あるいはＲＩＰデジタル処理部１４とも呼ぶ）である。以下、これらの詳細について順に説明する。
　光源部１１は、広帯域型の波長可変レーザ（ＴＷ－ＬＤ）であるＬＤ１とＬＤ１から出射されるレーザ光を制御する制御器を有し、自動温度制御および自動周波数制御を行いつつ、内蔵する光学局部発振器（ＯＬＯ）により、ＬＤ１の絶対周波数制御（ここで、絶対周波数制御とは、例えば、Gas Cellのように、分子の吸収波長を基準とする。その情報は公知であり、温度依存性も低く、計測精度に影響しない範囲で簡単に実現できる）、および周波数ステップ制御を行って、分割されたチャープパルスを出力として発振する。なお、波長可変レーザには外部変調型も含む。
　すなわち、光源部１１では、チャープ光パルスを利用することにより、波長可変のＬＤ１の周波数を制御し、この周波数を複数、ステップ状に変更しつつ（例えば、２００ＧＨｚで１２ステップに変更するなど）掃引する。この場合において、上記ＴＷ―ＬＤを活用して、ＴＧＤ（Time Gated Digital）法等を利用し、高精度光パルスを生成する。なお、ＴＧＤ法における掃引幅は、例えば、４ＧＨｚ（サンプルレートは８ＧＳ／ｓ）である。
　送受信光学部１２は、レーザ信号をパルス化する光スイッチ２、この光スイッチ２の出力信号を増幅するエルビウム添加光ファイバ増幅器３（Erbium-Doped Fiber Amplifier。以降、略称してＥＤＦＡ３と呼ぶ）、ＥＤＦＡ３からの信号を、参照用光ファイバ５を介して計測対象体（以下、被測定体（具体的には光ファイバ）とも呼ぶ）に伝送するとともに、伝送された信号の計測対象体からの散乱光である後方散乱光を、参照用光ファイバ５を介して受信後、後述の外部の受信器に送出する光カップラ４を有している。
　受信部１３は、上記送受信光学部１２の光カップラ４から送出された後方散乱光を受信するとともに、上記光源部１１からの光源光を受信する受信器を有する。この受信器は、受信感度の改善と位相情報の取得のため、受信の際に、偏波ダイバーシティ方式で受信するヘテロダイン受信器、あるいは位相ダイバーシティ方式・偏波ダイバーシティ方式で受信するホモダイン受信器で構成されている。
　ＲＩＰデジタル処理部１４は、デジタル手法で信号の強度、あるいは位相解析をする部分である。ここでは、上記受信部１３からの出力信号を高速に２ｃｈＡＤ変換器６でデジタルに変換して受信し、ＦＰＧＡ（Field Programmable Gate Array）、ＡＳＩＣ（Application Specific Integrated Circuit）などを用いて、後述する４種類の方法で信号を高速に処理する。また、ＲＩＰデジタル処理部１４は、デジタル信号を演算処理するプロセッサと、このプロセッサで処理された信号を記憶して保存するメモリ・データベース８（以下メモリ・ＤＢ８とも称する）を持ち、２ｃｈＡＤ変換器６から受信した信号を演算して（例えば、複数のチャープ信号を合成する演算を行う）、受信した信号とは異なる信号としてメモリ・データベース８に出力する第１の演算部７と、このメモリ・ＤＢ８に記憶された情報から、相関解析する２つの時刻（例えば、時刻ｔｉと時刻ｔｊ）のデータを選択し、これら異なる２つのデータを基に、初期のレイリー周波数シフトΔνＲ‐Iを求め、周波数軸での補間処理などのリサンプリング、距離補正、相関（相関判断）などの演算処理を行って、確定したレイリー周波数シフトΔνＲ‐Ｄとして出力する第２の演算部９（この第２の演算部９については、図７などを用いて以下で詳述する。以下、この第２の演算部９を単に演算部９と呼ぶこともある）、を有する。
　ここで、２ｃｈＡＤ変換器６は、例えば、信号の振幅成分と位相成分を、各ｃｈに、例えばＦＰＧＡチップなどの処理器で変換することが可能である。また、ＲＩＰデジタル処理部１４では、デジタル手法で信号の強度と位相解析を行うことができる。なお、図においては、上記確定したレイリー周波数シフトΔνＲ‐Ｄの出力が、ＲＩＰデジタル処理部１４の内部に含まれる態様になっているが、ＲＩＰデジタル処理部１４の外部にあっても良い。
　なお、上述の４つの主要な構成要素のうち、実線の枠で囲んだ、光源部１１、受信部１３、ＲＩＰデジタル処理部１４の３つは、特に、本実施の形態１に係るＲＩＰ計測装置の特徴となる構成要素である。一方、点線の枠で囲んだ送受信光学部１２については、既存のＯＴＤＲ技術を適用することが可能であり、特に目新しいものではない（例えば、特許文献３参照）。
　また、図中、方法１、方法２、方法３、方法４として示したものは、受信部１３での受信方法に応じて、ＲＩＰデジタル処理部１４での動作の違いを明確に説明するため、それぞれ、別の名称を与えたものであり、これらの詳細については、以下、図７を用いて詳しく説明する。
　そこで、以下、特に上記の第２の演算部９を中心として、本実施の形態１に係るＲＩＰ計測装置全体の動作について、図７を用いて上記方法１、方法２、方法３、方法４ごとに詳しく説明する。いずれの方法を用いた場合でも、信号の受信感度の改善と高密度波長多重化が可能となる。
　まず、方法１では、受信部１３でヘテロダイン受信した信号であるチャープ信号を、２ｃｈＡＤ変換器６では特別な処理は行わず（当該信号をスルーさせて）、そのまま、ＲＩＰデジタル処理部１４のメモリ・ＤＢ８に記憶する（ステップＳ１）。このメモリ・ＤＢ８では時系列ｔｋでデータが記憶されているため、このデータの中から、例えば、時刻ｔｉ、および時刻ｔｊの２つの時刻でのデータを選択する（ステップＳ４）。この選択した２つの時刻でのデータを基に、第２の演算部９において、以下の一連の動作を実行する。
　最初に、初期のレイリー周波数シフトΔνＲ－Ｉを演算で求める（ステップＳ５）。次に求めた初期のレイリー周波数シフトΔνＲ－Ｉについて、距離補正のために補間処理を含む周波数リサンプリングの処理を行う（ステップＳ６）。次に、周波数リサンプリングしたデータを基に、距離補正後の相互相関係数を求める（ステップＳ７）。次に、求めた相互相関係数の値が予め定めた閾値より大きいか否かを判断する（ステップＳ８）。この閾値は経験に基づき、例えば、０．３として予め定めることができる。
　上記判断の結果が「Ｙｅｓ」であれば、求めた相互相関係数の値を基に、レイリー周波数シフトΔνＲを求め、求めたレイリー周波数シフトΔνＲの値を確定したレイリー周波数シフトΔνＲ‐Ｄとして定める（ステップＳ１０）。
　一方、上記判断の結果が「Ｎｏ」であれば、不均一歪の検出適正化のための解析を行う。すなわち、距離補正後の相互相関係数を再計算するための基となるデータを解析により求める（ステップＳ９）、再度、距離補正後の相互相関係数を求める（ステップＳ７）。次に、求めた相互相関係数の値が予め定めた閾値より大きいか否かを判断する（ステップＳ８）。
　そして、この相互相関係数の値が予め定めた閾値より大きくなるまで、上記ステップＳ９→上記ステップＳ７→上記ステップＳ８の動作を繰り返し、相互相関係数の値が予め定めた閾値より大きくなった場合の相互相関係数の値を基にレイリー周波数シフトΔνＲを求め、求めたレイリー周波数シフトΔνＲの値を、確定したレイリー周波数シフトΔνＲ－Ｄとして定める（ステップＳ１０）。
　方法１では、２ｃｈＡＤ変換器６による処理がないため、処理データの高速処理が可能となる。
　次に、方法２について、方法１と異なる点を中心に説明する。方法２では、受信部１３でヘテロダイン受信した信号であるチャープ信号を、２ｃｈＡＤ変換器６でＡＤ変換し、ＦＰＧＡなどによりコムフィルタリング処理をした後の信号を、ＲＩＰデジタル処理部１４のメモリ・ＤＢ８に記憶する（ステップＳ２）。その後の処理は方法１と同じであるので、ここでは詳しい説明は省略する。
　方法２では、コムフィルタリング処理をすることで、信号に含まれる周期性ノイズを除去することができ、測定データの高精度化が図れる。
　次に、方法３について、方法１、方法２と異なる点を中心に説明する。方法３では、方法２の処理後、さらにこの処理した信号を２乗処理したデジタル処理後の信号を、ＲＩＰデジタル処理部１４のメモリ・ＤＢ８に記憶する（ステップＳ３）。その後の処理は方法１と同じであるので、ここでは詳しい説明は省略する。
　方法３では２乗した信号を用いるため、信号成分のうち、位相の影響を受けることが少なくなり、その分、測定データの高精度化が図れる。
　最後に、方法４について、方法１、方法２、方法３と異なる点を中心に説明する。方法１から方法３では、受信部１３でヘテロダイン受信した信号であるチャープ信号を用いたが、方法４では、これら３つの方法とは異なり、受信部１３で、位相情報を含む信号であるホモダイン受信した信号を用いる。また、位相情報を失わないようにするために、位相ダイバーシティ方式を採用している。このホモダイン受信した信号を２乗処理した後の信号は、ＲＩＰデジタル処理部１４のメモリ・ＤＢ８に記憶される（ステップＳ３）。上記以降の動作は方法１と同じなので、ここでは詳細の記載は省略する。
　方法４は、方法１から方法３のヘテロダイン受信信号を使用する場合に比べ、位相雑音の少ない狭線幅のＬＤと光ＰＬＬが必要となるが、最も高感度の受信が可能になり、また、高速でのデータ処理がより容易であるというメリットを持つ。
　実施の形態１に係るレイリー強度パターン計測装置１００は、以上のように構成されているので、実効的にパルスパワーを向上でき、また、複数の周波数における散乱波形を同時に取得できることから、送信する光パルスの数を大幅に減らすことができ、測定時間を短縮することが可能となる。なお、ＲＩＰの相関を解析するため、ＲＩＰデジタル処理部１４では大量の解析が行われる。
　上記方法１から方法４では、図７に示したような、検出適正化の解析法を用いて、適正に不均一歪の検出が可能になるよう工夫されている。そこで、次に、上記方法１から方法４で共通に使用される、上記不均一歪を適正に検出できる解析法について、以下詳しく説明する。この解析法には大きく分けて解析法Ｎ１と解析法Ｎ２の２種類がある。解析法Ｎ１から順に、以下説明する。
[解析法Ｎ１]
　解析法Ｎ１は、不均一歪がある場合でも、レイリー散乱光の参照データの電界から得られる参照スペクトル及びそれを変形して得られる種々の参照スペクトルと、観測データの電界から得られる観測スペクトルとの相互相関から、被測定体のひずみ分布を適正に検出することができる解析法である。参照スペクトルの変形は、歪の傾き（後ほど詳しく説明する）に対処するためであり、傾きの大きさにしたがって変形も異なる。以下、この解析法Ｎ１に用いられる処理方法（アルゴリズム）について詳しく説明する。
　上記アルゴリズムを作成するにあたり、まず、基礎となるモデルにより、パラメータを定める。ここでは、光ファイバの入力端から周波数ν、パルス幅Ｄのレーザ光を入射し、散乱されて戻ってくる光を観測した場合、特に、レイリー散乱光Ｐを観測した場合を考える。
　このような場合に、観測値として得られるのは、直接検波で得られる光の強度、あるいはヘテロダイン受信で得られる電界の強度であり、いずれも、Ｐ＝Ｐ（ｔ、ν）と、時間ｔと周波数νの２つのパラメータによって表される。例えば、レーザ光源の周波数νをスキャンしながら観測値を取得することで、レイリー散乱光のスペクトルデータが得られる。
　そこで、次に、光ファイバ中の歪、あるいは温度が変化した場合を考える。この場合、光路長が変化する。つまり、光ファイバの入力端から距離ｚまでの光の往復時間τ（ｚ）が変化する。この場合の観測値は、上記の時間ｔを光の往復時間τ（ｚ）で置き換えて、Ｐ＝Ｐ（τ（ｚ）、ν）と表すことができる。このＰ（τ（ｚ）、ν）は、記法を簡略化し、以下では、Ｐ（ｚ、ν）と記載する。
　そこで、以下で、歪、あるいは温度の変化前後のスペクトルをまず検討するため、歪、あるいは温度の変化前のスペクトルをＰ０（ｚ、ν）と記し、変化後のスペクトルをＰ１（ｚ、ν）と記す。
　いま、距離区間[ｚ－ｌｐ／２、ｚ＋ｌｐ／２]（以下、距離区間は、単に、区間と呼ぶ）において、歪、あるいは温度の変化が一定であると仮定し、レイリー周波数シフトνＲ（ｚ）と記すと、レイリー周波数シフトνＲ（ｚ）は、この区間で一定の値となる。
　そうすると、変化後のスペクトルＰ１（ｚ、ν）と変化前のスペクトルＰ０（ｚ、ν）との間には、以下の関係式が成立することが知られている。
Figure JPOXMLDOC01-appb-M000001
 つまり、変化後のスペクトルは変化前のスペクトルの周波数シフトで表される。
［歪、あるいは温度の変化がパルス内で一定でない場合］
　そこで、次に、課題となっている、歪、あるいは温度の変化がパルス内で一定でない場合を考える。この場合には、スペクトル間の関係は、上記式（１）の周波数シフトでは表されなくなり、検出性能の劣化をもたらすため、以下詳細に検討する。
［歪、あるいは温度の変化がパルス内で一定でない場合］
　計算上の都合により、周波数変数を上述の時間軸での周波数νから、対応する空間周波数（式（２）参照）に変更する。また、レイリー周波数シフトも、下記の式（３）とする。
Figure JPOXMLDOC01-appb-M000002
Figure JPOXMLDOC01-appb-M000003
　ここで、観測される、変化前と変化後の散乱光の電界をそれぞれ、Ｙ０（ｚ、ｕ）、Ｙ１（ｚ、ｕ）と置くと、これらの絶対値の２乗がスペクトルを与える（式（４）、式（５）参照）。
Figure JPOXMLDOC01-appb-M000004
Figure JPOXMLDOC01-appb-M000005
　ここで、Ｙ０（ｚ、ｕ）、Ｙ１（ｚ、ｕ）は、それぞれ、式（６）、式（７）で与えられる。
Figure JPOXMLDOC01-appb-M000006
Figure JPOXMLDOC01-appb-M000007
 なお、電界データを用いるためには、ヘテロダイン検波を行う必要がある。
　パルス幅に相当する区間内で、歪、あるいは温度の変化が線形に推移するとし、レイリー周波数シフトは、区間[ｚ－ｌｐ／２、ｚ＋ｌｐ／２]で、以下に示す式（８）となることを仮定する。
Figure JPOXMLDOC01-appb-M000008
 ここで、ａは区間の中心でのレイリー周波数シフトであり、ｂは区間内でのレイリー周波数シフトの傾きである。ここで、パラメータａ、ｂを、以下に示す空間周波数でのパラメータ表現α、βに変更する。すなわち、α＝２ａ／ｖｇ、β＝２ｂ／ｖｇとすると、レイリー周波数は、α、βを用いて、式（９）となる。
Figure JPOXMLDOC01-appb-M000009
　ここで、式（９）を活用すると、式（７）より次の式（１０）が成り立つ。
Figure JPOXMLDOC01-appb-M000010
　式（１０）の右辺は、２つの関数ρ（ｚ＋ζ）、－ｌｐ／２≦ζ≦ｌｐ／２とｅｘｐ（πｉβζ２）、－∞＜ζ＜∞の積のζに関するフーリエ変換になっている。ここで、ρ（ｚ）は、距離ｚにおけるレイリー後方散乱係数であり、２つの関数の積のフーリエ変換は、それぞれの関数のフーリエ変換のコンボリューション（畳み込み）に等しい。また、ρ（ｚ＋ζ）、－ｌｐ／２≦ζ≦ｌｐ／２のフーリエ変換は、式（６）からＹ０（ｚ、ｕ）に等しい。よって、定数係数を除いて以下の式（１１）が成り立つ。
Figure JPOXMLDOC01-appb-M000011
 ここで、右辺の＊はｕに関するコンボリューションを表す。
　ｑ（ｕ；β）は、空間軸上でのチャープｅｘｐ（πｉβζ２）のフーリエ変換であり、式（１２）で与えられる。
Figure JPOXMLDOC01-appb-M000012
 これは周波数軸上でのチャープになる。
　従って、変化後のスペクトルは、式（１３）で表せるが、区間内でｕＲが一定になる場合と異なり、元のスペクトルを用いて表現することができない。
Figure JPOXMLDOC01-appb-M000013
　式（４）、式（５）、式（６）、式（１０）において、簡単のためにα＝０とし、ｚは任意に一つに固定し、引数であるｚの表記を省略すると、散乱光のスペクトルは式（１４）、式（１５）のように簡単な表記で表される。
Figure JPOXMLDOC01-appb-M000014
Figure JPOXMLDOC01-appb-M000015
　ここで、レイリー後方散乱係数ρ（ｚ）が、空間的には白色ガウス過程とみなせるので、ディラックのδ関数を利用すると、以下の式（１６）が成立する。ここでＥ［・］は、期待値を表す。また、Ｑは定数である。
Figure JPOXMLDOC01-appb-M000016
　式（１６）を利用して、Ｐ０とＰ１に関する２次までの期待値を求めると、Ｐ０とＰ１の相互相関係数のピーク値は、式（１７）で表される。
Figure JPOXMLDOC01-appb-M000017
 この式（１７）に具体的な各期待値を代入すると、相互相関係数のピーク値は、式（１８）となる。
Figure JPOXMLDOC01-appb-M000018
　さらに、パルス幅ＤがＤ＝２ｌｐ／ｖｇ、パルス内でのレイリー周波数シフトの変動量（最大と最小の差）が、ΔνＲ＝ｂｌｐ＝ｖｇβｌｐ／２となることから、ＤΔνＲ＝β（ｌｐ）２
となるので、相互相関係数のピーク値は結局、式（１９）となる。
Figure JPOXMLDOC01-appb-M000019
　図８Ａ、図８Ｂは、式（１９）を用いて求めたｒｐｅａｋの値をプロットしたものである。図８Ａから分かるように、ｒｐｅａｋはＤΔνＲ＝０のとき１であり、ＤΔνＲの増加とともに、急速に減少して０（ゼロ）に近づく。スペクトル同士の相互相関でレイリー周波数シフトの大きさを検出する方法では、ピーク値ｒｐｅａｋが十分大きくなる必要がある。従って、このピーク値の低下を抑えることが必要であり、そのためには、図８Ｂに示したように、パルス幅を小さくすることが効果的であることが判る（実際にはパルス幅を小さくすると、比例してΔνＲも小さくなるので、両方の効果でピーク値の低下が抑制される）。
　そこで、例えば、ｒｐｅａｋ＝０．３のときのＤΔνＲの値を求めると３．１１となる（図８Ａ参照）。この値を適用して、ＤΔνＲ≦３．１１が適用限界の目安となる。レイリー周波数シフトの傾きに直すと、以下に示す式（２０）で表される。
Figure JPOXMLDOC01-appb-M000020
 例えば、空間分解能ｌｐ＝０．５ｍの場合には、ｂ≦１．２８ＧＨｚ／ｍである。
［歪、あるいは温度の変化が区間内で傾きを持つ場合のレイリー周波数シフトの検出法］
　歪、あるいは温度の変化が区間内で一定であればスペクトルのシフトとして変化を検出できるが、傾きを持っていれば、上述のように検出の限界がある。しかしながら、スペクトルの光の強度を求める前の電界データが利用できれば、傾きまで含めた検出が可能になる。以下その方式について説明する。
　すなわち、一言でいうと、チャープによるスペクトル変形を利用する方式である。まず、参照データが得られた段階で、チャープによるスペクトル変形でレイリー周波数シフトに種々の傾きがある場合のスペクトルを参照スペクトルとして保存する。次に、観測データが得られた段階で、その観測データから得られたスペクトル（以下、観測スペクトルとも言う）を先の参照スペクトルと比較し、相関を取ることによって、レイリー周波数のシフト量ａと傾きｂの推定値を得る。
［チャープによるスペクトル変形］
　各距離ｚ（０≦ｚ≦Ｌｆ、Ｌｆ：ファイバ長）で周波数νを走査したときの電界データＹ（ｚ、νｋ）、ｋ＝０、１、・・・、Ｋ－１が与えられるものとする。νｋ＝ｋΔνであり、Δνは周波数間隔である。Ｋは、周波数方向のデータの個数で、ＫΔνがデータの周波数帯域幅になる。レイリー周波数シフトの傾きｂも離散化して、ｂｊ、ｊ＝０、１、・・・、Ｊ－１とする。
　参照スペクトルは、変化前の電界データＹ０（ｚ、ν）を加工して作成する。ここでは、計算の都合上、周波数νを空間周波数ｕ＝２ν／ｖｇ（単位：ｃｙｃｌｅ／ｍ）へ変換し、レイリー周波数シフトの傾きｂをβ＝２ｂ／ｖｇ（単位：ｃｙｃｌｅ／ｍ２）へ変換する。そうすると、変化前の電界データは次の式（２１）で表される。
Figure JPOXMLDOC01-appb-M000021
 また、式（２１）を用いて参照スペクトルの電界は式（２２）で表される。
Figure JPOXMLDOC01-appb-M000022
　そこで次に、上式（２２）の右辺が、関数の積のフーリエ変換になっていることに着目して、上述の［歪、あるいは温度の変化がパルス内で一定でない場合］においては、フーリエ変換同士のコンボリューションとして、以下に示す式（２３）から式（２５）のように表現した。
Figure JPOXMLDOC01-appb-M000023
Figure JPOXMLDOC01-appb-M000024
Figure JPOXMLDOC01-appb-M000025
　式（２４）のｑ（ｕｋ、βｊ）は無限長のチャープのフーリエ変換であり、式（１２）に示したように、周波軸上での無限長のチャープとして表される。しかしながら、この無限長のチャープは、無限に大きい振動成分を持つことになるから、信号処理において離散化して扱うのが困難である。つまり、サンプリングレートをいくら大きくしてもエイリアシングが生じてしまう。
　そこで、空間分解能がｌｐであることを利用して、式（２３）、式（２４）の代わりに、等価な次の式（２６）、式（２７）を、それぞれ用いる。
Figure JPOXMLDOC01-appb-M000026
Figure JPOXMLDOC01-appb-M000027
　ここで、ｑｌｐ（ｕ；β）は、式（２８）で表される無限長のチャープのフーリエ変換であり、フレネル積分関数を用いて式（２９）で表される。
Figure JPOXMLDOC01-appb-M000028
Figure JPOXMLDOC01-appb-M000029
 ここで、｛￣｝は、複素共役を表し、Ｚ（・）はフレネル積分関数（下記の式（３２）～式（３４）参照）であり、引数Ｘ１、Ｘ２は、それぞれ、式（３０）、式（３１）で与えられる。
Figure JPOXMLDOC01-appb-M000030
Figure JPOXMLDOC01-appb-M000031
　上記フレネル積分関数は、以下に示す式（３２）から式（３４）で定義され、図９に示すような変化をする。
Figure JPOXMLDOC01-appb-M000032
Figure JPOXMLDOC01-appb-M000033
Figure JPOXMLDOC01-appb-M000034
 　ｘ→±∞のとき、Ｃ（ｘ）、Ｓ（ｘ）→±０．５となるが、これから、ｌｐ→∞のとき、ｑｌｐ（ｕ；β）→ｑ（ｕ；β）となることがわかる。
　計算に必要なｑｌｐ（ｕ；β）のｕ方向への長さは、ｓ（ｚ；β）の帯域幅で異なる。図１０Ａ、図１０Ｂ、図１１Ａ、図１１Ｂ、図１２Ａ、図１２Ｂには、それぞれ、β＝１、５０、２００ｃｙｃｌｅ／ｍ２（ｂ＝０．１、５、２０ＧＨｚ／ｍに相当）の場合の有限チャープ波形ｓ（ｚ；β）とそのフーリエ変換ｑｌｐ（ｕ；β）を示す。
　なお、ｌｐ＝０．５ｍとした。ｓ（ｚ；β）の空間周波数帯域幅は、ｌｐβである。想定されるβの絶対値の最大値をβｍａｘとおくと、帯域幅の最大は、ｌｐβｍａｘとなる。
［傾きまで含めたレイリー周波数シフトの検出法］
　次に、各距離ｚにおいて、距離方向への傾きまで含めたレイリー周波数の検出法を示す。この検出法においては、レイリー周波数シフトの傾きを考慮しなくても十分な相関が得られる場合には、傾きの検出は行わない。ここでも、計算の都合上、周波数νの代わりに空間周波数ｕ＝２ν／ｖｇ、レイリー周波数シフトのパラメータａ、ｂの代わりに、α＝２ａ／ｖｇ、β＝２ｂ／ｖｇを用いて以下表す。これらの関係は、散乱光における距離ｚと往復の時間ｔのｔ＝２ｚ／ｖｇの関係と同じであり、同様の換算ができる。
［検出法１］
　検出法１では、参照スペクトルは周波数領域で求める。
（１）参照データの電界Ｙ０（ｚ、ｕｋ）、ｋ＝０、１、・・・、Ｋ－１が与えられるものとする。これから、各β＝βｊ、ｊ＝０、１、・・・、Ｊ－１について、参照スペクトルＰ０（ｚ、ｕｋ;βｊ）を以下に説明する方法で求めておく。
　ａ）チャープのフーリエ変換ｑｌｐ（ｕ；β）を、上記の式（２９）～式（３１）から求める。なおＺ（・）は式（３２）で定義されるフレネル積分関数である。
　ｂ）チャープのフーリエ変換と参照データとの周波数領域でのコンボリューションを取る（式（３５）参照）。
Figure JPOXMLDOC01-appb-M000035
 　ｃ）参照スペクトルを次式（３６）で求める。
Figure JPOXMLDOC01-appb-M000036
 （２）観測データの電界Ｙ１（ｚ、ｕ）、０≦ｚ≦Ｌｆ、ｕｓｔ≦ｕ≦ｕｅｎｄが得られる都度、レイリースペクトルＰ１（ｚ、ｕ）を、式（５）を用いて求める。
（３）各距離ｚ、０≦ｚ≦Ｌｆについて、Ｐ０（ｚ、ｕ）とＰ１（ｚ、ｕ）のｕをずらしたときの相互相関係数を求め（式（３７）参照）、相関の最大値（式（３８）参照）が十分大きければ、そこでのずれをαの推定値（式（３９）参照）とする。
Figure JPOXMLDOC01-appb-M000037
Figure JPOXMLDOC01-appb-M000038
Figure JPOXMLDOC01-appb-M000039
 （４）この場合、空間座標の区間[ｚ－ｌｐ／２、ｚ＋ｌｐ／２]におけるレイリー周波数シフトｕＲの推定値を式（４０）で求める。
Figure JPOXMLDOC01-appb-M000040
 なお、空間サンプリング点が十分密ならば、推定値は式（４１）で評価してよい。
Figure JPOXMLDOC01-appb-M000041
 （５）相関係数の最大値ｒｍａｘ（ｚ）が十分大きくならないｚに対しては、各βｊ、ｊ＝０、１、・・・、Ｊ－１について、Ｐ０（ｚ、ｕ;βｊ）とＰ１（ｚ、ｕ）のｕをずらしたときの相互相関係数ｒ（ｚ、α；βｊ）を式（４２）によって求める。
Figure JPOXMLDOC01-appb-M000042
 そして式（４３）で表される相関係数の最大値ｒｍａｘ（ｚ、βｊ）が最大となるβｊを求め、βの推定値とする（式（４４））。
Figure JPOXMLDOC01-appb-M000043
Figure JPOXMLDOC01-appb-M000044
 また、αの推定値は、次の式（４５）で定める。
Figure JPOXMLDOC01-appb-M000045
 （６）この場合、空間座標の区間[ｚ－ｌｐ／２、ｚ＋ｌｐ／２]におけるレイリー周波数シフトの推定値を式（４６）で求める。
Figure JPOXMLDOC01-appb-M000046
 なお、空間サンプリング点が十分密ならば、推定値は式（４７）で評価してよい。
Figure JPOXMLDOC01-appb-M000047
［検出法２］
　上述の検出法１では、参照スペクトルを周波数領域でのコンボリューションで求めたが、検出法２では、フーリエ変換を用いて空間領域で積の演算を行ってから、逆フーリエ変換で周波数に戻す方法を用いる。この検出法２の方が計算時間は短縮できる。検出法２では、検出法１の（１）の代わりに、下記（１ａ）を用いる。なお、検出法２の（２）～（６）については、上記検出法１の（２）～（６）と同じであるので、ここでは詳細な説明は省略する。
（１ａ）参照データの電界Ｙ０（ｚ、ｕｋ）、ｋ＝０、１、・・・、Ｋ－１が与えられるものとする。ｕｋ＝ｋΔｕは空間周波数でΔｕはその刻み幅である。空間変数の刻み幅はΔζ＝１／（ＫΔｕ）とし、空間変数は、ζｎ＝ｎΔζ、ｎ＝０、１、・・・、Ｋ－１とする。これから、各β＝βｊ、ｊ＝０、１、・・・、Ｊ－１について、参照スペクトルＰ０（ｚ、ｕｋ;βｊ）を以下に説明する方法で求めておく。
　ａ）各β＝βｊ、ｊ＝０、１、・・・、Ｊ－１について空間軸上のチャープｓ（ζｎ；βｊ）を求める。そのため、まず、－Ｋ／２≦ｎ≦Ｋ／２－１の範囲でｓ（ζｎ；βｊ）を次式（４８）に従って求める。
Figure JPOXMLDOC01-appb-M000048
 　その後、循環的に延長して、式（４９）とし、後の計算で、ｓ（ζｎ；βｊ）、ｎ＝０、１、・・・、Ｋ－１を用いる。
Figure JPOXMLDOC01-appb-M000049
 　ｂ）Ｙ０（ｚ、ｕｋ）、ｋ＝０、１、・・・、Ｋ－１をｕに関して逆フーリエ変換し、式（５０）を計算する。
Figure JPOXMLDOC01-appb-M000050
 　ｃ）各β＝βｊ、ｊ＝０、１、・・・、Ｊ－１について空間軸上でチャープｓ（ζｎ；βｊ）、ｎ＝０、１、・・・、Ｋ－１との積を取り、フーリエ変換する（式（５１）参照）。
Figure JPOXMLDOC01-appb-M000051
 　ｄ）各β＝βｊ、ｊ＝０、１、・・・、Ｊ－１について、参照スペクトルを式（５２）に従って求める。
Figure JPOXMLDOC01-appb-M000052
 特に、βｊ＝０の場合には、参照スペクトルを次に示す式（５３）で定める。
Figure JPOXMLDOC01-appb-M000053
［シミュレーションによる検証］
　ここでは、チャープ光パルスを用いて受信側で多数の周波数に分割する方式を用いてシミュレーションによる検証を行う。なお、シミュレーションに用いたパラメータは以下のとおりである。
　ａ）対象ファイバ；
　光ファイバの長さは５０ｍである。想定した歪分布を図１３Ａ、図１３Ｂに示す。図１３Ａは、光ファイバの長さ方向における歪の大きさの分布を示した図である。図１３Ｂは、光ファイバの長さ方向における歪の大きさの変化の分布を示した図である。なお、いずれの図においても、横軸は光ファイバの長さ方向の位置（単位：ｍ）を示している。
これらの図から、光ファイバ中の２つの区間で歪があり、２０．７ｍ～２５．９ｍの区間では１０μεの一様な歪、３１．０ｍ～３４．１ｍの区間では一定の傾き３０με／ｍを持った１０με～１００μεの歪が存在することがわかる。
　ここで、歪とレイリー周波数シフトの関係は、以下に示す式（５４）で与えられる。
Figure JPOXMLDOC01-appb-M000054
 　式（５４）において、ＣＲεは歪からレイリー周波数シフトへの変換係数であり、－１５１ＭＨｚ／με程度の値を取る。図１４Ａ、図１４Ｂには、この値を用いて、歪変化からレイリー周波数シフトを求めた場合における、光ファイバの長さ方向におけるレイリー周波数シフトの分布を示した。図１３Ａ、図１３Ｂと同様、いずれの図においても、横軸は光ファイバの長さ方向の位置（単位：ｍ）を示している。
　ｂ）チャープ光パルス；
　パルス幅：８０ｎｓ、パルス波形：矩形、周波数掃引幅：４ＧＨｚ、パルスパワー：１００ｍＷ、光源線幅：０．３ＭＨｚ、送信パルス数：１６。
　ｃ）周波数配置
　個々のチャープ光パルスの周波数掃引範囲：νｎ～νｎ＋４ＧＨｚ、パルス間の周波数間隔：νｎ－νｎ―１＝３．８２５ＧＨｚ、個々のチャープ光パルスで設定可能なサブバンド（幅２００ＭＨｚ、間隔：２５ＭＨｚ）の数：１５３、１６のチャープ光パルスで設定可能なサブバンドの数：２４４８。
　ｄ)信号処理
　チャープパルスによるレイリー散乱光を受信して、チャープのサブバンド毎に分割復調することで、フーリエスペクトルを求める（参照スペクトル、観測スペクトルの両者）。変形した参照スペクトルと観測スペクトルの相関を取ることによりレイリー周波数シフトを求める。なお、シミュレーションにおける傾きのパラメータはΔｂ＝０．５ＧＨｚ／ｍ（Δβ＝４．８７cycle／ｍ２）、ｂｊ＝－６０～６０ＧＨｚ／ｍ（βｊ＝－５８４．４～５８４．４ｃｙｃｌｅ／ｍ２）、ｊ＝ｊ＝０、１、・・・、Ｊ－１、Ｊ＝２４１とした。
　以下、シミュレーション結果について説明する。比較のため、まず、傾きを考慮しないでスペクトルシフトを検出する従来の方法を適用した結果を図１５に示す。３１．０ｍ～３４．１ｍの区間では、－４．５ＧＨｚ／ｍ程度であるが、この区間では全く検出できていない。空間分解能は０．５ｍであるから、式（２０）から予想される傾きの検出限界は１．２８ＧＨｚ／ｍであり、３１．０ｍ～３４．１ｍの区間ではこの限界を超えている。
　次に、傾きを考慮した本検出方法を用いた場合に、スペクトル変形により、レイリー周波数シフトを検出した結果を図１６Ａ、図１６Ｂ、図１７Ａ、図１７Ｂ、図１８Ａ、図１８Ｂに示す。図１６Ａに、レイリー周波数シフトの検出結果を示す。レイリー周波数シフトがステップ状に変化する４点（２０．７ｍ、２５．９ｍ、３１．０ｍ、３４．０ｍ）以外の全点で検出がなされている。、図１６Ｂにレイリー周波数シフトの推定誤差を示す。２点（２５．２ｍ、３３．１ｍ）を除いて、量子化誤差の範囲内に収まっている。なお、量子化誤差は、フーリエスペクトルの周波数の刻みが２５ＭＨｚで相互相関を取る際のラグが２５ＭＨｚの整数倍になることから生じるものである。
　図１７Ａに、レイリー周波数シフトの傾きの推定結果を示す。１点（３３．１ｍ）を除いて傾きの大きさを旨く推定できていることがわかる。図１７Ｂには、相互相関係数のピーク値と閾値を示す。
　図１８Ａには、式（５４）に従ってレイリー周波数シフトから歪に変換された歪推定値を示す。また、図１８Ｂには、式（５４）に従ってレイリー周波数シフトから歪に変換された歪推定値の推定誤差を示す。
　以上で説明したシミュレーション結果から、本検出法を用いることにより、空間分解能が０．５ｍの場合には、歪の傾きが３０με／ｍ（レイリー周波数シフトの傾きが－４．５ＧＨｚ／ｍ）といった大きな値を持つ場合でも、レイリー周波数シフトを高い精度で検出できていることがわかる。
　以上、解析法Ｎ１の導入により、歪、あるいは温度の変化が局所的に変化する場合においても、レイリー周波数シフトの検出が可能になることを示した。すなわち、解析法Ｎ１を用いることで、参照スペクトルを種々のチャープレートのチャープで変形させておき、それらと観測スペクトルとの相互相関を取ることで、傾きを含むシフト量を計測できることを示した。なお、上記解析法Ｎ１においては、歪、あるいは温度の変化が傾きを持つときには、スペクトルがシフトするだけではなく、チャープによる変形が生じることも明らかにした。
　また、具体的な結果として、０．１５ｍの空間分解能の場合との比較で、従来の検出法では検出できなかった大きな傾き、例えば５０ＧＨｚ／ｍの傾きを持つ場合でもシフト量の検出が可能になることを明らかにした。
　以上のように、解析法Ｎ１により、不均一ひずみ分布に対する有効性を有するレイリー強度パターン計測装置が実現できることを示した。
　さらに、実施の形態１のレイリー強度パターン計測装置においては、多数チャープ光パルスをステップ状で走査する方法（各ステップでパルス化したチャープ信号とする）により、第２の課題である十分な計測範囲を実現している。ただし、このチャープ光パルスを発生させるためには、一般的には、狭い線幅のＬＤが必要である。例えば、先に説明した先行技術番号６のＴＧＤ―ＤＡＳ技術では、線幅が１００ＨｚのＬＤが必要である。本実施の形態１のレイリー強度パターン計測装置においては、この線幅の条件を緩和し、市販品のＬＤでも対応可能な線幅である３００ｋＨｚについても対応可能である。そこで、以下では、これらの内容について以下に説明する。
　先に説明した図６、図７に示したように、実施の形態１のレイリー強度パターン計測装置においては、方法１から方法３では、複数のチャープパルスを使用して受信システムを駆動し、ヘテロダイン検出を使用してレイリー後方散乱信号を取得している。
　さらに詳しく言うと、整合フィルタ法を使用し、複数のチャープパルスからシミュレーションデータを作成し、これを基にサブバンドのセットを抽出し、それらを組み合わせることにより、相互相関を求めるためのレイリー強度パターン（ＲＩＰ）を形成している。
　そこで、以下、本実施の形態１におけるレイリー強度パターンを形成するまでの詳細について、実施例をもとに具体的に説明する。
［レイリー強度パターン形成の実施例］
　まず、本シミュレーションに用いるひずみ分布のデータを図１９に示す。基本的には上記図１４で説明したひずみ分布に類似のものである。ここでは、区間２０ｍ～２５ｍでは一定の歪１０με、３０ｍ～３３ｍでは空間分解能空間内に傾き３０με／ｍの不均一な歪を持つ歪分布をモデルとする。なお、上記以外の区間では歪はない（歪はゼロ）。
　次に、ひずみ分布の計測に用いるレーザ発振光を模擬する、１６個の多数チャープパルスをステップ状に走査する方法について、図２０、図２１を用いて説明する。図２０には本シミュレーションで用いた主なパラメータを示した。また、図２１には、本シミュレーションで用いた多数チャープパルスの一部の態様について、具体的に示した。
　図２０、図２１からわかるように、チャープ光パルスの間には周波数がオーバーラップするように設定（隣接するパルス間のオーバーラップ周波数帯域は０．１７５ＧＨｚ）して、多数パルス間に周波数ステップのないＲＩＰを生成するようにするとともに、１６個の多数チャープパルスをステップ状に走査することにより、後程詳しく説明するように、十分な計測範囲である周波数範囲０．１ＧＨｚ～６１．２７５ＧＨｚを実現した。
　ここで、図２０に示したように、本実施例におけるレーザ線幅は３００ｋＨｚ、パルス幅は８０ｎｓである。なお、レーザ線幅は、レイリー後方散乱光信号の位相ノイズを誘発するため、最重要なパラメータの１つである。
　次に、本実施例の整合フィルタについて図２２を用いて説明する。サブバンド抽出のための整合フィルタは、図２２に示したような形態で定義される。単一パルスのチャープ信号は、４ＧＨｚの周波数帯域を掃引し、以下の式（５５）で表される。
Figure JPOXMLDOC01-appb-M000055
 ここで、ｆ０は開始周波数、κはチャープレート、ｔｃはパルス幅、ｒｅｃｔ（・）は矩形関数を示す。
　図２０、図２２から分かるように、サブバンド幅は２００ＭＨｚ、また、隣接するバンド間の間隔は２５ＭＨｚとした。従って、各サブバンドの整合フィルタｓｎ（ｔ）、ｎ＝１、２、・・・、１５３は、以上述べた定義に従って決定される。
　これから、距離に沿った特定の周波数を示すトレースは、整合フィルタと取得データの相互相関を実行することによって計算できる。この手順をサブバンド抽出と呼び、結果として、全てのシミュレーションデータから合計２４４８個のトレースを抽出できる（図２３参照）。このとき、サブバンドの中心周波数は、対応するトレースの周波数を表すために使用され、次の式（５６）で計算される。
Figure JPOXMLDOC01-appb-M000056
　すなわち、レイリー後方散乱スペクトルは、２４４８トレースで構成され、式（５６）で計算された０．１ＧＨｚから６１．２７５ＧＨｚまでの周波数をカバーする。ＲＩＰは対応するレーザ周波数でのレイリー後方散乱信号の強度として定義されるが、チャープ前のレーザ周波数が１９３５４８．４８７ＧＨｚであるので、ＲＩＰの絶対周波数は、１９３５４８．５８７ＧＨｚ～１９３６０９．７６２ＧＨｚと表される。よって、ＲＩＰは、絶対周波数で保存可能である。光ファイバの状態が変化しない限り、計測機器が変わっても、同じ絶対周波数を掃引させることにより、ＲＩＰを取得できる。距離１０ｍにおけるＲＩＰの例を図２４Ａ、図２４Ｂ、図２４Ｃに示す。ここで、図２４Ａにはレーザ線幅２ｋＨｚの場合のＲＩＰを、図２４Ｂにはレーザ線幅１００ｋＨｚの場合のＲＩＰを、図２４Ｃにはレーザ線幅３００ｋＨｚの場合のＲＩＰを、それぞれ示す。これらの図より、レーザ線幅の影響は、線幅２ｋＨｚ～３００ｋＨｚでは、ほとんど無視できることがわかった。
　以上説明したように、本実施の形態１のレイリー強度パターン計測装置においては、多数チャープ光パルスをステップ状で走査する方法により、十分な計測範囲を実現することが可能となった。また、レーザ線幅３００ｋＨｚでも、計測に使用可能なＲＩＰを得ることができることを示した。
　以上のように、本実施の形態１のレイリー強度パターン計測装置よれば、高い空間分解能、高速の計測、および十分な計測範囲の３つの特長を同時に実現した、高速・高精度の歪、あるいは温度の計測装置を提供することができることがわかった。
　また、数１０年以上の長年に渡り、ＲＩＰ情報の管理が可能となり、計測期間中は、高精度のＲＩＰ情報を連続的に取得することができる。
　また、空間分解能以下の区間に不均一な歪分布が生じても検出適正化解析法を用いて、周波数シフトを求めることができる。
　上記に加え、さらに、光ファイバを用いて、絶対周波数を保証したＲＩＰを計測するに際し、絶対周波数の精度はＬＤの線幅程度の保証ができ、計測したＲＩＰの範囲は実用上、十分な温度・歪変化の範囲を含み、空間分布上、不均一に生じた相関不良を是正する手段を備えることで、初期に計測したＲＩＰについて継続して信頼性のある相関値を得ることができ、動的歪変化にも対応できる高速での計測が可能となる。
　また、例えば３００ｋＨｚ程度の線幅が不十分なＬＤを用いても、動的歪を計測することが可能である。
実施の形態２．
　次に、実施の形態２に係るレイリー強度パターン計測装置について、実施の形態１と異なる点を中心に、以下、図を用いて説明する。実施の形態２に係るレイリー強度パターン計測装置では、上述した説明のうち、不均一歪を解消する解析法Ｎ１とは異なる解析方法である解析法Ｎ２を用いる点のみが異なる。そこで以下、この解析法Ｎ２について説明する。
［解析法Ｎ２］
　解析法Ｎ２においては、不均一歪が存在するような場合であっても、適切な相関が得られるようにするため、以下に説明する離散ウェーブレット変換手法に基づく近似成分復号法（以下、ＡＣＤ法とも呼ぶ。ＡＣＤ：Approximation Components Decoding）を用いる。以下、この方法の詳細内容と、レイリー周波数シフト（以下、簡略化してＲＦＳとも呼ぶ）の測定データを用いて、この方法の検証を行った例について説明する。
　上記ＡＣＤ法は、概略、以下に示す４つの工程で構成されている。
（ａ）データ準備工程（前処理工程。解読するデータを取得する工程）
（ｂ）データ分解工程（離散ウェーブレット変換により、データを所定の分解レベルで近似データ部と詳細データ部に分解する工程）
（ｃ）データ再構築工程（上記詳細データ部をゼロとおいて、ウェーブレット逆変換により、データを再構築する工程）
（ｄ）相関係数計算および決定工程（上記再構築データを基に相互相関係数を求め、求めた相関係数が閾値以上となる場合の相関係数を算出する工程）
　解析法Ｎ２をＡＣＤ法と呼ぶのは、特に、上記工程（ｂ）と（ｃ）とを併せた工程で使用される、ウェーブレット変換を用いて元のデータを近似データ部と詳細データ部に分解した後、詳細データ部をゼロとおき、近似データ部のみのデータをウェーブレット逆変換してデータの再構築を行う処理があることによる。
　このＡＣＤ法を用いれば、従来法では不均一な歪として検知されたレイリー散乱光の信号データを基にして相関係数を求めた場合であっても、適切なＲＦＳが算出できる。そして、この算出したＲＦＳを基に、目的とする歪を検知することができる。そこで、次に、上記各工程の詳細内容について、以下順に図を用いて説明する。
［データ準備工程］
　先に図４において点線の長丸で囲んだ領域中の不均一地点で得た生データから選択した、所定距離（例えば、５１２．５ｍ）における２つの周波数スペクトルデータ（データ１、データ２）（図２５Ａ、図２５Ｂ参照）を取得し、これらのデータを正規化した後、さらに小さい歪を検出する必要がある場合には、周波数ステップに関して複数倍の補間処理を行う。例えば、周波数ステップを０．２ＧＨｚから０．０４ＧＨｚに変更すると、検出できる最小歪は１．３２μεが０．２６４μεとなり、より小さな値まで検出可能となる。そして、その後、相互相関係数（図２５Ｃ参照）を求める。ここで、例えば、データ１が参照データ、データ２が計測データである。ここで、図２５Ｃに示したように、相互相関係数は、経験的に定めた閾値（０．３）以下の値を示しており、検知が旨く行われていないことがわかる。
［データ分解工程］
　上記周波数スペクトルデータは、数多くの周波数成分を含んだ波形と見なせるが、不均一歪分布は、適切な相関が取れなくなる要因となる、波形をゆがめる高周波成分を生じさせる。そこで、この波形をゆがめる要因である高周波成分を取り除くため、以下に述べる離散ウェーブレット変換（以下、離散ウェーブレット変換をＤＷＴとも呼ぶ。ここで、ＤＷＴ：Discrete wavelet transform）を用いて、上記信号データを分解する。この手法は、信号処理において信号を分解する通常の方法である。
　すなわち、図２６に示すように、上記信号データは、段階ごとに対応して（図中、分解レベル１、２、３，・・・、に従う段階ごとに対応して）、低周波数部分である、個々の近似データ部Ａｎ（ｎ＝１、２、３、・・・）と、高周波数部分である個々の詳細データ部Ｄｎ（ｎ＝１、２、３、・・・）とに、分けることができる。
　さらに詳しくは、図２６に示すように、各近似データ部Ａｎ（ｎ＝１、２、３、・・・）は、レベルが上がるごとに順に、さらに高いレベルの近似データ部Ａｎ＋１（ｎ＝１、２、３、・・・）と、詳細データ部Ｄｎ（ｎ＝１、２、３、・・・）に分解することができる。この場合において、最大レベルは、信号の長さと選択したウェーブレット関数に依存する（ウェーブレット変換においては使用するウェーブレット関数を選択する必要がある。ここでは、正規直交変換関数であるｃｏｉｆ２を選択して使用した）。
［データ再構築工程］
　そこで、次に、上記データの分解により（この際、ウェーブレット関数ｃｏｉｆ２を使用する）、求めた詳細データ部であるＤ１、Ｄ２、Ｄ３、…、を０（零）とおくことにより、詳細データ部を容易に取り除くことができる。そこで、詳細データ部が取り除かれた、上記詳細データ部との干渉がより少ない近似データ部のみを用いて、ウェーブレット逆変換することにより、詳細データ部が取り除かれたデータを再構築する。
　以上のアイデアを分析し拡張することにより、提案するＡＣＤ法の主要な処理フローを構築する。
［相関係数計算および決定工程］
　次に、上記再構築されたデータ（２つのレイリースペクトルデータ）を基に、これらのデータの相互相関を取り、相関係数を求める。この際、図２６に示した分解レベルにより、求めた相関係数の結果が異なることに注意する必要がある（図２７参照）。そして、分解レベルが高いほど、信号データ成分がより多く削除されることにより、相関係数が高くなる。一方で、実際のケースを考慮して、可能な限り、詳細な信号データを維持する必要がある。つまり分解レベルを低く抑える必要がある。
　そこで、このレベルを調整するため、相関関数の閾値を定め、相関係数がこの閾値を超えれば、相関係数の計算を終了するようにする。そして、これにより定めた相関係数を基にＲＦＳを求めることで、適切な歪の値を計測することが可能となる。
　以上の処理の全体をまとめて、図２８に、ＡＣＤ法におけるデータ処理フローチャートとして示す。
　まず、最初に、分解レベルｎをｎ＝１とする（ステップＳ２１）。次に、２つのレイリースペクトルのデータを取得する（ステップＳ２２）。次に、各データの分解レベルｎでのデータをＤＷＴにより分解し、各データの詳細データ部を零にセットしメモリ・データベースに保存する（ステップＳ２３）。次に、ウェーブレット逆変換により、先に保存したデータのレイリースペクトルデータを再構築する（ステップＳ２４）。そして、再構築された２つのデータの相互相関係数を計算する（ステップＳ２５）。計算した相互相関係数の値を閾値（例えば０．４）と比較し（ステップＳ２６）、閾値より大きければデータ処理を終える。一方、計算した相互相関係数の値が閾値より小さければ、ｎの値を１だけ増やし（ステップＳ２７）、上記処理フローのステップＳ２３に戻り、上記ステップＳ２３からステップＳ２７までの処理を、計算した相互相関係数の値が閾値より大きくなるまで続ける。
　この処理フローにより、データの分解をすすめるとともに、より詳細な成分を可能な限り維持できる。そして、現在の計測場所の歪の値を決定するため、隣接する場所の歪の値を使用することなく、すべての場所の歪の値を個別に処理することができる。なお、この処理は別途説明するＦＦＴ処理を前提としたＦＣＤ法（fundamental components decoding method）と比較して計算時間を節約できる。
　そこで、次に、このＡＣＤ法の検証を行うため、先に図４Ａで示した、不均一な歪が検知されたデータを処理した結果について、図２９Ａ、図２９Ｂを用いて説明する。
　図２９Ａは、このＡＣＤ法を使って求めた、光ファイバ位置（横軸に示した距離を参照）に対応する光ファイバに生じた歪の分布を示したグラフである。なお、歪の分布を求めるに当たって、ここでは閾値として０．４を与えた（図２８参照）。
　また、図２９Ｂは、このＡＣＤ法を使って求めた、光ファイバ位置に対応する相関係数の分布である。なお、参考のために、先に図４Ａで示した従来の方法で求めた歪の分布のグラフを図２９Ｃに、また、これに対応する、相関係数の分布を図２９Ｄに示した。
　これらの図（図２９Ａ～図２９Ｄ）より、不均一歪があった場合でも、ＡＣＤ法が有効であることが判る。なお、図２９Ｂ、図２９Ｄにおいて、閾値として経験的に用いられる値０．３を点線で表示した。
　ＡＣＤ法が有効であることを別の観点から調べるため、次に、上記図２９Ａに示した歪分布について、ブリルアン散乱法を用いて得たデータと比較した。結果として得たデータを図３０Ａに示す。また、同様に、図３０Ａに示した２種類のデータを、移動平均法により平滑化したグラフを図３０Ｂに示す。
　これらのグラフから２種類のデータを比較すると、レイリー散乱光を基にして上記ＡＣＤ法で処理した歪のデータは、ブリルアン散乱法を用いて得たデータと比較して、同等に歪分布を測定できていると判断される。
　なお、このＡＣＤ法を不均一歪がある場所に適用して、処理が旨くできた点の数と全処理点の数（処理が旨くできた点と処理が旨くできなかった点を合計した点の数）との割合を成功率（単位：％）として、その成功率を調べた結果、９０．４６％の高い値を示した。
　以上により、ＡＣＤ法が有効であることが判る。
［ＦＣＤ法］
　上記解析法Ｎ２として、上記ＡＣＤ法と類似の解析法である基本成分複合法（ＦＣＤ法とも呼ぶ。ＦＣＤ：Fundamental Components Decoding）について以下説明する。すなわち、ＡＣＤ法に代えて、ＦＣＤ法を解析法Ｎ２として用いてもよい。以下、ＦＣＤ法の内容について、ＡＣＤ法と異なる点を中心に説明する。
　ＡＣＤ法では、基本的な解析の方式として、離散ウェーブレット変換（ＤＷＴ：Discrete Wavelet Transform）を用いたが、ＦＣＤ法では、ＦＦＴ（Fast Fourier Transform）を用いる。また、原理的には、ＡＣＤ法ではＲＩＰの近似データ部を抽出したが、ＦＣＤ法ではローパスフィルタ（以下、ＬＰＦとも呼ぶ）を利用して、ＲＩＰの低周波領域の要素を抽出する。一方、ＡＣＤ法では、使用するウェーブレット関数の選択が必要であるが、ＦＣＤ法ではその必要がない。逆に、ＡＣＤ法では初期値は不要であるが、ＦＣＤ法では初期値を与える必要がある。そこで、以下では、上記原理的な面から順にＦＣＤ法について詳しく説明する。
　被測定体に不均一な歪分布があると、そのレイリー周波数スペクトルは、高周波成分を持ち、そのことが原因で歪波形を歪ませることにつながる。一方、２つのレイリー周波数のデータに明らかに相関があるときには、図３１に一例を示したように、そのスペクトルは共通する基本成分（基本となる周波数成分。具体的には、図中、４個の矢印を付した箇所の周波数成分などを参照）を持っていると考えられる。
　そこで、これらの共通する基本成分、特に低周波領域の基本成分により、相関関係を求める手法である。このような共通する基本成分を得るために、ここではＬＰＦを使ってレイリー周波数信号の高周波成分をカットする（この際、最適なカットオフ周波数を選ぶ必要がある）。このとき、カットオフ周波数が低いほど、２つのレイリー周波数スペクトルの相関係数は高くなり、カットオフ周波数が高いほど、２つのレイリー周波数スペクトルの相関係数は低くなる。最適なカットオフ周波数を選ぶ（適切な値を見つける）ために、以下の２つの条件を考慮した。１つ目は、相関係数が所定の閾値を超えることであり、２つ目は、現在の観測場所でのシフト周波数はその前の観測場所のシフト周波数に近い必要がある。
　以上の概括的に述べた処理内容をデータ処理フローで図３２に示した。以下では、図３２を用いて、ＦＣＤ法におけるデータ処理フローについて説明する。
　図３２では、与えられた閾値Ａよりも高い相関係数を持つ一連のレイリー周波数シフトを計算する処理フローとなっている。そして、それらのレイリー周波数シフトを前の場所（位置番号：ｎ－１）のレイリー周波数シフトと比較して差の絶対値が最小となるものを選択して、現在の場所（位置番号：ｎ）のレイリー周波数シフトとする。
　具体的には、図３２に示すように、２つのレイリースペクトルデータ（データ１、データ２）をＬＰＦ（ローパスフィルタ）に通過させて（ステップＳ３１）、フィルタ通過後のデータ１とフィルタ通過後のデータ２をそれぞれ求め、これら２つのデータの相関を取る（ステップＳ３２）。その後、相関を取って得られた相関係数と所定の閾値Ａとの大小を比較する（ステップＳ３３）。そして、この相関係数が所定の閾値Ａより大きければ、上記ＬＰＦのカットオフ周波数を変更して（ステップＳ３４）、再度、２つのレイリースペクトルデータ（データ１、データ２）をＬＰＦ（ローパスフィルタ）に通過させ、再度、フィルタ通過後のデータ１とデータ２の相関を取り、閾値Ａとの比較を行う。ここで得られた相関係数が所定の閾値Ａ以上であれば、この時の相関係数からレイリー周波数シフトの値を求め（ステップＳ３５）、レイリー周波数シフトデータとしてメモリ等に保存する。そして、事前に得ていた上記前の場所のレイリー周波数シフトデータと、上記ステップＳ３１～ステップＳ３５の処理により得たレイリー周波数シフトデータとの差の絶対値を求め、その値が最小となるものを選択して（ステップＳ３６）、現在の場所（位置番号：ｎ）のレイリー周波数シフトとする。
　なお、初期値（つまりｎ＝１としたときの位置番号０での値）を定めるため、最初の場所を、既知のレイリー周波数シフトの場所から処理を開始する。このＦＣＤ法による処理を経て得られた相関係数を図３３Ａに示す。なお、比較のため、ＦＣＤ法による処理を経ずに得たデータを図３３Ｂに示した。ＦＣＤ法による処理を経たデータでは、明らかに相関係数が改善されていることが分かる。
　次に、このＦＣＤ法による処理により得られた周波数シフトデータを基に、先に不均一歪を示した距離を含む距離位置での歪分布データを求め、図３４Ａ、図３４Ｂに示した。図３４Ａには、このＦＣＤ法による処理により得られた歪分布データをそのまま、プロットして示した。また、図３４Ｂには、図３４Ａのデータを移動平均して平滑化処理したデータを示す。いずれのデータも、先に図４Ａに示した歪分布データに比較して、ＦＣＤ法の使用により、著しく改善されていることが分かる。これは先に説明したＡＣＤ法と同様である。
　最後に、ＦＣＤ法の使用による改善効果を示すため、ＦＣＤ法を適用した場合に求めたスペクトルの大きさと相関係数のスペクトル分布のデータと、何の処理も行わなかった（ＦＣＤ法を適用しなかった）場合に求めたスペクトルの大きさと相関係数のスペクトル分布のデータを、それぞれ、図３５Ａ、図３５Ｂ、図３５Ｃ、図３５Ｄにそれぞれ示した。
　以上、解析法Ｎ２は、実施の形態１で説明した解析法Ｎ１と同等の効果を持つことが分かる。また、解析法Ｎ２においては、ＡＣＤ法、およびＦＣＤ法のいずれの方法も有効であることを説明した。どちらの方法を選ぶかは、使用目的、各方法に用いる（例えばＦＣＤ法ではＦＦＴなどの）所有ツールの有無などの都合により、決定することができる。
　そして、この解析法Ｎ２を用いた場合でも、実施の形態１と同様の効果を得ることができる。
　なお、本願は、様々な例示的な実施の形態及び実施例が記載されているが、１つ、または複数の実施の形態に記載された様々な特徴、態様、及び機能は特定の実施の形態の適用に限られるのではなく、単独で、または様々な組み合わせで実施の形態に適用可能である。
従って、例示されていない無数の変形例が、本願明細書に開示される技術の範囲内において想定される。例えば、以上においては、計測される物理量として、主に歪を例に説明したが、歪を温度に置き換えても同様の議論ができるなど、少なくとも１つの構成要素を変形する場合、追加する場合または省略する場合、さらには、少なくとも１つの構成要素を抽出し、他の実施の形態の構成要素と組み合わせる場合が含まれるものとする。
　１　ＬＤ、２　光スイッチ、３　エルビウム添加光ファイバ増幅器（ＥＤＦＡ）、４　光カップラ、５　参照用光ファイバ、６　２ｃｈＡＤ変換器、７　第１の演算部、８　メモリ・データベース（メモリ・ＤＢ）、９　第２の演算部（演算部）、１１　光源部、１２　送受信光学部、１３　受信部、１４　レイリー強度パターンデジタル処理部（ＲＩＰデジタル処理部）、１００　レイリー強度パターン計測装置（ＲＩＰ計測装置）

Claims (11)

  1. 波長可変ＬＤと、前記波長可変ＬＤから発振されるレーザ光の周波数を変更可能な制御器と、を有する光源部、
    前記光源部から発振されたレーザ光を光ファイバに入射するとともに、前記光ファイバからのレイリー散乱光を前記レーザ光の入射経路とは別の経路で出射する光カップラ、
    前記光源部からのローカル光と、前記光カップラからの前記レイリー散乱光と、をコヒーレント受信する受信部、
    前記受信部から出力された信号をＡＤ変換するＡＤ変換器と、前記ＡＤ変換器で変換された信号を所定の方法で演算する第１の演算部と、前記第１の演算部で演算された信号が入力され保存されるメモリ・データベースと、前記メモリ・データベースに保存されたデータを基に所定の演算を行う第２の演算部と、を有し、前記レイリー散乱光の電界信号から得られる、予め定められた測定位置での異なる２つのレイリー強度パターンから、前記第２の演算部で前記測定位置について所定の補正をした後、相互相関係数を求め、求めた相互相関係数が予め定められた閾値以上となった場合の相互相関係数を前記メモリ・データベースに保存するレイリー強度パターンデジタル処理部、
    を備え、
    予め定められた閾値以上となった相互相関係数を基にレイリー周波数シフトを求めることにより、被測定体の歪分布、あるいは温度分布を定めることを特徴とするレイリー強度パターン計測装置。
  2. 前記光源部は局部発振器を備え、
    前記波長可変ＬＤの絶対周波数を制御するとともに、前記波長可変ＬＤから発振されるレーザの周波数をステップ状に走査し、各ステップでパルス化したチャープ信号を前記受信部で受信した後、前記第１の演算部で合成することを特徴とする請求項１に記載のレイリー強度パターン計測装置。
  3. 整合フィルタを有し、
    前記チャープ信号が前記整合フィルタにより抽出されたサブバンドにより、複数の窓関数を用いて設定されているとともに、前記抽出されたサブバンドの中心周波数を用いることにより、前記サブバンドを再度組み合わせた複数サブバンドを含む前記レイリー強度パターンが前記絶対周波数で保存されることを特徴とする請求項２に記載のレイリー強度パターン計測装置。
  4. 前記レイリー強度パターンデジタル処理部はＦＰＧＡあるいはＡＳＩＣを備え、
    前記第１の演算部で前記サブバンドの周波数成分をデジタル処理により取得するとともに、前記第２の演算部で周波数軸での信号の補間処理であるリサンプリング、および前記相互相関係数の演算を含む演算処理を行うことを特徴とする請求項３に記載のレイリー強度パターン計測装置。
  5. 前記第１の演算部で取得される前記サブバンドの周波数成分は、位相情報を含む非２乗処理のレイリー強度パターンの信号を基に形成されていることを特徴とする請求項４に記載のレイリー強度パターン計測装置。
  6. 前記第２の演算部は、
    前記２つのレイリー強度パターンを基に、各計測において、距離、周波数の双方を変数として求まる相互相関により、測定場所および前記レイリー周波数シフトを解析することを特徴とする請求項１から５のいずれか1項に記載のレイリー強度パターン計測装置。
  7. 請求項１に記載のレイリー強度パターン計測装置を用いて、レイリー強度パターン計測を行うレイリー強度パターン計測方法であって、
    前記求めた相互相関係数が予め定められた閾値より小さい場合には、前記第２の演算部で演算したレイリー強度パターンとは異なるレイリー強度パターンを用いて、与えられた閾値以上の相互相関係数が得られるまで繰り返し演算を行うことを特徴とするレイリー強度パターン計測方法。
  8. 請求項１に記載のレイリー強度パターン計測装置を用いて、レイリー強度パターン計測を行うレイリー強度パターン計測方法であって、
    前記求めた相互相関係数が予め定められた閾値より小さい場合に、
    レイリー散乱光の参照データの電界から得られる参照スペクトルと、観測データの電界から得られる観測スペクトルとを基に、各スペクトルの周波数をずらして各スペクトルの相互相関係数を求め、求めた相互相関係数が予め定められた閾値以上である相互相関係数を基に、被測定体のひずみ分布を検出することを特徴とするレイリー強度パターン計測方法。
  9. 請求項２に記載のレイリー強度パターン計測装置を用いて、レイリー強度パターン計測を行うレイリー強度パターン計測方法であって、
    前記求めた相互相関係数が予め定められた閾値より小さい場合に、
    レイリー散乱光の参照データの電界から得られる参照スペクトルを種々のチャープレートのチャープ信号で変形させておき、変形された前記参照スペクトルと観測スペクトルとの相互相関を取ることで、被測定体のレイリー散乱光の周波数シフトを計測することを特徴とするレイリー強度パターン計測方法。
  10. 請求項１に記載のレイリー強度パターン計測装置を用いて、レイリー強度パターン計測を行うレイリー強度パターン計測方法であって、
    前記求めた相互相関係数が予め定められた閾値より小さい場合に、
    データを所定の分解レベルで近似データ部と詳細データ部に分解する工程、
    前記詳細データ部をゼロとおいて、ウェーブレット逆変換により、データを再構築する工程、
    再構築されたデータを基に相互相関係数を求め、求めた相関係数が前記閾値以上となる場合の相関係数を算出する工程、
    を含む離散ウェーブレット変換手法を用いて傾きを含むレイリー散乱光の周波数シフト量を計測することを特徴とするレイリー強度パターン計測方法。
  11. 請求項１に記載のレイリー強度パターン計測装置を用いて、レイリー強度パターン計測を行うレイリー強度パターン計測方法であって、
    前記求めた相互相関係数が予め定められた閾値より小さい場合に、
    規定の測定場所に関して、２つのレイリースペクトルデータをローパスフィルタに通し、当該ローパスフィルタを通過後の２つのレイリースペクトルデータの相関を取った後、得られた相関係数と前記予め定められた閾値とは異なる別の閾値との大小を比較し、当該相関係数が前記別の閾値より小さければ、前記ローパスフィルタのカットオフ周波数を変更して、再度、２つのレイリースペクトルデータをローパスフィルタに通し、当該ローパスフィルタを通過後の２つのレイリースペクトルデータの相関を取り、前記別の閾値との比較を行い、得られた相関係数が前記別の閾値以上になった場合に、当該得られた相関係数からレイリー周波数シフトの値を求めて前記メモリ・データベースに保存するとともに、今回保存したレイリー周波数シフトのデータより前に保存していたレイリー周波数シフトのデータと、今回保存したレイリー周波数シフトのデータとの差の絶対値を複数個求め、求めた絶対値のうち、値が最小となるものを選択して、前記規定の測定場所のレイリー周波数シフトとすることを特徴とするレイリー強度パターン計測方法。
PCT/JP2021/037005 2021-10-06 2021-10-06 レイリー強度パターン計測装置およびレイリー強度パターン計測方法 WO2023058160A1 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
PCT/JP2021/037005 WO2023058160A1 (ja) 2021-10-06 2021-10-06 レイリー強度パターン計測装置およびレイリー強度パターン計測方法
JP2023552604A JPWO2023058160A1 (ja) 2021-10-06 2021-10-06

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2021/037005 WO2023058160A1 (ja) 2021-10-06 2021-10-06 レイリー強度パターン計測装置およびレイリー強度パターン計測方法

Publications (1)

Publication Number Publication Date
WO2023058160A1 true WO2023058160A1 (ja) 2023-04-13

Family

ID=85803314

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2021/037005 WO2023058160A1 (ja) 2021-10-06 2021-10-06 レイリー強度パターン計測装置およびレイリー強度パターン計測方法

Country Status (2)

Country Link
JP (1) JPWO2023058160A1 (ja)
WO (1) WO2023058160A1 (ja)

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7440087B2 (en) 2004-02-24 2008-10-21 Luna Innovations Incorporated Identifying optical fiber segments and determining characteristics of an optical device under test based on fiber segment scatter pattern data
JP2009042005A (ja) * 2007-08-07 2009-02-26 Nippon Telegr & Teleph Corp <Ntt> 光ファイバを用いた歪・温度の分布測定方法及び測定装置
WO2010061718A1 (ja) * 2008-11-27 2010-06-03 ニューブレクス株式会社 分布型光ファイバセンサ
JP2011232138A (ja) * 2010-04-27 2011-11-17 Neubrex Co Ltd 分布型光ファイバセンサ
WO2014083989A1 (ja) 2012-11-30 2014-06-05 ニューブレクス株式会社 ３次元位置計測装置
JP6552983B2 (ja) 2016-02-29 2019-07-31 ニューブレクス株式会社 ブリルアン散乱測定方法およびブリルアン散乱測定装置

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7440087B2 (en) 2004-02-24 2008-10-21 Luna Innovations Incorporated Identifying optical fiber segments and determining characteristics of an optical device under test based on fiber segment scatter pattern data
JP2009042005A (ja) * 2007-08-07 2009-02-26 Nippon Telegr & Teleph Corp <Ntt> 光ファイバを用いた歪・温度の分布測定方法及び測定装置
WO2010061718A1 (ja) * 2008-11-27 2010-06-03 ニューブレクス株式会社 分布型光ファイバセンサ
JP2011232138A (ja) * 2010-04-27 2011-11-17 Neubrex Co Ltd 分布型光ファイバセンサ
WO2014083989A1 (ja) 2012-11-30 2014-06-05 ニューブレクス株式会社 ３次元位置計測装置
JP6552983B2 (ja) 2016-02-29 2019-07-31 ニューブレクス株式会社 ブリルアン散乱測定方法およびブリルアン散乱測定装置

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
D. CHEN ET AL.: "Fiber-optic distributed acoustic sensor based on a chirped pulse and a non-matched filter", OPTICS EXPRESS, vol. 27, no. 20, 30 September 2019 (2019-09-30)
J. PASTOR-GRAELLS ET AL.: "Single-shot distributed temperature and strain tracking using direct detection phase-sensitive OTDR with chirped pulses", OPTICS EXPRESS, vol. 24, no. 12, 13 June 2016 (2016-06-13), XP055587402, DOI: 10.1364/OE.24.013121
S. LIEHR ET AL.: "Wavelength-Scanning Distributed Acoustic Sensing for Structual Monitoring and Seismic Application", PROCEEDINGS, vol. 15, 2019, pages 30
Y KOYAMADA: "Novel Technique to Improve Spatial Resolution in Brillouin Optical Time-Domain Reflectometry", IEEE PHOTONICS TECHNOLOGY LETTERS, vol. 19, no. 23, 1 December 2007 (2007-12-01), XP011197164, DOI: 10.1109/LPT.2007.908651
Y. WANG ET AL.: "Distributed Fiber-Optic Dynamic-Strain Sensor With Sub-Meter Spatial Resolution and Single-Shot Measurement", IEEE PHOTONIC JOURNAL, vol. 11, no. 6, December 2019 (2019-12-01), XP011755228, DOI: 10.1109/JPHOT.2019.2946293

Also Published As

Publication number Publication date
JPWO2023058160A1 (ja) 2023-04-13

Similar Documents

Publication Publication Date Title
CN104990620B (zh) 基于布拉格光纤光栅阵列的相敏光时域反射装置及方法
CN101470030B (zh) 用于将物理特性作为位置的函数确定的方法和系统
EP3477266B1 (en) Distributed acoustic sensing device using different coherent interrogating light patterns, and corresponding sensing method
US6788419B2 (en) Determination of properties of an optical device
US7515273B2 (en) Method for measuring the brillouin shift distribution along an optical fiber based on the optical demodulation of the signals, and relevant apparatus
US9958300B2 (en) Time division multiplexing (TDM) and wavelength division multiplexing (WDM) fast-sweep interrogator
US9551809B2 (en) Arrayed wave division multiplexing to improve spatial resolution of IOFDR fiber Bragg sensing system
JP6647420B2 (ja) ブリルアン散乱測定方法およびブリルアン散乱測定装置
AU2012321275B2 (en) Enhancing functionality of reflectometry based systems using parallel mixing operations
US9441947B2 (en) N-wavelength interrogation system and method for multiple wavelength interferometers
JP5043714B2 (ja) 光ファイバ特性測定装置及び方法
Lu et al. Numerical modeling of Fcy OTDR sensing using a refractive index perturbation approach
US7679729B2 (en) Light wave radar apparatus
CN108700401B (zh) 用干涉频谱法测量腔体
US9404831B2 (en) Arrayed wave division multiplex to extend range of IOFDR fiber bragg sensing system
US11549860B2 (en) Method and system for interrogating optical fibers
WO2023058160A1 (ja) レイリー強度パターン計測装置およびレイリー強度パターン計測方法
WO2022044174A1 (ja) 振動分布測定装置及びその方法
WO2020158033A1 (ja) 光パルス試験装置、及び光パルス試験方法
EP3150969B1 (en) Sensor for measuring the distribution of physical magnitudes in an optical fibre and associated measuring method
US20220291043A1 (en) Spectrum measurement method and spectrum measurement device
US20220349860A1 (en) Optical fiber distribution measurement system and signal processing method for optical fiber distribution measurement
CN116793239A (zh) 基于相位解调的光纤应变确定方法、装置和电子设备
Sahoo et al. Measurements and Signal Processing for Distributed Sensing with Centimeter-Scale Spatial Resolution by BOFDA

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: 21959902

Country of ref document: EP

Kind code of ref document: A1

WWE Wipo information: entry into national phase

Ref document number: 2023552604

Country of ref document: JP