WO2016103849A1 - 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム - Google Patents

超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム Download PDF

Info

Publication number
WO2016103849A1
WO2016103849A1 PCT/JP2015/078548 JP2015078548W WO2016103849A1 WO 2016103849 A1 WO2016103849 A1 WO 2016103849A1 JP 2015078548 W JP2015078548 W JP 2015078548W WO 2016103849 A1 WO2016103849 A1 WO 2016103849A1
Authority
WO
WIPO (PCT)
Prior art keywords
attenuation rate
ultrasonic
feature amount
unit
attenuation
Prior art date
Application number
PCT/JP2015/078548
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 EP15872418.7A priority Critical patent/EP3238632B1/en
Priority to CN201580004506.2A priority patent/CN105916448B/zh
Priority to JP2016505644A priority patent/JP5927367B1/ja
Publication of WO2016103849A1 publication Critical patent/WO2016103849A1/ja
Priority to US15/207,839 priority patent/US20160317126A1/en

Links

Images

Classifications

    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/52Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/5207Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving processing of raw data to produce diagnostic data, e.g. for generating an image
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/12Diagnosis using ultrasonic, sonic or infrasonic waves in body cavities or body tracts, e.g. by using catheters
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/44Constructional features of the ultrasonic, sonic or infrasonic diagnostic device
    • A61B8/4444Constructional features of the ultrasonic, sonic or infrasonic diagnostic device related to the probe
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/52Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/5269Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving detection or reduction of artifacts
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/52Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
    • G01S7/52017Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00 particularly adapted to short-range imaging
    • G01S7/52023Details of receivers
    • G01S7/52033Gain control of receivers
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/52Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
    • G01S7/52017Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00 particularly adapted to short-range imaging
    • G01S7/52023Details of receivers
    • G01S7/52036Details of receivers using analysis of echo signal for target characterisation
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/52Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
    • G01S7/52017Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00 particularly adapted to short-range imaging
    • G01S7/52053Display arrangements
    • G01S7/52057Cathode ray tube displays
    • G01S7/52071Multicolour displays; using colour coding; Optimising colour or information content in displays, e.g. parametric imaging
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/13Tomography
    • A61B8/14Echo-tomography
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/44Constructional features of the ultrasonic, sonic or infrasonic diagnostic device
    • A61B8/4483Constructional features of the ultrasonic, sonic or infrasonic diagnostic device characterised by features of the ultrasound transducer
    • AHUMAN NECESSITIES
    • A61MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61BDIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/52Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/5215Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving processing of medical diagnostic data
    • A61B8/5223Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving processing of medical diagnostic data for extracting a diagnostic or physiological parameter from medical diagnostic data

Definitions

  • the present invention relates to an ultrasonic observation apparatus that observes a tissue to be observed using ultrasonic waves, an operation method of the ultrasonic observation apparatus, and an operation program of the ultrasonic observation apparatus.
  • the present invention has been made in view of the above, and provides an ultrasonic observation apparatus, an operation method of the ultrasonic observation apparatus, and an ultrasonic observation apparatus capable of accurately distinguishing tissue characteristics of an observation target having a nonuniform attenuation rate, and An object is to provide an operation program of an ultrasonic observation apparatus.
  • an ultrasonic observation apparatus includes an ultrasonic transducer that transmits ultrasonic waves to an observation target and receives ultrasonic waves reflected by the observation target.
  • An ultrasonic observation apparatus that generates an ultrasonic image based on an ultrasonic signal acquired by an ultrasonic probe provided, and analyzing the frequency of the ultrasonic signal in the reception depth and reception direction of the ultrasonic signal.
  • a frequency analysis unit that calculates a plurality of frequency spectra in response, a feature amount calculation unit that calculates a feature amount of each of the plurality of frequency spectra, and each of the divided regions obtained by dividing the ultrasound image into a plurality of regions.
  • each frequency spectrum is used.
  • a pre-correction feature amount of each frequency spectrum for each attenuation rate candidate value is calculated by performing attenuation correction that eliminates the influence of the ultrasonic wave on the feature amount of the spectrum, and based on the calculation result,
  • An attenuation rate setting unit that sets an optimal attenuation rate for the observation target from among the attenuation rate candidate values, and among the optimal attenuation rates that the attenuation rate setting unit sets for each of the divided regions, the ultrasonic transducer
  • a cumulative attenuation rate per unit frequency at the sampling point is calculated using an optimum attenuation factor of the divided region existing between the surface and the sampling point, and attenuation correction of the feature amount is performed using the cumulative attenuation factor.
  • a feature amount correcting unit that calculates a corrected feature amount by performing the operation.
  • the feature amount correction unit may calculate the optimum attenuation rate for each divided region existing between the surface position of the ultrasonic transducer and the sampling point.
  • the cumulative attenuation rate at the sampling point is calculated by cumulatively adding the weights of the reciprocating distances in the depth direction in the divided areas.
  • the ultrasonic observation apparatus is the ultrasonic observation device according to the above aspect, wherein the two divided regions adjacent to each other in the depth direction have a length in the depth direction far from the ultrasonic transducer. It is more than the length of the said depth direction of the nearer side.
  • the attenuation rate setting unit calculates a statistical variation of the preliminary correction feature value for each attenuation rate candidate value, and the statistical variation is minimized.
  • a certain attenuation rate candidate value is set as the optimum attenuation rate.
  • the ultrasonic observation apparatus is characterized in that, in the above-described invention, the ultrasonic observation apparatus further includes a feature amount image data generation unit that generates feature amount image data for displaying information on the correction feature amount together with the ultrasound image. .
  • the feature amount calculation unit calculates the feature amount by performing a process of approximating each frequency spectrum by an n-order equation (n is a positive integer). It is characterized by that.
  • the feature amount calculation unit approximates a predetermined frequency band in the frequency spectrum by a linear expression, an intercept and an inclination of the linear expression, and an intermediate between the frequency bands.
  • One or a plurality of midband fits which is one of the values of the linear expression in frequency and including any one of the slope and the midband fit, are calculated as the feature amount, and the attenuation rate setting unit includes the slope and The optimum attenuation rate is set based on any one of the midband fits.
  • the attenuation rate setting unit sets the optimum attenuation rate based on the inclination when the inclination is the feature amount, and the midband fit is performed.
  • the optimum attenuation rate is set based on the midband fit.
  • An operation method of an ultrasonic observation apparatus is an ultrasonic signal acquired by an ultrasonic probe including an ultrasonic transducer that transmits ultrasonic waves to an observation target and receives ultrasonic waves reflected by the observation target.
  • a method of operating an ultrasonic observation apparatus that generates an ultrasonic image based on the frequency, wherein a frequency analysis unit analyzes a frequency of the ultrasonic signal to thereby analyze a plurality of frequencies according to a reception depth and a reception direction of the ultrasonic signal.
  • a frequency analysis step for calculating the frequency spectrum, a feature amount calculation unit for calculating a feature amount for each of the plurality of frequency spectra, and an attenuation rate setting unit for converting the ultrasonic image into a plurality of regions.
  • each of a plurality of unit lengths and attenuation rate candidate values per unit frequency giving different attenuation characteristics when the ultrasonic wave propagates through the observation target Is used to calculate a preliminary correction feature quantity of each frequency spectrum for each attenuation rate candidate value by performing attenuation correction that eliminates the influence of the ultrasonic wave on the feature quantity of each frequency spectrum.
  • An attenuation rate setting step for setting an optimum attenuation rate for the observation target from the plurality of attenuation rate candidate values based on the result, and an optimum attenuation rate set for each of the divided regions by the feature amount correction unit.
  • the operation program of the ultrasonic observation apparatus is an ultrasonic signal acquired by an ultrasonic probe including an ultrasonic transducer that transmits ultrasonic waves to an observation target and receives ultrasonic waves reflected by the observation target.
  • the frequency analysis unit calculates a plurality of frequency spectra corresponding to the reception depth and reception direction of the ultrasonic signal by analyzing the frequency of the ultrasonic signal.
  • Each of a plurality of unit lengths and attenuation rate candidate values per unit frequency giving different attenuation characteristics when the ultrasonic wave propagates through the observation target is used.
  • the preliminary correction feature amount of each frequency spectrum for each attenuation rate candidate value is calculated, and the calculation result is
  • An attenuation factor setting step for setting an optimum attenuation factor for the observation target from among the plurality of attenuation factor candidate values based on the optimum attenuation factor set for each of the divided regions.
  • a cumulative attenuation rate per unit frequency at the sampling point is calculated using an optimum attenuation factor of the divided region existing between the surface of the ultrasonic transducer and the sampling point, and the cumulative attenuation factor is used to calculate the cumulative attenuation factor.
  • a feature amount correction step of calculating a corrected feature amount by performing attenuation correction of the feature amount is performed.
  • the optimal attenuation rate for the observation target is set for each divided area based on the calculation result, and exists between the surface of the ultrasonic transducer and the sampling point.
  • FIG. 1 is a block diagram showing a functional configuration of an ultrasonic diagnostic system including an ultrasonic observation apparatus according to an embodiment of the present invention.
  • FIG. 2 is a diagram illustrating a relationship between the reception depth and the amplification factor in the amplification processing performed by the signal amplification unit of the ultrasonic observation apparatus according to the embodiment of the present invention.
  • FIG. 3 is a diagram illustrating a relationship between the reception depth and the amplification factor in the amplification correction process performed by the amplification correction unit of the ultrasonic observation apparatus according to the embodiment of the present invention.
  • FIG. 4 is a diagram schematically showing a data array in one sound ray of the ultrasonic signal.
  • FIG. 5 is a diagram illustrating an example of a frequency spectrum calculated by the frequency analysis unit of the ultrasonic observation apparatus according to the embodiment of the present invention.
  • FIG. 6 is a diagram schematically illustrating an example of setting divided areas in an ultrasound image.
  • FIG. 7 is a diagram showing a straight line having as a parameter the preliminary correction feature value corrected by the attenuation rate setting unit of the ultrasonic observation apparatus according to the embodiment of the present invention.
  • FIG. 8 is a diagram schematically illustrating an example of distribution of pre-correction feature quantities that have been subjected to attenuation correction based on two different attenuation rate candidate values for the same observation target.
  • FIG. 9 is a flowchart showing an outline of processing performed by the ultrasonic observation apparatus according to the embodiment of the present invention.
  • FIG. 10 is a flowchart showing an outline of processing executed by the frequency analysis unit of the ultrasonic observation apparatus according to the embodiment of the present invention.
  • FIG. 11 is a diagram showing an outline of processing performed by the attenuation rate setting unit of the ultrasonic observation apparatus according to the embodiment of the present invention.
  • FIG. 12 is a diagram illustrating another example of setting divided areas in an ultrasound image.
  • FIG. 13 is a diagram showing an outline of another optimum attenuation rate setting method performed by the attenuation rate setting unit of the ultrasonic observation apparatus according to the embodiment of the present invention.
  • FIG. 1 is a block diagram showing a functional configuration of an ultrasonic diagnostic system provided with an ultrasonic observation apparatus according to an embodiment of the present invention.
  • An ultrasonic diagnostic system 1 shown in FIG. 1 transmits an ultrasonic wave to a subject to be observed and receives an ultrasonic wave reflected by the subject, and an ultrasonic endoscope 2.
  • an ultrasonic observation device 3 that generates an ultrasonic image based on the ultrasonic signal acquired by the
  • a display device 4 that displays the ultrasonic image generated by the ultrasonic observation device 3.
  • the ultrasonic endoscope 2 converts an electrical pulse signal received from the ultrasonic observation device 3 into an ultrasonic pulse (acoustic pulse) and irradiates the subject at the tip thereof, and is reflected by the subject.
  • the ultrasonic transducer 21 converts the ultrasonic echo into an electrical echo signal expressed by a voltage change and outputs it.
  • the ultrasonic transducer 21 may be a convex transducer, a linear transducer, or a radial transducer.
  • the ultrasonic endoscope 2 may be one that mechanically scans the ultrasonic transducer 21, or a plurality of elements are provided in an array as the ultrasonic transducer 21, and the elements involved in transmission and reception are electronically arranged. Electronic scanning may be performed by switching or delaying transmission / reception of each element.
  • the ultrasonic transducer 21 is a linear transducer for convenience of explanation.
  • the ultrasonic endoscope 2 usually has an imaging optical system and an imaging device, and is inserted into the digestive tract (esophagus, stomach, duodenum, large intestine) or respiratory organ (trachea, bronchi) of the subject, and the digestive tract. It is possible to image the respiratory organ and surrounding organs (pancreas, gallbladder, bile duct, biliary tract, lymph node, mediastinal organ, blood vessel, etc.).
  • the ultrasonic endoscope 2 has a light guide that guides illumination light to be irradiated onto the subject during imaging.
  • the light guide has a distal end portion that reaches the distal end of the insertion portion of the ultrasonic endoscope 2 into the subject, and a proximal end portion that is connected to a light source device that generates illumination light.
  • the ultrasonic observation device 3 is electrically connected to the ultrasonic endoscope 2 and transmits a transmission signal (pulse signal) including a high voltage pulse to the ultrasonic transducer 21 based on a predetermined waveform and transmission timing.
  • a transmitter / receiver 31 that receives an echo signal that is an electrical reception signal from the ultrasonic transducer 21 and generates and outputs digital radio frequency (RF) data (hereinafter referred to as RF data);
  • RF data digital radio frequency
  • a signal processing unit 32 that generates digital B-mode reception data based on the RF data received from the unit 31; a calculation unit 33 that performs predetermined calculations on the RF data received from the transmission / reception unit 31;
  • An image processing unit 34 that generates data, and an input unit that is implemented using a user interface such as a keyboard, a mouse, and a touch panel, and receives input of various types of information It includes a 5, a control unit 36 for controlling the entire ultrasound diagnostic system 1, a storage unit 37 that stores various information
  • the transmission / reception unit 31 includes a signal amplification unit 311 that amplifies the echo signal.
  • the signal amplification unit 311 performs STC (Sensitivity Time Control) correction in which an echo signal having a larger reception depth is amplified with a higher amplification factor.
  • FIG. 2 is a diagram illustrating the relationship between the reception depth and the amplification factor in the amplification process performed by the signal amplification unit 311.
  • the reception depth z shown in FIG. 2 is an amount calculated based on the elapsed time from the reception start point of the ultrasonic wave. As shown in FIG.
  • the amplification factor ⁇ (dB) increases linearly from ⁇ 0 to ⁇ th (> ⁇ 0 ) as the reception depth z increases.
  • the amplification factor ⁇ (dB) takes a constant value ⁇ th when the reception depth z is equal to or greater than the threshold value z th .
  • the value of the threshold value z th is such a value that the ultrasonic signal received from the observation target is almost attenuated and the noise becomes dominant. More generally, when the reception depth z is smaller than the threshold value z th , the amplification factor ⁇ may increase monotonously as the reception depth z increases.
  • the relationship shown in FIG. 2 is stored in the storage unit 37 in advance.
  • the transmission / reception unit 31 performs processing such as filtering on the echo signal amplified by the signal amplification unit 311 and then performs A / D conversion to generate time domain RF data, and the signal processing unit 32 and the calculation unit To 33.
  • the transmission / reception unit 31 is a multi-channel circuit for beam synthesis corresponding to the plurality of elements.
  • the frequency band of the pulse signal transmitted by the transmission / reception unit 31 may be a wide band that substantially covers the linear response frequency band of the electroacoustic conversion of the pulse signal to the ultrasonic pulse in the ultrasonic transducer 21.
  • the various processing frequency bands of the echo signal in the signal amplifying unit 311 may be a wide band that substantially covers the linear response frequency band of the acoustoelectric conversion of the ultrasonic transducer 21 into the echo signal of the ultrasonic echo. Accordingly, it is possible to perform accurate approximation when performing frequency spectrum approximation processing, which will be described later.
  • the transmission / reception unit 31 transmits various control signals output from the control unit 36 to the ultrasonic endoscope 2 and receives various types of information including an identification ID from the ultrasonic endoscope 2 and receives the control unit 36. It also has a function to transmit to.
  • the signal processing unit 32 performs known processing such as band-pass filter, envelope detection, and logarithmic conversion on the RF data to generate digital B-mode reception data. In logarithmic conversion, a common logarithm of an amount obtained by dividing RF data by a reference voltage is taken and expressed as a decibel value.
  • the signal processing unit 32 outputs the generated B-mode reception data to the image processing unit 34.
  • the signal processing unit 32 is realized using a CPU (Central Processing Unit), various arithmetic circuits, and the like.
  • the calculation unit 33 performs amplification correction on the RF data output from the transmission / reception unit 31 so as to make the amplification factor constant regardless of the reception depth, and fast Fourier transform (RF Fourier transform) on the RF data on which the amplification correction has been performed.
  • RF Fourier transform fast Fourier transform
  • a frequency analysis unit 332 that calculates a plurality of frequency spectra according to the reception depth and reception direction of an ultrasonic signal by performing frequency analysis by applying FFT (Fast Fourier Transform), and a feature that calculates a feature amount of each frequency spectrum Attenuation per unit length and unit frequency giving different attenuation characteristics when the ultrasonic wave propagates through the observation object in each divided region when the amount calculation unit 333 divides the ultrasonic image into a plurality of regions By using each of the candidate rate values, attenuation is performed by applying attenuation correction that eliminates the influence of ultrasonic waves on the feature quantity of each frequency spectrum.
  • FFT Fast Fourier Transform
  • An attenuation rate setting unit 334 that calculates a preliminary correction feature amount of each frequency spectrum for each rate candidate value, and sets an optimum attenuation rate for an observation target from a plurality of attenuation rate candidate values based on the calculation result; Of the optimum attenuation factors set by the attenuation factor setting unit 334 for each divided region, the optimum attenuation factor of the divided region existing between the surface of the ultrasonic transducer and the sampling point is used to obtain the per unit frequency at the sampling point. And a feature amount correction unit 335 that calculates a corrected feature amount by calculating an accumulated attenuation rate and performing attenuation correction of the feature amount using the accumulated attenuation rate.
  • the calculation unit 33 is realized using a CPU, various calculation circuits, and the like.
  • FIG. 3 is a diagram illustrating a relationship between the reception depth and the amplification factor in the amplification correction process performed by the amplification correction unit 331.
  • the amplification rate ⁇ (dB) in the amplification process performed by the amplification correction unit 331 takes the maximum value ⁇ th ⁇ 0 when the reception depth z is zero, and the reception depth z is zero to the threshold value z th. Decreases linearly until reaching 0 and is zero when the reception depth z is greater than or equal to the threshold z th .
  • the relationship shown in FIG. 3 is stored in the storage unit 37 in advance.
  • the amplification correction unit 331 amplifies and corrects the digital RF signal based on the relationship shown in FIG.
  • the relationship between the reception depth z and the amplification factor ⁇ in the amplification correction process performed by the amplification correction unit 331 varies depending on the relationship between the reception depth and the amplification factor in the amplification process performed by the signal amplification unit 311.
  • STC correction is a correction process that eliminates the influence of attenuation from the amplitude of the analog signal waveform by amplifying the amplitude of the analog signal waveform uniformly over the entire frequency band and with a gain that monotonously increases with respect to the depth. is there. For this reason, when generating a B-mode image to be displayed by converting the amplitude of the echo signal into luminance, and when scanning a uniform tissue, the luminance value is constant regardless of the depth by performing STC correction. become. That is, an effect of eliminating the influence of attenuation from the luminance value of the B-mode image can be obtained.
  • a reception signal subjected to STC correction is generated when generating a B-mode image, while transmission for generating a B-mode image is generated when generating an image based on a frequency spectrum. It is conceivable to perform a new transmission different from the above and output a reception signal not subjected to STC correction. However, in this case, there is a problem that the frame rate of the image data generated based on the received signal is lowered.
  • the amplification correction unit 331 Correct the gain.
  • the frequency analysis unit 332 samples the RF data (line data) of each sound ray amplified and corrected by the amplification correction unit 331 at a predetermined time interval to generate sample data.
  • the frequency analysis unit 332 calculates a frequency spectrum at a plurality of locations (data positions) on the RF data by performing FFT processing on the sample data group.
  • FIG. 4 is a diagram schematically showing a data array in one sound ray of the ultrasonic signal.
  • a white or black rectangle means data at one sample point.
  • the data located on the right side is sample data from a deeper location when measured from the ultrasonic transducer 21 along the sound ray SR k (see the arrow in FIG. 4).
  • the sound ray SR k is discretized at a time interval corresponding to a sampling frequency (for example, 50 MHz) in A / D conversion performed by the transmission / reception unit 31.
  • FIG. 1 a sampling frequency
  • a sample data group needs to have a power number of 2 data.
  • a process for generating a normal sample data group is performed by inserting zero data in an insufficient amount. This point will be described in detail when the processing of the frequency analysis unit 332 is described (see FIG. 10).
  • FIG. 5 is a diagram illustrating an example of a frequency spectrum calculated by the frequency analysis unit 332.
  • the “frequency spectrum” means “frequency distribution of intensity at a certain reception depth z” obtained by performing FFT processing on a sample data group.
  • intensity refers to parameters such as the voltage of the echo signal, the power of the echo signal, the sound pressure of the ultrasonic echo, the acoustic energy of the ultrasonic echo, the amplitude and time integral value of these parameters, and combinations thereof. Points to either.
  • the horizontal axis represents the frequency f.
  • the reception depth z is constant. It will be described later linear L 10 shown in FIG. In the present embodiment, the curve and the straight line are composed of a set of discrete points.
  • the lower limit frequency f L and the upper limit frequency f H of the frequency band used for the subsequent calculation are the frequency band of the ultrasonic transducer 21 and the frequency band of the pulse signal transmitted by the transmitting / receiving unit 31.
  • f L 3 MHz
  • f H 10 MHz.
  • the frequency band determined by the lower limit frequency f L and the upper limit frequency f H is referred to as “frequency band U”.
  • the frequency spectrum shows a tendency that varies depending on the properties of the living tissue scanned with ultrasonic waves. This is because the frequency spectrum has a correlation with the size, number density, acoustic impedance, and the like of the scatterer that scatters ultrasonic waves.
  • the “characteristics of the biological tissue” referred to here includes, for example, malignant tumor (cancer), benign tumor, endocrine tumor, mucinous tumor, normal tissue, cyst, vascular vessel and the like.
  • the feature quantity calculation unit 333 calculates a feature quantity that characterizes the approximated primary expression by performing regression analysis on the frequency spectrum in the predetermined frequency band and approximating it with the linear expression. For example, in the case of the frequency spectrum C 1 illustrated in FIG. 5, the feature amount calculation unit 333 obtains the approximate straight line L 10 by performing regression analysis in the frequency band U.
  • the feature amount calculation unit 333 uses the slope a 0 , the intercept b 0 , and the frequency band as the feature amount corresponding to the straight line L 10.
  • the feature amount calculation unit 333 may approximate the frequency spectrum with a second-order or higher polynomial by regression analysis.
  • the slope a 0 has a correlation with the size of the ultrasonic scatterer, and it is generally considered that the larger the scatterer, the smaller the slope.
  • the intercept b 0 has a correlation with the size of the scatterer, the difference in acoustic impedance, the number density (concentration) of the scatterer, and the like. Specifically, the intercept b 0 has a larger value as the scatterer is larger, has a larger value as the difference in acoustic impedance is larger, and has a larger value as the number density of the scatterers is larger.
  • the mid-band fit c 0 is an indirect parameter derived from the slope a 0 and the intercept b 0 and gives the intensity of the spectrum at the center in the effective frequency band. Therefore, the midband fit c 0 is considered to have a certain degree of correlation with the brightness of the B-mode image in addition to the size of the scatterer, the difference in acoustic impedance, and the number density of the scatterers.
  • the attenuation rate setting unit 334 corrects the attenuation of the feature amount using the attenuation rate that gives the attenuation amount of the ultrasonic wave per unit length and unit frequency in each divided region when the ultrasonic image is divided into a plurality of regions.
  • FIG. 6 is a diagram schematically illustrating an example of setting divided areas in an ultrasound image.
  • the broken line of FIG. 6 represents the sound ray.
  • all the divided areas P ij have the same area.
  • the length in the depth direction from the surface location 102 of the ultrasonic transducer 21 in the divided region P ij H (hereinafter, referred to as the height H of the divided area P ij) is constant.
  • Information regarding the divided area P ij is stored in a divided area information storage unit 371 included in the storage unit 37.
  • FIG. 6 is an ultrasonic image when the ultrasonic transducer 21 is a linear transducer, and thus has a rectangular shape.
  • the outer shape of the ultrasonic image is illustrated. Needless to say, they are different.
  • the ultrasonic attenuation A (f, z) is attenuation that occurs while the ultrasonic waves reciprocate between the reception depth 0 and the reception depth z, and the intensity change before and after the reciprocation (difference in decibel expression). ).
  • the attenuation amount A (f, z) is empirically known to be proportional to the frequency in a uniform tissue, and is expressed by the following equation (1).
  • a (f, z) 2 ⁇ zf (1)
  • the proportional constant ⁇ is an amount called an attenuation rate, and gives an attenuation amount of ultrasonic waves per unit length and unit frequency.
  • Z is the ultrasonic reception depth, and f is the frequency.
  • a specific value of the attenuation rate ⁇ is determined according to the part of the living body.
  • the unit of the attenuation rate ⁇ is, for example, dB / cm / MHz.
  • the attenuation rate setting unit 334 sets an optimal attenuation rate from a plurality of attenuation rate candidate values. At this time, the attenuation rate setting unit 334 uses the attenuation rate candidate value ⁇ to perform the following with respect to the feature amounts (slope a 0 , intercept b 0 , midband fit c 0 ) calculated by the feature amount calculation unit 333. Preliminary correction feature values a, b, and c are calculated by performing attenuation correction according to equations (2) to (4).
  • the attenuation rate setting unit 334 performs correction with a larger correction amount as the ultrasonic reception depth z increases.
  • the correction related to the intercept is an identity transformation. This is because the intercept is a frequency component corresponding to a frequency of 0 (Hz) and is not affected by attenuation.
  • FIG. 7 is a diagram illustrating a straight line having the preliminary correction feature amounts a, b, and c corrected by the attenuation rate setting unit 334 as parameters.
  • the attenuation rate setting unit 334 sets the attenuation rate candidate value having the smallest statistical variation of the preliminary correction feature value calculated for each attenuation rate candidate value as the optimum attenuation rate in each divided region. In this embodiment, dispersion is applied as an amount indicating statistical variation. In this case, the attenuation rate setting unit 334 sets the attenuation rate candidate value that minimizes the variance as the optimal attenuation rate. Two of the three preliminary correction feature values a, b, and c described above are independent. In addition, the preliminary correction feature amount b does not depend on the attenuation rate. Therefore, when setting an optimal attenuation rate for the preliminary correction feature amounts a and c, the attenuation rate setting unit 334 may calculate the variance of one of the preliminary correction feature amounts a and c.
  • the attenuation rate setting unit 334 sets an optimal attenuation rate using the preliminary correction feature amount a
  • the dispersion of the preliminary correction feature amount a is applied and the optimal attenuation rate is set using the preliminary correction feature amount c.
  • a (f, z) 2 ⁇ zf + 2 ⁇ 1 z (6) ⁇ 1 in the second term on the right-hand side of equation (6) is a coefficient representing the magnitude of the change in signal intensity in proportion to the ultrasonic reception depth z.
  • Equation (6) It is a coefficient representing a change in signal intensity caused by a change in the number of channels at the time of synthesis. Since the second term on the right side of Equation (6) exists, when setting the optimal attenuation rate using the preliminary correction feature quantity c, the attenuation is more accurately corrected by applying the variance of the preliminary correction feature quantity c. (See equation (4)). On the other hand, when the optimum attenuation rate is set using the preliminary correction feature quantity a that is a coefficient proportional to the frequency f, the second term on the right side of the equation (6) is more suitable when the variance of the preliminary correction feature quantity a is applied. It is possible to correct the attenuation accurately by eliminating the influence of.
  • the reason why the optimum attenuation rate can be set based on statistical variation will be described.
  • the feature amount is converged to a value unique to the observation target regardless of the distance between the observation target and the ultrasonic transducer 21, and the statistical variation is considered to be small.
  • the attenuation rate candidate value that does not match the observation target is set as the optimal attenuation rate, the attenuation correction is excessive or insufficient, and thus the feature amount is shifted depending on the distance from the ultrasonic transducer 21. It is considered that the statistical variation of the feature amount is increased. Therefore, it can be said that the attenuation rate candidate value having the smallest statistical variation is the optimum attenuation rate for the observation target.
  • FIG. 8 is a diagram schematically illustrating an example of distribution of pre-correction feature quantities that have been subjected to attenuation correction based on two different attenuation rate candidate values for the same observation target.
  • the horizontal axis represents the preliminary correction feature amount
  • the vertical axis represents the frequency.
  • the two distribution curves N 1 and N 2 shown in FIG. 8 have the same total frequency.
  • the distribution curve N 1 has a small statistical variation (small variance) in the feature amount compared to the distribution curve N 2, and has a steep mountain shape.
  • the attenuation rate setting unit 334 sets an optimum attenuation rate from the two attenuation rate candidate values corresponding to the two distribution curves N 1 and N 2 . Therefore, when the attenuation rate setting unit 334 sets the attenuation rate candidate value corresponding to the distribution curve N 1. Set as the optimum attenuation rate.
  • the feature amount correction unit 335 calculates a cumulative attenuation rate per unit frequency (hereinafter also simply referred to as a cumulative attenuation rate) at a sampling point using an optimum attenuation rate for each divided region, and uses the accumulated attenuation rate to perform a feature. Attenuate the amount of attenuation.
  • the cumulative attenuation rate at an arbitrary sampling point is calculated by using the distance from the surface of the ultrasonic transducer 21 and the optimum attenuation rate in the divided area existing therebetween.
  • Equation (7) The first term on the right-hand side of Equation (7) is the sum of the ultrasonic round-trip distance 2H for each divided region multiplied by the optimum attenuation rate ⁇ (P Ij ).
  • the feature amount correction unit 335 calculates the cumulative attenuation rate ⁇ IJ (h) by accumulating the attenuation rate from the surface of the ultrasonic transducer 21.
  • the unit of cumulative attenuation rate is dB / MHz.
  • the image processing unit 34 converts the amplitude of the echo signal into luminance and displays a B-mode image data 341 that generates B-mode image data that is an ultrasonic image to be displayed, and a feature amount calculated by the feature amount calculation unit 333.
  • a feature amount image data generation unit 342 that generates feature amount image data for displaying information.
  • the B-mode image data generation unit 341 performs signal processing using known techniques such as gain processing and contrast processing on the B-mode reception data received from the signal processing unit 32, and displays an image on the display device 4.
  • B-mode image data is generated by thinning out data according to the data step width determined according to the range.
  • the B-mode image is a grayscale image in which values of R (red), G (green), and B (blue), which are variables when the RGB color system is adopted as a color space, are matched.
  • the B-mode image data generation unit 341 performs coordinate conversion on the B-mode reception data so that the scanning range can be spatially represented correctly, and then performs interpolation processing between the B-mode reception data to perform B-mode image data. A gap between received data is filled to generate B-mode image data.
  • the B-mode image data generation unit 341 outputs the generated B-mode image data to the feature amount image data generation unit 342.
  • the feature amount image data generation unit 342 generates feature amount image data by associating a hue as visual information with any one of the above-described inclination, intercept, and midband fit, for example.
  • the feature amount image data generation unit 342 generates the feature amount image data by associating the hue with one of the two feature amounts selected from the inclination, the intercept, and the midband fit, and by associating the light and dark with the other. May be.
  • a predetermined color system such as saturation, luminance value, R (red), G (green), and B (blue) is configured. List color space variables.
  • the control unit 36 is realized using a CPU (Central Processing Unit) having various calculation and control functions, various arithmetic circuits, and the like.
  • the control unit 36 controls the ultrasonic observation apparatus 3 in an integrated manner by reading information stored and stored in the storage unit 37 from the storage unit 37 and executing various arithmetic processes related to the operation method of the ultrasonic observation apparatus 3. To do.
  • the control unit 36 may be configured using a CPU or the like common to the signal processing unit 32 and the calculation unit 33.
  • the storage unit 37 includes a divided region information storage unit 371 that stores information related to divided regions, a spectrum information storage unit 372 that stores information on the frequency spectrum calculated by the frequency analysis unit 332 together with a reception depth and a reception direction, and a feature amount calculation.
  • the feature amount information storage unit 373 that stores the feature amount calculated by the unit 333 and the corrected feature amount corrected by the feature amount correction unit 335, and the optimum attenuation rate and feature that the attenuation rate setting unit 334 sets for each divided region
  • An attenuation rate information storage unit 374 that stores information on the accumulated attenuation rate for each sampling point calculated by the amount correction unit 335.
  • the storage unit 37 has information necessary for amplification processing (relationship between the amplification factor and the reception depth shown in FIG. 2) and information necessary for amplification correction processing (the amplification factor and the reception depth shown in FIG. 3). ), Information necessary for the attenuation correction process (see equation (1)), information on window functions (Hamming, Hanning, Blackman, etc.) necessary for the frequency analysis process, and the like are stored.
  • the storage unit 37 stores various programs including an operation program for executing the operation method of the ultrasonic observation apparatus 3.
  • the operation program can be recorded on a computer-readable recording medium such as a hard disk, a flash memory, a CD-ROM, a DVD-ROM, or a flexible disk and widely distributed.
  • the various programs described above can also be obtained by downloading via a communication network.
  • the communication network here is realized by, for example, an existing public line network, LAN (Local Area Network), WAN (Wide Area Network), etc., and may be wired or wireless.
  • the storage unit 37 having the above configuration is realized using a ROM (Read Only Memory) in which various programs are installed in advance, and a RAM (Random Access Memory) that stores calculation parameters and data of each process. .
  • ROM Read Only Memory
  • RAM Random Access Memory
  • FIG. 9 is a flowchart showing an outline of processing performed by the ultrasonic observation apparatus 3 having the above configuration. Specifically, it is a flowchart showing an outline of processing after the ultrasound observation apparatus 3 receives an echo signal from the ultrasound endoscope 2.
  • the ultrasonic observation device 3 receives an echo signal as a measurement result of an observation target by the ultrasonic transducer 21 from the ultrasonic endoscope 2 (step S1).
  • the signal amplifying unit 311 that has received the echo signal from the ultrasonic transducer 21 amplifies the echo signal (step S2).
  • the signal amplifying unit 311 performs amplification (STC correction) of the echo signal based on the relationship between the amplification factor and the reception depth shown in FIG. 2, for example.
  • the B-mode image data generation unit 341 generates B-mode image data using the echo signal amplified by the signal amplification unit 311 and outputs the B-mode image data to the display device 4 (step S3).
  • the display device 4 that has received the B-mode image data displays a B-mode image corresponding to the B-mode image data.
  • the amplification correction unit 331 performs amplification correction on the RF data output from the transmission / reception unit 31 so that the amplification factor is constant regardless of the reception depth (step S4).
  • the amplification correction unit 331 performs amplification correction based on, for example, the relationship between the amplification factor and the reception depth shown in FIG.
  • the frequency analysis unit 332 calculates frequency spectra for all the sample data groups by performing frequency analysis by FFT on the RF data of each sound ray after amplification correction, and stores the frequency spectrum in the spectrum information storage unit 372. (Step S5).
  • FIG. 10 is a flowchart showing an outline of the processing performed by the frequency analysis unit 332 in step S5.
  • the frequency analysis processing will be described in detail with reference to the flowchart shown in FIG.
  • the frequency analysis unit 332 sets a counter k for identifying a sound ray to be analyzed as k 0 (step S21).
  • the frequency analysis unit 332 sets an initial value Z (k) 0 of a data position (corresponding to a reception depth) Z (k) that represents a series of data groups (sample data group) generated for the FFT calculation.
  • FIG. 4 shows a case where the eighth data position of the sound ray SR k is set as the initial value Z (k) 0 as described above.
  • the frequency analysis unit 332 acquires a sample data group (step S23), and causes the window function stored in the storage unit 37 to act on the acquired sample data group (step S24).
  • the window function By applying the window function to the sample data group in this way, it is possible to avoid the sample data group from becoming discontinuous at the boundary and to prevent the occurrence of artifacts.
  • the frequency analysis unit 332 determines whether or not the sample data group at the data position Z (k) is a normal data group (step S25).
  • the sample data group needs to have the number of powers of two.
  • the number of data in the normal sample data group is 2 n (n is a positive integer).
  • step S25 If the result of determination in step S25 is that the sample data group at data position Z (k) is normal (step S25: Yes), the frequency analysis unit 332 proceeds to step S27 described later.
  • step S25 If the result of determination in step S25 is that the sample data group at the data position Z (k) is not normal (step S25: No), the frequency analysis unit 332 inserts zero data as much as the deficient amount into a normal sample data group. Generate (step S26). A window function is applied to the sample data group determined to be not normal in step S25 (for example, the sample data group F K in FIG. 4) before adding zero data. For this reason, even if zero data is inserted into the sample data group, discontinuity of data does not occur. After step S26, the frequency analysis unit 332 proceeds to step S27 described later.
  • step S27 the frequency analysis unit 332 obtains a frequency spectrum that is a frequency distribution of the amplitude by performing an FFT operation using the sample data group (step S27).
  • the frequency analysis unit 332 changes the data position Z (k) by the step width D (step S28). It is assumed that the step width D is stored in advance in the storage unit 37.
  • the step width D is desirably matched with the data step width used when the B-mode image data generation unit 341 generates B-mode image data.
  • the step width D A value larger than the data step width may be set as the width D.
  • the frequency analysis unit 332 determines whether or not the data position Z (k) is larger than the maximum value Z (k) max in the sound ray SR k (step S29).
  • the frequency analysis unit 332 increases the counter k by 1 (step S30). This means that the processing is shifted to the next sound ray.
  • the frequency analysis unit 332 returns to step S23.
  • the frequency analysis unit 332 determines whether or not the counter k is greater than the maximum value k max (step S31). When the counter k is larger than k max (step S31: Yes), the frequency analysis unit 332 ends the series of frequency analysis processing. On the other hand, when the counter k is equal to or less than k max (step S31: No), the frequency analysis unit 332 returns to step S22.
  • the maximum value k max is a value arbitrarily input by a user such as an operator through the input unit 35 or a value preset in the storage unit 37.
  • the frequency analysis unit 332 performs the FFT operation a plurality of times for each of (k max ⁇ k 0 +1) sound rays in the analysis target region.
  • the frequency spectrum obtained as a result of the FFT operation is stored in the spectrum information storage unit 372 together with the reception depth and the reception direction.
  • the frequency analysis unit 332 performs frequency analysis processing on all the areas where the ultrasonic signal is received.
  • the input unit 35 is divided by a specific depth width and sound ray width. It is also possible to accept the setting input of the region, and to perform the frequency analysis process only in the set partial region.
  • the attenuation rate setting unit 334 sets a value of the attenuation rate candidate value ⁇ to be applied when performing attenuation correction, which will be described later, to a predetermined initial value ⁇ 0 (step S7).
  • the initial value ⁇ 0 may be stored in advance in the attenuation rate information storage unit 374.
  • the attenuation rate setting unit 334 calculates a preliminary correction feature amount by performing attenuation correction on the feature amount approximated to each frequency spectrum by the feature amount calculation unit 333 using the attenuation rate candidate value as ⁇ . Then, it is stored in the feature amount information storage unit 373 together with the attenuation rate candidate value ⁇ (step S8).
  • a straight line L 1 illustrated in FIG. 7 is an example of a straight line obtained by the attenuation rate setting unit 334 performing attenuation correction processing.
  • f sp is the data sampling frequency
  • v s is the sound speed
  • D is the data step width
  • n is the number of data steps from the first data of the sound ray up to the data position of the amplitude data group to be processed.
  • the sampling frequency f sp data and 50 MHz, the sound velocity v s and 1530 m / sec, when a 15 step width D employs a data sequence shown in FIG. 4, a z 0.2295n (mm).
  • the attenuation rate setting unit 334 calculates the variance of representative preliminary correction feature amounts among a plurality of preliminary correction feature amounts obtained by the attenuation rate setting unit 334 performing attenuation correction on each frequency spectrum, and sets the attenuation rate.
  • the feature value is stored in the feature amount information storage unit 373 in association with the candidate value ⁇ (step S9). For example, when the preliminary correction feature amount is the slope a and the midband fit c, as described above, the attenuation rate setting unit 334 calculates the variance of one of the preliminary correction feature amounts a and c.
  • this step S9 when the feature amount image is generated using the preliminary correction feature amount a, the distribution of the preliminary correction feature amount a in each divided region is applied, and the feature amount image is generated using the preliminary correction feature amount c. In this case, it is preferable to apply the variance of the preliminary correction feature quantity c in each divided area.
  • the attenuation rate setting unit 334 increases the value of the attenuation rate candidate value ⁇ by ⁇ (step S10), and compares the increased attenuation rate candidate value ⁇ with a predetermined maximum value ⁇ max (step S10). S11). If the attenuation rate candidate value ⁇ is greater than the maximum value ⁇ max as a result of the comparison in step S11 (step S11: Yes), the ultrasound observation apparatus 3 proceeds to step S12. On the other hand, when the attenuation rate candidate value ⁇ is equal to or less than the maximum value ⁇ max as a result of the comparison in step S11 (step S11: No), the ultrasound observation apparatus 3 returns to step S8.
  • the attenuation rate setting unit 334 refers to the variance of the preliminary correction feature amount for each attenuation rate candidate value stored in the feature amount information storage unit 373 for each divided region, and the attenuation rate candidate having the minimum variance The value is set as the optimum attenuation rate in the divided area (step S12).
  • FIG. 11 is a diagram illustrating an outline of processing performed by the attenuation rate setting unit 334.
  • 0.2 (dB / cm / MHz)
  • the feature amount correction unit 335 calculates the cumulative attenuation rate at the sampling point using the optimal attenuation rate set by each attenuation region by the attenuation rate setting unit 334 (step S13).
  • the cumulative attenuation rate ⁇ IJ (h) at the sampling point S IJ (h) in the divided region P IJ shown in FIG. 6 is expressed by Expression (7).
  • the feature amount correcting unit 335 calculates a corrected feature amount by performing attenuation correction of the feature amount using the cumulative attenuation rate (step S14).
  • the correction feature amounts a IJ (h), b IJ (h), c IJ of the slope a 0 , the intercept b 0 , and the midband fit c 0 at the sampling point S IJ (h) of the divided region P IJ shown in FIG. (H) is calculated using equations (9) to (11), respectively.
  • the feature amount image data generation unit 342 superimposes visual information (for example, hue) associated with the corrected feature amount calculated in step S14 on each pixel in the B mode image data generated by the B mode image data generation unit 341. Thus, feature amount image data is generated (step S15).
  • the feature amount image data generation unit 342 transmits the generated feature amount image data to the display device 4.
  • the display device 4 that has received the feature amount image data displays a feature amount image corresponding to the received feature amount image data.
  • step S15 the ultrasonic observation apparatus 3 ends a series of processes. Note that the ultrasound observation apparatus 3 periodically repeats the processes of steps S1 to S15.
  • a plurality of unit lengths and attenuation rate candidate values per unit frequency giving different attenuation characteristics are obtained. Based on the result of this calculation, a pre-correction feature value for the frequency spectrum is calculated, and from among the multiple attenuation rate candidate values, the optimal attenuation rate for the observation target is set for each divided region, and the surface of the ultrasonic transducer is set.
  • Correction is performed by calculating the accumulated attenuation rate per unit frequency at the sampling point using the optimum attenuation rate of the divided region existing between the sampling point and the sampling point, and correcting the attenuation of the feature value using the accumulated attenuation rate Since the feature amount is calculated, the corrected feature amount considering the non-uniformity can be calculated even for an observation target with a non-uniform attenuation rate. Therefore, according to the present invention, it is possible to accurately identify the tissue properties of the observation target having a non-uniform attenuation rate.
  • the optimum attenuation rate for each divided region existing between the surface position of the ultrasonic transducer and the sampling point is accumulated by weighting the round-trip distance in the depth direction in each divided region. Since the cumulative attenuation rate at the sampling point is calculated by adding, the attenuation rate in the middle from the surface of the ultrasonic transducer to the sampling point can be set appropriately.
  • FIG. 12 is a diagram illustrating another example of setting divided areas in an ultrasound image.
  • regions having the same depth, that is, regions having the same value of j all have the same area.
  • the S / N ratio in the remote area can be improved by increasing the area of the divided area far from the ultrasonic transducer 21.
  • This expression (12) differs from expression (7) in that the length H j in the depth direction of the divided area Q Ij in the first term on the right side depends on the area.
  • FIG. 13 is a diagram illustrating an outline of another method for setting the optimum attenuation rate performed by the attenuation rate setting unit 334.
  • the attenuation rate setting unit 334 calculates a minimum value S ′ ( ⁇ ) min at 0 (dB / cm / MHz) ⁇ ⁇ ⁇ 1.0 (dB / cm / MHz) for the curve R, and
  • the value ⁇ ′ of the attenuation rate candidate value at that time is set as the optimum attenuation rate. Accordingly, in the case shown in FIG. 13, the optimum attenuation rate ⁇ ′ is a value between 0 (dB / cm / MHz) and 0.2 (dB / cm / MHz).
  • the attenuation rate setting unit 334 calculates an optimum attenuation rate equivalent value corresponding to the optimum attenuation rate for all frames of the ultrasonic image, and a predetermined number of optimum attenuation values including the optimum attenuation rate equivalent value in the latest frame.
  • An average value, median value, or mode value of the rate equivalent values may be set as the optimum attenuation rate. In this case, compared with the case where the optimum attenuation rate is set in each frame, the change in the optimum attenuation rate is reduced, and the value can be stabilized.
  • the attenuation rate setting unit 334 may set an optimal attenuation rate at a predetermined frame interval of the ultrasonic image. Thereby, the amount of calculation can be reduced significantly. In this case, the most recently set optimum attenuation value may be used until the next optimum attenuation rate is set.
  • the input unit 35 may be configured to accept an input of setting change of the initial value ⁇ 0 of the attenuation rate candidate value.
  • an amount giving statistical variation for example, any one of standard deviation, a difference between the maximum value and minimum value of the feature amount in the population, and a half-value width of the distribution of the feature amount can be applied.
  • a difference between the maximum value and minimum value of the feature amount in the population, and a half-value width of the distribution of the feature amount can be applied.
  • distribution is applied as an amount which gives statistical dispersion
  • the attenuation rate setting unit 334 can calculate statistical variations of a plurality of types of preliminary correction feature amounts, and can set an attenuation rate candidate value when the statistical variation is minimum as an optimal attenuation rate. It is.
  • the attenuation rate setting unit 334 may attenuate the frequency spectrum using a plurality of attenuation rate candidate values, and perform a regression analysis on the frequency spectrum after the attenuation correction to calculate the preliminary correction feature amount. Good.
  • ultrasonic probes other than the ultrasonic endoscope.
  • a thin ultrasonic miniature probe without an optical system may be applied.
  • Ultrasonic miniature probes are usually inserted into the biliary tract, bile duct, pancreatic duct, trachea, bronchi, urethra, ureter, and used to observe surrounding organs (pancreas, lung, prostate, bladder, lymph nodes, etc.).
  • an external ultrasonic probe that irradiates ultrasonic waves from the body surface of the subject may be applied as the ultrasonic probe.
  • the extracorporeal ultrasonic probe is usually used for observing an abdominal organ (liver, gallbladder, bladder), breast (particularly mammary gland), and thyroid gland.
  • the ultrasonic observation apparatus, the operation method of the ultrasonic observation apparatus, and the operation program of the ultrasonic observation apparatus according to the present invention are used to accurately distinguish the tissue characteristics of the observation target having a nonuniform attenuation rate. Useful.

Abstract

 本発明にかかる超音波観測装置は、超音波信号の周波数を解析することによって得られた複数の周波数スペクトルの特徴量をそれぞれ算出する特徴量算出部と、超音波画像を複数の領域に分割した各分割領域において、互いに異なる減衰特性を与える複数の単位長さおよび単位周波数あたりの減衰率候補値の各々を用いることにより、各周波数スペクトルの特徴量に対して超音波の影響を排除する減衰補正を施すことによって減衰率候補値ごとの各周波数スペクトルの予備補正特徴量を算出し、複数の減衰率候補値の中から観測対象に最適な減衰率を設定する減衰率設定部と、超音波振動子の表面とサンプリング点との間に存在する分割領域の最適な減衰率を用いてサンプリング点における単位周波数あたりの累積減衰率を算出し、該累積減衰率を用いて特徴量の減衰補正を行うことにより補正特徴量を算出する特徴量補正部とを備える。

Description

超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
 本発明は、超音波を用いて観測対象の組織を観測する超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラムに関する。
 従来、超音波を用いて観測対象の組織を観察する超音波観測装置において、受信信号に対して生体組織の周波数依存減衰を補償するために補正を行う技術が知られている(例えば、特許文献1を参照)。この技術では、被検体からの反射波に対して受信点の深さに応じた動的な補正処理およびパルス圧縮処理を順次行った受信信号を用いて超音波画像を形成する。
特開2010-246640号公報
 しかしながら、上述した特許文献1に記載の技術では、超音波を送受信するプローブ表面から離れた関心領域の受信信号を補正する際、その関心領域までの減衰率の変化は考慮されていない。このため、観測対象の減衰率が均一である場合にしか適用できず、減衰率が不均一な観測対象の組織性状を精度よく鑑別することは困難であった。
 本発明は、上記に鑑みてなされたものであって、減衰率が不均一である観測対象の組織性状を精度よく鑑別することを可能にする超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラムを提供することを目的とする。
 上述した課題を解決し、目的を達成するために、本発明に係る超音波観測装置は、観測対象へ超音波を送信し、該観測対象で反射された超音波を受信する超音波振動子を備えた超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波観測装置であって、前記超音波信号の周波数を解析することによって前記超音波信号の受信深度および受信方向に応じた複数の周波数スペクトルを算出する周波数解析部と、前記複数の周波数スペクトルの特徴量をそれぞれ算出する特徴量算出部と、前記超音波画像を複数の領域に分割した各分割領域において、前記超音波が前記観測対象を伝播する際の互いに異なる減衰特性を与える複数の単位長さおよび単位周波数あたりの減衰率候補値の各々を用いることにより、各周波数スペクトルの特徴量に対して前記超音波の影響を排除する減衰補正を施すことによって前記減衰率候補値ごとの前記各周波数スペクトルの予備補正特徴量を算出し、この算出結果をもとに前記複数の減衰率候補値の中から前記観測対象に最適な減衰率を設定する減衰率設定部と、前記減衰率設定部が前記分割領域ごとに設定した最適な減衰率のうち前記超音波振動子の表面とサンプリング点との間に存在する前記分割領域の最適な減衰率を用いて前記サンプリング点における単位周波数あたりの累積減衰率を算出し、該累積減衰率を用いて前記特徴量の減衰補正を行うことによって補正特徴量を算出する特徴量補正部と、を備えたことを特徴とする。
 本発明に係る超音波観測装置は、上記発明において、前記特徴量補正部は、前記超音波振動子の表面位置と前記サンプリング点との間に存在する分割領域ごとの前記最適な減衰率を各分割領域における深度方向の往復距離で重み付けしたものを累積加算することによって前記サンプリング点における前記累積減衰率を算出することを特徴とする。
 本発明に係る超音波観測装置は、上記発明において、深度方向に沿って隣接する2つの前記分割領域は、前記超音波振動子から遠い方における前記深度方向の長さが、前記超音波振動子から近い方の前記深度方向の長さ以上であることを特徴とする。
 本発明に係る超音波観測装置は、上記発明において、前記減衰率設定部は、前記減衰率候補値ごとに前記予備補正特徴量の統計的なばらつきを算出し、該統計的なばらつきが最小である減衰率候補値を前記最適な減衰率とすることを特徴とする。
 本発明に係る超音波観測装置は、上記発明において、前記補正特徴量に関する情報を前記超音波画像とともに表示する特徴量画像データを生成する特徴量画像データ生成部をさらに備えたことを特徴とする。
 本発明に係る超音波観測装置は、上記発明において、前記特徴量算出部は、前記各周波数スペクトルをn次式(nは正の整数)で近似する処理を行うことによって前記特徴量を算出することを特徴とする。
 本発明に係る超音波観測装置は、上記発明において、前記特徴量算出部は、前記周波数スペクトルにおける所定の周波数帯域を一次式で近似し、前記一次式の切片および傾き、ならびに前記周波数帯域の中間周波数における前記一次式の値であるミッドバンドフィットのうち、前記傾きおよび前記ミッドバンドフィットのいずれか一方を含む一つまたは複数を前記特徴量として算出し、前記減衰率設定部は、前記傾きおよび前記ミッドバンドフィットのいずれか一方に基づいて前記最適な減衰率を設定することを特徴とする。
 本発明に係る超音波観測装置は、上記発明において、前記減衰率設定部は、前記傾きを前記特徴量とする場合は前記傾きに基づいて前記最適な減衰率を設定し、前記ミッドバンドフィットを前記特徴量とする場合は前記ミッドバンドフィットに基づいて前記最適な減衰率を設定することを特徴とする。
 本発明に係る超音波観測装置の作動方法は、観測対象へ超音波を送信し、該観測対象で反射された超音波を受信する超音波振動子を備えた超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波観測装置の作動方法であって、周波数解析部が、前記超音波信号の周波数を解析することによって前記超音波信号の受信深度および受信方向に応じた複数の周波数スペクトルを算出する周波数解析ステップと、特徴量算出部が、前記複数の周波数スペクトルの特徴量をそれぞれ算出する特徴量算出ステップと、減衰率設定部が、前記超音波画像を複数の領域に分割した各分割領域において、前記超音波が前記観測対象を伝播する際の互いに異なる減衰特性を与える複数の単位長さおよび単位周波数あたりの減衰率候補値の各々を用いることにより、各周波数スペクトルの特徴量に対して前記超音波の影響を排除する減衰補正を施すことによって前記減衰率候補値ごとの前記各周波数スペクトルの予備補正特徴量を算出し、この算出結果をもとに前記複数の減衰率候補値の中から前記観測対象に最適な減衰率を設定する減衰率設定ステップと、特徴量補正部が、前記分割領域ごとに設定された最適な減衰率のうち前記超音波振動子の表面とサンプリング点との間に存在する前記分割領域の最適な減衰率を用いて前記サンプリング点における単位周波数あたりの累積減衰率を算出し、該累積減衰率を用いて前記特徴量の減衰補正を行うことによって補正特徴量を算出する特徴量補正ステップと、を有することを特徴とする。
 本発明に係る超音波観測装置の作動プログラムは、観測対象へ超音波を送信し、該観測対象で反射された超音波を受信する超音波振動子を備えた超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波観測装置に、周波数解析部が、前記超音波信号の周波数を解析することによって前記超音波信号の受信深度および受信方向に応じた複数の周波数スペクトルを算出する周波数解析ステップと、特徴量算出部が、前記複数の周波数スペクトルの特徴量をそれぞれ算出する特徴量算出ステップと、減衰率設定部が、前記超音波画像を複数の領域に分割した各分割領域において、前記超音波が前記観測対象を伝播する際の互いに異なる減衰特性を与える複数の単位長さおよび単位周波数あたりの減衰率候補値の各々を用いることにより、各周波数スペクトルの特徴量に対して前記超音波の影響を排除する減衰補正を施すことによって前記減衰率候補値ごとの前記各周波数スペクトルの予備補正特徴量を算出し、この算出結果をもとに前記複数の減衰率候補値の中から前記観測対象に最適な減衰率を設定する減衰率設定ステップと、特徴量補正部が、前記分割領域ごとに設定された最適な減衰率のうち前記超音波振動子の表面とサンプリング点との間に存在する前記分割領域の最適な減衰率を用いて前記サンプリング点における単位周波数あたりの累積減衰率を算出し、該累積減衰率を用いて前記特徴量の減衰補正を行うことによって補正特徴量を算出する特徴量補正ステップと、を実行させることを特徴とする。
 本発明によれば、超音波画像を複数の領域に分割した各分割領域において、互いに異なる減衰特性を与える複数の単位長さおよび単位周波数あたりの減衰率候補値に応じた周波数スペクトルの予備補正特徴量を算出し、この算出結果をもとに複数の減衰率候補値の中から観測対象に最適な減衰率を分割領域ごとに設定し、超音波振動子の表面とサンプリング点との間に存在する分割領域の最適な減衰率を用いてサンプリング点における単位周波数あたりの累積減衰率を算出し、該累積減衰率を用いて特徴量の減衰補正を行うことによって補正特徴量を算出するため、減衰率が不均一な観測対象であっても不均一さを考慮した補正特徴量を算出することができる。したがって、本発明によれば、減衰率が不均一である観測対象の組織性状を精度よく鑑別することが可能となる。
図1は、本発明の一実施の形態に係る超音波観測装置を備えた超音波診断システムの機能構成を示すブロック図である。 図2は、本発明の一実施の形態に係る超音波観測装置の信号増幅部が行う増幅処理における受信深度と増幅率との関係を示す図である。 図3は、本発明の一実施の形態に係る超音波観測装置の増幅補正部が行う増幅補正処理における受信深度と増幅率との関係を示す図である。 図4は、超音波信号の1つの音線におけるデータ配列を模式的に示す図である。 図5は、本発明の一実施の形態に係る超音波観測装置の周波数解析部が算出した周波数スペクトルの例を示す図である。 図6は、超音波画像における分割領域の設定例を模式的に示す図である。 図7は、本発明の一実施の形態に係る超音波観測装置の減衰率設定部が補正した予備補正特徴量をパラメータとして有する直線を示す図である。 図8は、同じ観測対象に対して2つの異なる減衰率候補値に基づいてそれぞれ減衰補正された予備補正特徴量の分布例を模式的に示す図である。 図9は、本発明の一実施の形態に係る超音波観測装置が行う処理の概要を示すフローチャートである。 図10は、本発明の一実施の形態に係る超音波観測装置の周波数解析部が実行する処理の概要を示すフローチャートである。 図11は、本発明の一実施の形態に係る超音波観測装置の減衰率設定部が行う処理の概要を示す図である。 図12は、超音波画像における分割領域の別な設定例を示す図である。 図13は、本発明の一実施の形態に係る超音波観測装置の減衰率設定部が行う最適な減衰率の別な設定方法の概要を示す図である。
 以下、添付図面を参照して、本発明を実施するための形態(以下、「実施の形態」という)を説明する。
 図1は、本発明の一実施の形態に係る超音波観測装置を備えた超音波診断システムの機能構成を示すブロック図である。同図に示す超音波診断システム1は、観測対象である被検体へ超音波を送信し、該被検体で反射された超音波を受信する超音波内視鏡2と、超音波内視鏡2が取得した超音波信号に基づいて超音波画像を生成する超音波観測装置3と、超音波観測装置3が生成した超音波画像を表示する表示装置4と、を備える。
 超音波内視鏡2は、その先端部に、超音波観測装置3から受信した電気的なパルス信号を超音波パルス(音響パルス)に変換して被検体へ照射するとともに、被検体で反射された超音波エコーを電圧変化で表現する電気的なエコー信号に変換して出力する超音波振動子21を有する。超音波振動子21は、コンベックス振動子、リニア振動子およびラジアル振動子のいずれでも構わない。超音波内視鏡2は、超音波振動子21をメカ的に走査させるものであってもよいし、超音波振動子21として複数の素子をアレイ状に設け、送受信にかかわる素子を電子的に切り替えたり、各素子の送受信に遅延をかけたりすることで、電子的に走査させるものであってもよい。以下、本実施の形態では、説明の便宜上、超音波振動子21がリニア振動子であるものとする。
 超音波内視鏡2は、通常は撮像光学系および撮像素子を有しており、被検体の消化管(食道、胃、十二指腸、大腸)または呼吸器(気管、気管支)へ挿入され、消化管、呼吸器やその周囲臓器(膵臓、胆嚢、胆管、胆道、リンパ節、縦隔臓器、血管等)を撮像することが可能である。また、超音波内視鏡2は、撮像時に被検体へ照射する照明光を導くライトガイドを有する。このライトガイドは、先端部が超音波内視鏡2の被検体への挿入部の先端まで達している一方、基端部が照明光を発生する光源装置に接続されている。
 超音波観測装置3は、超音波内視鏡2と電気的に接続され、所定の波形および送信タイミングに基づいて高電圧パルスからなる送信信号(パルス信号)を超音波振動子21へ送信するとともに、超音波振動子21から電気的な受信信号であるエコー信号を受信してデジタルの高周波(RF:Radio Frequency)信号のデータ(以下、RFデータという)を生成、出力する送受信部31と、送受信部31から受信したRFデータをもとにデジタルのBモード用受信データを生成する信号処理部32と、送受信部31から受信したRFデータに対して所定の演算を施す演算部33と、各種画像データを生成する画像処理部34と、キーボード、マウス、タッチパネル等のユーザインタフェースを用いて実現され、各種情報の入力を受け付ける入力部35と、超音波診断システム1全体を制御する制御部36と、超音波観測装置3の動作に必要な各種情報を記憶する記憶部37と、を備える。
 送受信部31は、エコー信号を増幅する信号増幅部311を有する。信号増幅部311は、受信深度が大きいエコー信号ほど高い増幅率で増幅するSTC(Sensitivity Time Control)補正を行う。図2は、信号増幅部311が行う増幅処理における受信深度と増幅率との関係を示す図である。図2に示す受信深度zは、超音波の受信開始時点からの経過時間に基づいて算出される量である。図2に示すように、増幅率β(dB)は、受信深度zが閾値zthより小さい場合、受信深度zの増加に伴ってβからβth(>β0)へ線型に増加する。また、増幅率β(dB)は、受信深度zが閾値zth以上である場合、一定値βthをとる。閾値zthの値は、観測対象から受信する超音波信号がほとんど減衰してしまい、ノイズが支配的になるような値である。より一般に、増幅率βは、受信深度zが閾値zthより小さい場合、受信深度zの増加に伴って単調増加すればよい。図2に示す関係は、予め記憶部37に記憶されている。
 送受信部31は、信号増幅部311によって増幅されたエコー信号に対してフィルタリング等の処理を施した後、A/D変換することによって時間ドメインのRFデータを生成し、信号処理部32および演算部33へ出力する。超音波内視鏡2が複数の素子をアレイ状に設けた超音波振動子21を電子的に走査させる構成を有する場合、送受信部31は、複数の素子に対応したビーム合成用の多チャンネル回路を有する。
 送受信部31が送信するパルス信号の周波数帯域は、超音波振動子21におけるパルス信号の超音波パルスへの電気音響変換の線型応答周波数帯域をほぼカバーする広帯域にするとよい。また、信号増幅部311におけるエコー信号の各種処理周波数帯域は、超音波振動子21による超音波エコーのエコー信号への音響電気変換の線型応答周波数帯域をほぼカバーする広帯域にするとよい。これらにより、後述する周波数スペクトルの近似処理を実行する際、精度のよい近似を行うことが可能となる。
 送受信部31は、制御部36が出力する各種制御信号を超音波内視鏡2に対して送信するとともに、超音波内視鏡2から識別用のIDを含む各種情報を受信して制御部36へ送信する機能も有する。
 信号処理部32は、RFデータに対してバンドパスフィルタ、包絡線検波、対数変換など公知の処理を施し、デジタルのBモード用受信データを生成する。対数変換では、RFデータを基準電圧で除した量の常用対数をとってデシベル値で表現する。信号処理部32は、生成したBモード用受信データを、画像処理部34へ出力する。信号処理部32は、CPU(Central Processing Unit)や各種演算回路等を用いて実現される。
 演算部33は、送受信部31が出力したRFデータに対して受信深度によらず増幅率を一定とするよう増幅補正を行う増幅補正部331と、増幅補正を行ったRFデータに高速フーリエ変換(FFT:Fast Fourier Transform)を施して周波数解析を行うことにより超音波信号の受信深度および受信方向に応じた複数の周波数スペクトルを算出する周波数解析部332と、各周波数スペクトルの特徴量を算出する特徴量算出部333と、超音波画像を複数の領域に分割したときの各分割領域において、超音波が観測対象を伝播する際の互いに異なる減衰特性を与える複数の単位長さおよび単位周波数あたりの減衰率候補値の各々を用いることにより、各周波数スペクトルの特徴量に対して超音波の影響を排除する減衰補正を施すことによって減衰率候補値ごとの各周波数スペクトルの予備補正特徴量を算出し、この算出結果をもとに複数の減衰率候補値の中から観測対象に最適な減衰率を設定する減衰率設定部334と、減衰率設定部334が分割領域ごとに設定した最適な減衰率のうち超音波振動子の表面とサンプリング点との間に存在する分割領域の最適な減衰率を用いてサンプリング点における単位周波数あたりの累積減衰率を算出し、該累積減衰率を用いて特徴量の減衰補正を行うことによって補正特徴量を算出する特徴量補正部335と、を有する。演算部33は、CPUや各種演算回路等を用いて実現される。
 図3は、増幅補正部331が行う増幅補正処理における受信深度と増幅率との関係を示す図である。図3に示すように、増幅補正部331が行う増幅処理における増幅率β(dB)は、受信深度zがゼロのとき最大値βth-β0をとり、受信深度zがゼロから閾値zthに達するまで線型に減少し、受信深度zが閾値zth以上のときゼロである。なお、図3に示す関係は、予め記憶部37に記憶されている。増幅補正部331が図3に示す関係に基づいてデジタルRF信号を増幅補正することにより、信号増幅部311におけるSTC補正の影響を相殺し、一定の増幅率βthの信号を出力することができる。なお、増幅補正部331が行う増幅補正処理における受信深度zと増幅率βの関係は、信号増幅部311が行う増幅処理における受信深度と増幅率の関係に応じて異なることは勿論である。
 このような増幅補正を行う理由を説明する。STC補正は、アナログ信号波形の振幅を全周波数帯域にわたって均一に、かつ、深度に対しては単調増加する増幅率で増幅させることで、アナログ信号波形の振幅から減衰の影響を排除する補正処理である。このため、エコー信号の振幅を輝度に変換して表示するBモード画像を生成する場合、かつ、一様な組織を走査した場合には、STC補正を行うことによって深度によらず輝度値が一定になる。すなわち、Bモード画像の輝度値から減衰の影響を排除する効果を得ることができる。
 一方、本実施の形態のように超音波の周波数スペクトルを算出して解析した結果を利用する場合、STC補正でも超音波の伝播に伴う減衰の影響を正確に排除できるとは限らない、という問題がある。なぜなら、一般に減衰量は周波数によって異なるが(後述する式(1)を参照)、STC補正の増幅率は距離だけに応じて変化し、周波数依存性がないためである。
 上述した問題を解決するには、Bモード画像を生成する際にSTC補正を施した受信信号を出力する一方、周波数スペクトルに基づいた画像を生成する際に、Bモード画像を生成するための送信とは異なる新たな送信を行い、STC補正を施していない受信信号を出力することが考えられる。ところがこの場合には、受信信号に基づいて生成される画像データのフレームレートが低下してしまうという問題がある。
 そこで、本実施の形態では、生成される画像データのフレームレートを維持しつつ、Bモード画像用にSTC補正を施した信号に対してSTC補正の影響を排除するために、増幅補正部331によって増幅率の補正を行う。
 周波数解析部332は、増幅補正部331が増幅補正した各音線のRFデータ(ラインデータ)を所定の時間間隔でサンプリングし、サンプルデータを生成する。周波数解析部332は、サンプルデータ群にFFT処理を施すことにより、RFデータ上の複数の箇所(データ位置)における周波数スペクトルを算出する。
 図4は、超音波信号の1つの音線におけるデータ配列を模式的に示す図である。同図に示す音線SRkにおいて、白または黒の長方形は、1つのサンプル点におけるデータを意味している。また、音線SRkにおいて、右側に位置するデータほど、超音波振動子21から音線SRkに沿って計った場合の深い箇所からのサンプルデータである(図4の矢印を参照)。音線SRkは、送受信部31が行うA/D変換におけるサンプリング周波数(例えば50MHz)に対応した時間間隔で離散化されている。図4では、番号kの音線SRkの8番目のデータ位置を受信深度zの方向の初期値Z(k) 0として設定した場合を示しているが、初期値の位置は任意に設定することができる。周波数解析部332による算出結果は複素数で得られ、記憶部37に格納される。
 図4に示すデータ群Fj(j=1、2、・・・、K)は、FFT処理の対象となるサンプルデータ群である。一般に、FFT処理を行うためには、サンプルデータ群が2のべき乗のデータ数を有している必要がある。この意味で、サンプルデータ群Fj(j=1、2、・・・、K-1)はデータ数が16(=24)で正常なデータ群である一方、サンプルデータ群FKは、データ数が12であるため異常なデータ群である。異常なデータ群に対してFFT処理を行う際には、不足分だけゼロデータを挿入することにより、正常なサンプルデータ群を生成する処理を行う。この点については、周波数解析部332の処理を説明する際に詳述する(図10を参照)。
 図5は、周波数解析部332が算出した周波数スペクトルの例を示す図である。ここでいう「周波数スペクトル」とは、サンプルデータ群にFFT処理を施すことによって得られた「ある受信深度zにおける強度の周波数分布」を意味する。また、ここでいう「強度」とは、例えばエコー信号の電圧、エコー信号の電力、超音波エコーの音圧、超音波エコーの音響エネルギー等のパラメータ、これらパラメータの振幅や時間積分値やその組み合わせのいずれかを指す。
 図5では、横軸に周波数fを取っている。また、図5では、縦軸に、強度I0を基準強度Ic(定数)で除した量の常用対数(デシベル表現)I=10log10(I0/Ic)を取っている。図5において、受信深度zは一定である。図5に示す直線L10については後述する。なお、本実施の形態において、曲線および直線は、離散的な点の集合からなる。
 図5に示す周波数スペクトルC1において、以後の演算に使用する周波数帯域の下限周波数fLおよび上限周波数fHは、超音波振動子21の周波数帯域、送受信部31が送信するパルス信号の周波数帯域などをもとに決定されるパラメータであり、例えばfL=3MHz、fH=10MHzである。以下、図5において、下限周波数fLおよび上限周波数fHによって定まる周波数帯域を「周波数帯域U」という。
 一般に、周波数スペクトルは、観測対象が生体組織である場合、超音波が走査された生体組織の性状によって異なる傾向を示す。これは、周波数スペクトルが、超音波を散乱する散乱体の大きさ、数密度、音響インピーダンス等と相関を有しているためである。ここでいう「生体組織の性状」とは、例えば悪性腫瘍(癌)、良性腫瘍、内分泌腫瘍、粘液性腫瘍、正常組織、嚢胞、脈管などのことである。
 特徴量算出部333は、所定周波数帯域における周波数スペクトルに対して回帰分析を行って一次式で近似することにより、この近似した一次式を特徴付ける特徴量を算出する。例えば、図5に示す周波数スペクトルC1の場合、特徴量算出部333は、周波数帯域Uで回帰分析を行うことによって近似直線L10を得る。以下、近似直線L10を周波数fの一次式I=a0f+b0で表すと、特徴量算出部333は、直線L10に対応する特徴量として、傾きa0、切片b0、および周波数帯域Uの中心周波数fM=(fs+fe)/2における強度Iの値であるミッドバンドフィット(Mid-band fit)c0=a0M+b0を算出する。なお、特徴量算出部333は、回帰分析によって二次以上の多項式で周波数スペクトルを近似するようにしてもよい。
 3つの補正前特徴量のうち、傾きa0は、超音波の散乱体の大きさと相関を有し、一般に散乱体が大きいほど傾きが小さな値を有すると考えられる。また、切片b0は、散乱体の大きさ、音響インピーダンスの差、散乱体の数密度(濃度)等と相関を有している。具体的には、切片b0は、散乱体が大きいほど大きな値を有し、音響インピーダンスの差が大きいほど大きな値を有し、散乱体の数密度が大きいほど大きな値を有すると考えられる。ミッドバンドフィットc0は、傾きa0と切片b0から導出される間接的なパラメータであり、有効な周波数帯域内の中心におけるスペクトルの強度を与える。このため、ミッドバンドフィットc0は、散乱体の大きさ、音響インピーダンスの差、散乱体の数密度に加えて、Bモード画像の輝度とある程度の相関を有していると考えられる。
 減衰率設定部334は、超音波画像を複数の領域に分割したときの各分割領域において、単位長さおよび単位周波数あたりの超音波の減衰量を与える減衰率を用いて特徴量の減衰補正を行う。図6は、超音波画像における分割領域の設定例を模式的に示す図である。図6では、矩形状をなす超音波画像101を30個の分割領域Pij(i=1~5、j=1~6)に分割した場合を示している。なお、図6の破線は音線を表している。図6に示す場合、分割領域Pijは、すべて同じ面積を有している。このため、分割領域Pijにおける超音波振動子21の表面位置102からの深度方向の長さH(以下、分割領域Pijの高さHという)は一定である。分割領域Pijに関する情報は、記憶部37が有する分割領域情報記憶部371に格納されている。なお、音線の数はあくまでも一例に過ぎず、図6に示す場合に限定されるものではない。また、図6は超音波振動子21がリニア振動子である場合の超音波画像であるため矩形状をなしているが、他のタイプの超音波振動子21の場合に超音波画像の外形状が異なることはいうまでもない。
 一般に、超音波の減衰量A(f,z)は、超音波が受信深度0と受信深度zとの間を往復する間に生じる減衰であり、往復する前後の強度変化(デシベル表現での差)として定義される。減衰量A(f,z)は、一様な組織内では周波数に比例することが経験的に知られており、以下の式(1)で表現される。
  A(f,z)=2αzf  ・・・(1)
ここで、比例定数αは減衰率と呼ばれる量であり、単位長さおよび単位周波数あたりの超音波の減衰量を与える。また、zは超音波の受信深度であり、fは周波数である。減衰率αの具体的な値は、観測対象が生体である場合、生体の部位に応じて定まる。減衰率αの単位は、例えばdB/cm/MHzである。
 減衰率設定部334は、複数の減衰率候補値の中から最適な減衰率を設定する。この際、減衰率設定部334は、減衰率候補値αを用いて、特徴量算出部333が算出した特徴量(傾きa0、切片b0、ミッドバンドフィットc0)に対し、以下に示す式(2)~(4)にしたがって減衰補正を行うことにより、予備補正特徴量a、b、cを算出する。
  a=a0+2αz  ・・・(2)
  b=b0  ・・・(3)
  c=c0+A(fM,z)=c0+2αzfM(=afM+b) ・・・(4)
式(2)、(4)からも明らかなように、減衰率設定部334は、超音波の受信深度zが大きいほど、補正量が大きい補正を行う。また、式(3)によれば、切片に関する補正は恒等変換である。これは、切片が周波数0(Hz)に対応する周波数成分であって減衰の影響を受けないためである。
 図7は、減衰率設定部334が補正した予備補正特徴量a、b、cをパラメータとして有する直線を示す図である。直線L1の式は、
  I=af+b=(a0+2αz)f+b0  ・・・(5)
で表される。この式(5)からも明らかなように、直線L1は、減衰補正前の直線L10と比較して、傾きが大きく(a>a0)、かつ切片が同じ(b=b0)である。
 減衰率設定部334は、各分割領域において、減衰率候補値ごとに算出した予備補正特徴量の統計的なばらつきが最小である減衰率候補値を最適な減衰率として設定する。本実施の形態では、統計的なばらつきを示す量として分散を適用する。この場合、減衰率設定部334は、分散が最小となる減衰率候補値を最適な減衰率として設定する。上述した3つの予備補正特徴量a、b、cのうち独立なのは2つである。加えて、予備補正特徴量bは減衰率に依存しない。したがって、予備補正特徴量a、cに対して最適な減衰率を設定する場合、減衰率設定部334は、予備補正特徴量aおよびcのいずれか一方の分散を算出すればよい。
 ただし、減衰率設定部334が、予備補正特徴量aを用いて最適な減衰率を設定する場合は予備補正特徴量aの分散を適用し、予備補正特徴量cを用いて最適な減衰率を設定する場合は予備補正特徴量cの分散を適用するのがより好ましい。これは、減衰量A(f,z)を与える式(1)があくまで理想的なものに過ぎず、現実には以下の式(6)の方が適切であることによる。
  A(f,z)=2αzf+2α1z  ・・・(6)
式(6)の右辺第2項のα1は、超音波の受信深度zに比例して信号強度が変化する大きさを表す係数であり、観測対象の組織が不均一であることや、ビーム合成時のチャンネル数の変更などに起因して発生する信号強度の変化を表す係数である。式(6)の右辺第2項が存在するため、予備補正特徴量cを用いて最適な減衰率を設定する場合は、予備補正特徴量cの分散を適用した方が正確に減衰を補正することができる(式(4)を参照)。一方、周波数fに比例する係数である予備補正特徴量aを用いて最適な減衰率を設定する場合は、予備補正特徴量aの分散を適用した方が、式(6)の右辺第2項の影響を排除して正確に減衰を補正することができる。
 ここで、統計的なばらつきに基づいて最適な減衰率を設定することができる理由を説明する。観測対象に最適な減衰率を適用した場合、観測対象と超音波振動子21との距離に関わらず、特徴量は観測対象に固有の値へ収束し、統計的なばらつきが小さくなると考えられる。その一方で、観測対象に適合しない減衰率候補値を最適な減衰率とした場合、減衰補正が過剰であるかまたは不足するため、超音波振動子21との距離に応じて特徴量にずれが生じ、特徴量の統計的なばらつきが大きくなると考えられる。したがって、統計的なばらつきが最も小さい減衰率候補値が、観測対象にとって最適な減衰率であるということができる。
 図8は、同じ観測対象に対して2つの異なる減衰率候補値に基づいてそれぞれ減衰補正された予備補正特徴量の分布例を模式的に示す図である。図8では、横軸を予備補正特徴量とし、縦軸を頻度としている。図8に示す2つの分布曲線N1、N2は、頻度の総和が同じである。図8に示す場合、分布曲線N1は、分布曲線N2と比較して特徴量の統計的なばらつきが小さく(分散が小さく)、山が急峻な形状をなす。したがって、減衰率設定部334は、この2つの分布曲線N1、N2に対応する2つの減衰率候補値から最適な減衰率を設定する場合、分布曲線N1に対応する減衰率候補値を最適な減衰率として設定する。
 特徴量補正部335は、分割領域ごとの最適な減衰率を用いてサンプリング点における単位周波数あたりの累積減衰率(以下、単に累積減衰率ともいう)を算出し、該累積減衰率を用いて特徴量の減衰補正を行う。任意のサンプリング点における累積減衰率は、超音波振動子21の表面からの距離とその間に存在する分割領域における最適な減衰率を用いて算出される。分割領域PIJ(1≦I≦imax、1≦J≦jmax;I,Jは定数)で超音波振動子21に近い側の境界からの距離がh(0<h≦H)であるようなサンプリング点SIJ(h)における累積減衰率γIJ(h)は、
  γIJ(h)=[Σj=1,2,…,J-1(2H・α(PIj))]+2h・α(PIJ
                              ・・・(7)
と表される。ここで、Σj=1,2,…,J-1はj=1~J-1の和を取ることを意味し、α(PIj)は分割領域PIjにおける最適な減衰率を表している。式(7)の右辺の第1項は、分割領域ごとの超音波の往復距離2Hに対して最適な減衰率α(PIj)を乗じたものの和である。また、式(7)の右辺の第2項は、分割領域PIJにおけるサンプリング点SIJ(h)までの超音波の往復距離2hに分割領域PIJにおける最適な減衰率α(PIJ)を乗じたものである。このようにして、特徴量補正部335は、超音波振動子21の表面からの減衰率を累積していくことによって累積減衰率γIJ(h)を算出する。最適な減衰率の単位をdB/cm/MHzとする場合、累積減衰率の単位はdB/MHzである。
 例えば図6では、式(7)でI=2、J=3(ただしimax=5、jmax=6)の場合を例示している。この場合、サンプリング点S23(h)における累積減衰率γ23(h)は、式(7)より、
  γ23(h)=2H・α(P21)+2H・α(P22)+2h・α(P23
                                                        ・・・(8)
となる。
 特徴量補正部335は、累積減衰率γIJ(h)を用いてサンプリング点SIJ(h)における特徴量を以下のように減衰補正する。
  aIJ(h)=a0+2γIJ(h)  ・・・(9)
  bIJ(h)=b0  ・・・(10)
  cIJ(h)=c0+2fMγIJ(h)  ・・・(11)
 画像処理部34は、エコー信号の振幅を輝度に変換して表示する超音波画像であるBモード画像データを生成するBモード画像データ生成部341と、特徴量算出部333が算出した特徴量に関する情報を表示する特徴量画像データを生成する特徴量画像データ生成部342と、を有する。
 Bモード画像データ生成部341は、信号処理部32から受信したBモード用受信データに対してゲイン処理、コントラスト処理等の公知の技術を用いた信号処理を行うとともに、表示装置4における画像の表示レンジに応じて定まるデータステップ幅に応じたデータの間引き等を行うことによってBモード画像データを生成する。Bモード画像は、色空間としてRGB表色系を採用した場合の変数であるR(赤)、G(緑)、B(青)の値を一致させたグレースケール画像である。
 Bモード画像データ生成部341は、Bモード用受信データに走査範囲を空間的に正しく表現できるよう並べ直す座標変換を施した後、Bモード用受信データ間の補間処理を施すことによってBモード用受信データ間の空隙を埋め、Bモード画像データを生成する。Bモード画像データ生成部341は、生成したBモード画像データを特徴量画像データ生成部342へ出力する。
 特徴量画像データ生成部342は、特徴量算出部333が算出した特徴量に関連する視覚情報をBモード画像データにおける画像の各画素に対して重畳することによって特徴量画像データを生成する。特徴量画像データ生成部342は、例えば図4に示す1つの振幅データ群Fj(j=1、2、・・・、K)のデータ量に対応する画素領域に対し、その振幅データ群Fjから算出される周波数スペクトルの特徴量に対応する視覚情報を割り当てる。特徴量画像データ生成部342は、例えば上述した傾き、切片、ミッドバンドフィットのいずれか一つに視覚情報としての色相を対応付けることによって特徴量画像データを生成する。なお、特徴量画像データ生成部342が、傾き、切片、ミッドバンドフィットから選択される2つの特徴量の一方に色相を対応付けるとともに、他方に明暗を対応付けることによって特徴量画像データを生成するようにしてもよい。特徴量に関連する視覚情報としては、色相や明暗(明度)のほか、例えば彩度、輝度値、R(赤)、G(緑)、B(青)などの所定の表色系を構成する色空間の変数を挙げることができる。
 制御部36は、演算および制御機能を有するCPU(Central Processing Unit)や各種演算回路等を用いて実現される。制御部36は、記憶部37が記憶、格納する情報を記憶部37から読み出し、超音波観測装置3の作動方法に関連した各種演算処理を実行することによって超音波観測装置3を統括して制御する。なお、制御部36を信号処理部32および演算部33と共通のCPU等を用いて構成することも可能である。
 記憶部37は、分割領域に関する情報を記憶する分割領域情報記憶部371と、周波数解析部332が算出した周波数スペクトルの情報を受信深度および受信方向とともに記憶するスペクトル情報記憶部372と、特徴量算出部333が算出する特徴量および特徴量補正部335が補正した補正特徴量に関する情報を記憶する特徴量情報記憶部373と、減衰率設定部334が分割領域ごとに設定する最適な減衰率および特徴量補正部335が算出するサンプリング点毎の累積減衰率に関する情報を記憶する減衰率情報記憶部374と、を有する。
 記憶部37は、上記以外にも、例えば増幅処理に必要な情報(図2に示す増幅率と受信深度との関係)、増幅補正処理に必要な情報(図3に示す増幅率と受信深度との関係)、減衰補正処理に必要な情報(式(1)参照)、周波数解析処理に必要な窓関数(Hamming、Hanning、Blackman等)の情報等を記憶する。
 また、記憶部37は、超音波観測装置3の作動方法を実行するための作動プログラムを含む各種プログラムを記憶する。作動プログラムは、ハードディスク、フラッシュメモリ、CD-ROM、DVD-ROM、フレキシブルディスク等のコンピュータ読み取り可能な記録媒体に記録して広く流通させることも可能である。なお、上述した各種プログラムは、通信ネットワークを介してダウンロードすることによって取得することも可能である。ここでいう通信ネットワークは、例えば既存の公衆回線網、LAN(Local Area Network)、WAN(Wide Area Network)などによって実現されるものであり、有線、無線を問わない。
 以上の構成を有する記憶部37は、各種プログラム等が予めインストールされたROM(Read Only Memory)、および各処理の演算パラメータやデータ等を記憶するRAM(Random Access Memory)等を用いて実現される。
 図9は、以上の構成を有する超音波観測装置3が行う処理の概要を示すフローチャートである。具体的には、超音波観測装置3が超音波内視鏡2からエコー信号を受信する以降の処理の概要を示すフローチャートである。以下、図9を参照して、超音波観測装置3が行う処理を説明する。まず、超音波観測装置3は、超音波内視鏡2から超音波振動子21による観測対象の測定結果としてのエコー信号を受信する(ステップS1)。
 超音波振動子21からエコー信号を受信した信号増幅部311は、そのエコー信号の増幅を行う(ステップS2)。ここで、信号増幅部311は、例えば図2に示す増幅率と受信深度との関係に基づいてエコー信号の増幅(STC補正)を行う。
 続いて、Bモード画像データ生成部341は、信号増幅部311が増幅したエコー信号を用いてBモード画像データを生成して、表示装置4へ出力する(ステップS3)。Bモード画像データを受信した表示装置4は、そのBモード画像データに対応するBモード画像を表示する。
 増幅補正部331は、送受信部31から出力されたRFデータに対して受信深度によらず増幅率が一定となるように増幅補正を行う(ステップS4)。ここで、増幅補正部331は、例えば図3に示す増幅率と受信深度との関係に基づいて増幅補正を行う。
 この後、周波数解析部332は、増幅補正後の各音線のRFデータに対してFFTによる周波数解析を行うことによって全てのサンプルデータ群に対する周波数スペクトルを算出し、スペクトル情報記憶部372へ格納する(ステップS5)。図10は、ステップS5において周波数解析部332が行う処理の概要を示すフローチャートである。以下、図10に示すフローチャートを参照して、周波数解析処理を詳細に説明する。
 まず、周波数解析部332は、解析対象の音線を識別するカウンタkをk0とする(ステップS21)。
 続いて、周波数解析部332は、FFT演算用に生成した一連のデータ群(サンプルデータ群)を代表するデータ位置(受信深度に相当)Z(k)の初期値Z(k) 0を設定する(ステップS22)。例えば、図4では、上述したように、音線SRkの8番目のデータ位置を初期値Z(k) 0として設定した場合を示している。
 その後、周波数解析部332は、サンプルデータ群を取得し(ステップS23)、取得したサンプルデータ群に対し、記憶部37が記憶する窓関数を作用させる(ステップS24)。このようにサンプルデータ群に対して窓関数を作用させることにより、サンプルデータ群が境界で不連続になることを回避し、アーチファクトが発生するのを防止することができる。
 続いて、周波数解析部332は、データ位置Z(k)のサンプルデータ群が正常なデータ群であるか否かを判定する(ステップS25)。図4を参照した際に説明したように、サンプルデータ群は、2のべき乗のデータ数を有している必要がある。以下、正常なサンプルデータ群のデータ数を2n(nは正の整数)とする。本実施の形態では、データ位置Z(k)が、できるだけZ(k)が属するサンプルデータ群の中心になるよう設定される。具体的には、サンプルデータ群のデータ数は2nであるので、Z(k)はそのサンプルデータ群の中心に近い2n/2(=2n-1)番目の位置に設定される。この場合、サンプルデータ群が正常であるとは、データ位置Z(k)より浅部側に2n-1-1(=Nとする)個のデータがあり、データ位置Z(k)より深部側に2n-1(=Mとする)個のデータがあることを意味する。図4に示す場合、サンプルデータ群Fj(j=1、2、・・・、K-1)は正常である。なお、図4ではn=4(N=7,M=8)の場合を例示している。
 ステップS25における判定の結果、データ位置Z(k)のサンプルデータ群が正常である場合(ステップS25:Yes)、周波数解析部332は、後述するステップS27へ移行する。
 ステップS25における判定の結果、データ位置Z(k)のサンプルデータ群が正常でない場合(ステップS25:No)、周波数解析部332は、不足分だけゼロデータを挿入することによって正常なサンプルデータ群を生成する(ステップS26)。ステップS25において正常でないと判定されたサンプルデータ群(例えば図4のサンプルデータ群FK)は、ゼロデータを追加する前に窓関数が作用されている。このため、サンプルデータ群にゼロデータを挿入してもデータの不連続は生じない。ステップS26の後、周波数解析部332は、後述するステップS27へ移行する。
 ステップS27において、周波数解析部332は、サンプルデータ群を用いてFFT演算を行うことにより、振幅の周波数分布である周波数スペクトルを得る(ステップS27)。
 続いて、周波数解析部332は、データ位置Z(k)をステップ幅Dで変化させる(ステップS28)。ステップ幅Dは、記憶部37が予め記憶しているものとする。図4では、D=15の場合を例示している。ステップ幅Dは、Bモード画像データ生成部341がBモード画像データを生成する際に利用するデータステップ幅と一致させることが望ましいが、周波数解析部332における演算量を削減したい場合には、ステップ幅Dとしてデータステップ幅より大きい値を設定してもよい。
 その後、周波数解析部332は、データ位置Z(k)が音線SRkにおける最大値Z(k) maxより大きいか否かを判定する(ステップS29)。データ位置Z(k)が最大値Z(k) maxより大きい場合(ステップS29:Yes)、周波数解析部332はカウンタkを1増加させる(ステップS30)。これは、処理をとなりの音線へ移すことを意味する。一方、データ位置Z(k)が最大値Z(k) max以下である場合(ステップS29:No)、周波数解析部332はステップS23へ戻る。
 ステップS30の後、周波数解析部332は、カウンタkが最大値kmaxより大きいか否かを判定する(ステップS31)。カウンタkがkmaxより大きい場合(ステップS31:Yes)、周波数解析部332は一連の周波数解析処理を終了する。一方、カウンタkがkmax以下である場合(ステップS31:No)、周波数解析部332はステップS22に戻る。この最大値kmaxは、術者等のユーザが入力部35を通じて任意に指示入力した値、もしくは、記憶部37にあらかじめ設定された値とする。
 このようにして、周波数解析部332は、解析対象領域内の(kmax-k0+1)本の音線の各々について複数回のFFT演算を行う。FFT演算の結果として得られる周波数スペクトルは、受信深度、受信方向とともにスペクトル情報記憶部372に格納される。
 なお、以上の説明では、周波数解析部332が超音波信号を受信したすべての領域に対して周波数解析処理を行うものとしたが、入力部35が特定の深度幅および音線幅で区切られる部分領域の設定入力を受け付け可能な構成とし、設定された部分領域内においてのみ周波数解析処理を行うようにすることも可能である。
 以上説明したステップS5の周波数解析処理に続いて、特徴量算出部333は、各分割領域に含まれるサンプリング点で周波数スペクトルの特徴量を算出する(ステップS6)。具体的には、特徴量算出部333は、所定の周波数帯域の周波数スペクトルに対して回帰分析を行うことによって一次式I=a0f+b0で近似し、特徴量として傾きa0、切片b0、ミッドバンドフィットc0を算出する。例えば、図5に示す直線L10は、特徴量算出部333が周波数帯域Uの周波数スペクトルC1に対し回帰分析によって近似した回帰直線である。
 この後、減衰率設定部334は、後述する減衰補正を行う際に適用する減衰率候補値αの値を所定の初期値α0に設定する(ステップS7)。この初期値α0の値は、予め減衰率情報記憶部374が記憶しておくようにすればよい。
 続いて、減衰率設定部334は、特徴量算出部333が各周波数スペクトルに対して近似した特徴量に対し、減衰率候補値をαとして減衰補正を行うことにより、予備補正特徴量を算出し、減衰率候補値αとともに特徴量情報記憶部373に格納する(ステップS8)。図7に示す直線L1は、減衰率設定部334が減衰補正処理を行うことによって得られる直線の例である。
 ステップS8において、減衰率設定部334は、上述した式(2)、(4)における受信深度zに、超音波信号の音線のデータ配列を用いて得られるデータ位置Z=(fsp/2vs)Dnを代入することによって算出する。ここで、fspはデータのサンプリング周波数、vsは音速、Dはデータステップ幅、nは処理対象の振幅データ群のデータ位置までの音線の1番目のデータからのデータステップ数である。例えば、データのサンプリング周波数fspを50MHzとし、音速vsを1530m/secとし、図4に示すデータ配列を採用してステップ幅Dを15とすると、z=0.2295n(mm)となる。
 減衰率設定部334は、減衰率設定部334が各周波数スペクトルに対して減衰補正することによって得られた複数の予備補正特徴量のうち代表となる予備補正特徴量の分散を算出し、減衰率候補値αと対応づけて特徴量情報記憶部373へ格納する(ステップS9)。例えば、予備補正特徴量が傾きa、ミッドバンドフィットcである場合、上述したように、減衰率設定部334は、予備補正特徴量aおよびcのいずれか一方の分散を算出する。このステップS9において、予備補正特徴量aを用いて特徴量画像を生成する場合は各分割領域における予備補正特徴量aの分散を適用し、予備補正特徴量cを用いて特徴量画像を生成する場合は各分割領域における予備補正特徴量cの分散を適用するのが好ましい。
 この後、減衰率設定部334は、減衰率候補値αの値をΔαだけ増加させ(ステップS10)、増加後の減衰率候補値αと所定の最大値αmaxとの大小を比較する(ステップS11)。ステップS11における比較の結果、減衰率候補値αが最大値αmaxより大きい場合(ステップS11:Yes)、超音波観測装置3はステップS12へ移行する。一方、ステップS11における比較の結果、減衰率候補値αが最大値αmax以下である場合(ステップS11:No)、超音波観測装置3はステップS8へ戻る。
 ステップS12において、減衰率設定部334は、分割領域ごとに、特徴量情報記憶部373が記憶する減衰率候補値ごとの予備補正特徴量の分散を参照し、該分散が最小である減衰率候補値を分割領域における最適な減衰率として設定する(ステップS12)。
 図11は、減衰率設定部334が行う処理の概要を示す図である。α0=0(dB/cm/MHz)、αmax=1.0(dB/cm/MHz)、Δα=0.2(dB/cm/MHz)とした場合の減衰率候補値αと分散S(α)との関係の例を示す図である。図11に示す場合、減衰率候補値αが0.2(dB/cm/MHz)のときに分散が最小値S(α)minをとる。したがって、図11に示す場合、減衰率設定部334は、α=0.2(dB/cm/MHz)を最適な減衰率として設定する。
 この後、特徴量補正部335は、減衰率設定部334が各分割領域で設定した最適な減衰率を用いて、サンプリング点における累積減衰率を算出する(ステップS13)。例えば、図6に示す分割領域PIJのサンプリング点SIJ(h)における累積減衰率γIJ(h)は、式(7)で表される。
 続いて、特徴量補正部335は、累積減衰率を用いて特徴量を減衰補正することによって補正特徴量を算出する(ステップS14)。例えば、図6に示す分割領域PIJのサンプリング点SIJ(h)における傾きa0、切片b0、ミッドバンドフィットc0の補正特徴量aIJ(h)、bIJ(h)、cIJ(h)は、式(9)~(11)を用いてそれぞれ算出される。
 特徴量画像データ生成部342は、Bモード画像データ生成部341が生成したBモード画像データにおける各画素に対し、ステップS14で算出された補正特徴量と関連づけた視覚情報(例えば色相)を重畳することによって特徴量画像データを生成する(ステップS15)。特徴量画像データ生成部342は、生成した特徴量画像データを表示装置4へ送信する。特徴量画像データを受信した表示装置4は、受信した特徴量画像データに対応する特徴量画像を表示する。
 ステップS15の後、超音波観測装置3は一連の処理を終了する。なお、超音波観測装置3は、ステップS1~S15の処理を周期的に繰り返し実行する。
 以上説明した本発明の一実施の形態によれば、超音波画像を複数の領域に分割した各分割領域において、互いに異なる減衰特性を与える複数の単位長さおよび単位周波数あたりの減衰率候補値に応じた周波数スペクトルの予備補正特徴量を算出し、この算出結果をもとに複数の減衰率候補値の中から観測対象に最適な減衰率を分割領域ごとに設定し、超音波振動子の表面とサンプリング点との間に存在する分割領域の最適な減衰率を用いてサンプリング点における単位周波数あたりの累積減衰率を算出し、該累積減衰率を用いて特徴量の減衰補正を行うことによって補正特徴量を算出するため、減衰率が不均一な観測対象であっても不均一さを考慮した補正特徴量を算出することができる。したがって、本発明によれば、減衰率が不均一である観測対象の組織性状を精度よく鑑別することが可能となる。
 また、本実施の形態によれば、超音波振動子の表面位置とサンプリング点との間に存在する分割領域ごとの最適な減衰率を各分割領域における深度方向の往復距離で重み付けしたものを累積加算することによってサンプリング点における累積減衰率を算出するため、超音波振動子の表面からサンプリング点に至るまでの途中の減衰率を適切に設定することができる。
 ここまで、本発明を実施するための形態を説明してきたが、本発明は上述した実施の形態によってのみ限定されるべきものではない。例えば、超音波画像の分割領域の設定方法は、図6に示すものに限られるわけではない。図12は、超音波画像における分割領域の別な設定例を示す図である。図12では、超音波画像201を30個の分割領域Qij(i=1~5、j=1~6)に分割した場合を示している。分割領域Qijは、深度が同じ領域すなわちjの値が共通な領域はすべて同じ面積を有している。一方、分割領域Qijは、深度が異なる領域では、超音波振動子21の表面位置202からの深度が大きいほど高さHj(j=1~6)が大きくなっている(H1<H2<・・・<H6)。一般に、超音波の減衰は深度が大きいほど大きくなるので、超音波振動子21から遠方の分割領域の面積を大きくすることにより、遠方におけるS/N比を向上させることができる。
 この場合、分割領域QIJ(1≦I≦imax,1≦J≦jmax;I,Jは定数)で超音波振動子21に近い側の境界からの距離がh(0<h≦HJ)であるようなサンプリング点SIJ(h)における累積減衰率γIJ(h)は、
 γIJ(h)=[Σj=1,2,…,J-1(2Hj・α(QIj))]+2h・α(QIJ
                           ・・・(12)
と表される。この式(12)は、式(7)と比較して右辺第1項の分割領域QIjの深度方向の長さHjが領域に依存している点が異なる。
 また、減衰率設定部334が行う最適な減衰率の設定方法は上述した方法に限定されるわけではない。図13は、減衰率設定部334が行う最適な減衰率の別な設定方法の概要を示す図である。図13では、α0=0(dB/cm/MHz)、αmax=1.0(dB/cm/MHz)、Δα=0.2(dB/cm/MHz)とした場合の減衰率候補値αと分散S(α)との関係の例を示しており、減衰率候補値α=0、0.2、0.4、0.6、0.8、1.0(いずれもdB/cm/MHz)における分散S(α)の値は、図11とそれぞれ同じである。この場合には、減衰率設定部334が最適な減衰率を設定する前に、特徴量算出部333が回帰分析を行うことによって減衰率候補値αにおける分散S(α)の値を補間する曲線Rを算出する。その後、減衰率設定部334は、この曲線Rに対し、0(dB/cm/MHz)≦α≦1.0(dB/cm/MHz)における最小値S’(α)minを算出し、そのときの減衰率候補値の値α’を最適な減衰率として設定する。したがって、図13に示す場合、最適な減衰率α’は、0(dB/cm/MHz)と0.2(dB/cm/MHz)の間の値となる。
 また、減衰率設定部334は、超音波画像の全てのフレームで最適な減衰率に相当する最適減衰率相当値をそれぞれ算出し、最新のフレームにおける最適減衰率相当値を含む所定数の最適減衰率相当値の平均値、中央値または最頻値を最適な減衰率として設定してもよい。この場合には、各フレームで最適な減衰率を設定する場合と比較して、最適な減衰率の変化が少なくなってその値を安定させることができる。
 また、減衰率設定部334は、超音波画像の所定のフレーム間隔で最適な減衰率を設定するようにしてもよい。これにより、計算量を大幅に削減することができる。この場合には、次に最適な減衰率を設定するまでの間、最後に設定した最適な減衰率の値を使用すればよい。
 また、入力部35が減衰率候補値の初期値α0の設定変更の入力を受け付け可能な構成としてもよい。
 また、統計的なばらつきを与える量として、例えば標準偏差、母集団における特徴量の最大値と最小値の差、特徴量の分布の半値幅のいずれかを適用することも可能である。なお、統計的なばらつきを与える量として分散の逆数を適用する場合も考えられるが、この場合には、その値が最大となる減衰率候補値が最適な減衰率となることはいうまでもない。
 また、減衰率設定部334が複数種類の予備補正特徴量の統計的なばらつきをそれぞれ算出し、統計的なばらつきが最小である場合の減衰率候補値を最適な減衰率として設定することも可能である。
 また、減衰率設定部334が複数の減衰率候補値を用いて周波数スペクトルを減衰補正し、減衰補正後の周波数スペクトルに対して回帰分析を行うことによって予備補正特徴量を算出するようにしてもよい。
 また、超音波内視鏡以外の超音波プローブに対しても適用することが可能である。超音波プローブとして、例えば光学系のない細径の超音波ミニチュアプローブを適用してもよい。超音波ミニチュアプローブは、通常、胆道、胆管、膵管、気管、気管支、尿道、尿管へ挿入され、その周囲臓器(膵臓、肺、前立腺、膀胱、リンパ節等)を観察する際に用いられる。また、超音波プローブとして、被検体の体表から超音波を照射する体外式超音波プローブを適用してもよい。体外式超音波プローブは、通常、腹部臓器(肝臓、胆嚢、膀胱)、乳房(特に乳腺)、甲状腺を観察する際に用いられる。
 このように、本発明は、請求の範囲に記載した技術的思想を逸脱しない範囲内において、様々な実施の形態を含みうるものである。
 以上のように、本発明にかかる超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラムは、減衰率が不均一である観測対象の組織性状を精度よく鑑別するのに有用である。
 1 超音波診断システム
 2 超音波内視鏡
 3 超音波観測装置
 4 表示装置
 21 超音波振動子
 31 送受信部
 32 信号処理部
 33 演算部
 34 画像処理部
 35 入力部
 36 制御部
 37 記憶部
 101、201 超音波画像
 102、202 表面位置
 311 信号増幅部
 331 増幅補正部
 332 周波数解析部
 333 特徴量算出部
 334 減衰率設定部
 335 特徴量補正部
 341 Bモード画像データ生成部
 342 特徴量画像データ生成部
 371 分割領域情報記憶部
 372 スペクトル情報記憶部
 373 特徴量情報記憶部
 374 減衰率情報記憶部

Claims (10)

  1.  観測対象へ超音波を送信し、該観測対象で反射された超音波を受信する超音波振動子を備えた超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波観測装置であって、
     前記超音波信号の周波数を解析することによって前記超音波信号の受信深度および受信方向に応じた複数の周波数スペクトルを算出する周波数解析部と、
     前記複数の周波数スペクトルの特徴量をそれぞれ算出する特徴量算出部と、
     前記超音波画像を複数の領域に分割した各分割領域において、前記超音波が前記観測対象を伝播する際の互いに異なる減衰特性を与える複数の単位長さおよび単位周波数あたりの減衰率候補値の各々を用いることにより、各周波数スペクトルの特徴量に対して前記超音波の影響を排除する減衰補正を施すことによって前記減衰率候補値ごとの前記各周波数スペクトルの予備補正特徴量を算出し、この算出結果をもとに前記複数の減衰率候補値の中から前記観測対象に最適な減衰率を設定する減衰率設定部と、
     前記減衰率設定部が前記分割領域ごとに設定した最適な減衰率のうち前記超音波振動子の表面とサンプリング点との間に存在する前記分割領域の最適な減衰率を用いて前記サンプリング点における単位周波数あたりの累積減衰率を算出し、該累積減衰率を用いて前記特徴量の減衰補正を行うことによって補正特徴量を算出する特徴量補正部と、
     を備えたことを特徴とする超音波観測装置。
  2.  前記特徴量補正部は、
     前記超音波振動子の表面位置と前記サンプリング点との間に存在する分割領域ごとの前記最適な減衰率を各分割領域における深度方向の往復距離で重み付けしたものを累積加算することによって前記サンプリング点における前記累積減衰率を算出することを特徴とする請求項1に記載の超音波観測装置。
  3.  深度方向に沿って隣接する2つの前記分割領域は、前記超音波振動子から遠い方における前記深度方向の長さが、前記超音波振動子から近い方の前記深度方向の長さ以上であることを特徴とする請求項1に記載の超音波観測装置。
  4.  前記減衰率設定部は、
     前記減衰率候補値ごとに前記予備補正特徴量の統計的なばらつきを算出し、該統計的なばらつきが最小である減衰率候補値を前記最適な減衰率とすることを特徴とする請求項1に記載の超音波観測装置。
  5.  前記補正特徴量に関する情報を前記超音波画像とともに表示する特徴量画像データを生成する特徴量画像データ生成部をさらに備えたことを特徴とする請求項1に記載の超音波観測装置。
  6.  前記特徴量算出部は、
     前記各周波数スペクトルをn次式(nは正の整数)で近似する処理を行うことによって前記特徴量を算出することを特徴とする請求項1に記載の超音波観測装置。
  7.  前記特徴量算出部は、
     前記周波数スペクトルにおける所定の周波数帯域を一次式で近似し、前記一次式の切片および傾き、ならびに前記周波数帯域の中間周波数における前記一次式の値であるミッドバンドフィットのうち、前記傾きおよび前記ミッドバンドフィットのいずれか一方を含む一つまたは複数を前記特徴量として算出し、
     前記減衰率設定部は、
     前記傾きおよび前記ミッドバンドフィットのいずれか一方に基づいて前記最適な減衰率を設定することを特徴とする請求項6に記載の超音波観測装置。
  8.  前記減衰率設定部は、
     前記傾きを前記特徴量とする場合は前記傾きに基づいて前記最適な減衰率を設定し、前記ミッドバンドフィットを前記特徴量とする場合は前記ミッドバンドフィットに基づいて前記最適な減衰率を設定することを特徴とする請求項7に記載の超音波観測装置。
  9.  観測対象へ超音波を送信し、該観測対象で反射された超音波を受信する超音波振動子を備えた超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波観測装置の作動方法であって、
     周波数解析部が、前記超音波信号の周波数を解析することによって前記超音波信号の受信深度および受信方向に応じた複数の周波数スペクトルを算出する周波数解析ステップと、
     特徴量算出部が、前記複数の周波数スペクトルの特徴量をそれぞれ算出する特徴量算出ステップと、
     減衰率設定部が、前記超音波画像を複数の領域に分割した各分割領域において、前記超音波が前記観測対象を伝播する際の互いに異なる減衰特性を与える複数の単位長さおよび単位周波数あたりの減衰率候補値の各々を用いることにより、各周波数スペクトルの特徴量に対して前記超音波の影響を排除する減衰補正を施すことによって前記減衰率候補値ごとの前記各周波数スペクトルの予備補正特徴量を算出し、この算出結果をもとに前記複数の減衰率候補値の中から前記観測対象に最適な減衰率を設定する減衰率設定ステップと、
     特徴量補正部が、前記分割領域ごとに設定された最適な減衰率のうち前記超音波振動子の表面とサンプリング点との間に存在する前記分割領域の最適な減衰率を用いて前記サンプリング点における単位周波数あたりの累積減衰率を算出し、該累積減衰率を用いて前記特徴量の減衰補正を行うことによって補正特徴量を算出する特徴量補正ステップと、
     を有することを特徴とする超音波観測装置の作動方法。
  10.  観測対象へ超音波を送信し、該観測対象で反射された超音波を受信する超音波振動子を備えた超音波プローブが取得した超音波信号に基づいて超音波画像を生成する超音波観測装置に、
     周波数解析部が、前記超音波信号の周波数を解析することによって前記超音波信号の受信深度および受信方向に応じた複数の周波数スペクトルを算出する周波数解析ステップと、
     特徴量算出部が、前記複数の周波数スペクトルの特徴量をそれぞれ算出する特徴量算出ステップと、
     減衰率設定部が、前記超音波画像を複数の領域に分割した各分割領域において、前記超音波が前記観測対象を伝播する際の互いに異なる減衰特性を与える複数の単位長さおよび単位周波数あたりの減衰率候補値の各々を用いることにより、各周波数スペクトルの特徴量に対して前記超音波の影響を排除する減衰補正を施すことによって前記減衰率候補値ごとの前記各周波数スペクトルの予備補正特徴量を算出し、この算出結果をもとに前記複数の減衰率候補値の中から前記観測対象に最適な減衰率を設定する減衰率設定ステップと、
     特徴量補正部が、前記分割領域ごとに設定された最適な減衰率のうち前記超音波振動子の表面とサンプリング点との間に存在する前記分割領域の最適な減衰率を用いて前記サンプリング点における単位周波数あたりの累積減衰率を算出し、該累積減衰率を用いて前記特徴量の減衰補正を行うことによって補正特徴量を算出する特徴量補正ステップと、
     を実行させることを特徴とする超音波観測装置の作動プログラム。
PCT/JP2015/078548 2014-12-22 2015-10-07 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム WO2016103849A1 (ja)

Priority Applications (4)

Application Number Priority Date Filing Date Title
EP15872418.7A EP3238632B1 (en) 2014-12-22 2015-10-07 Ultrasound observation apparatus, method for operating ultrasound observation apparatus, and program for operating ultrasound observation apparatus
CN201580004506.2A CN105916448B (zh) 2014-12-22 2015-10-07 超声波观测装置以及超声波观测装置的工作方法
JP2016505644A JP5927367B1 (ja) 2014-12-22 2015-10-07 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
US15/207,839 US20160317126A1 (en) 2014-12-22 2016-07-12 Ultrasound observation apparatus, method for operating ultrasound observation apparatus, and computer-readable recording medium

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2014259474 2014-12-22
JP2014-259474 2014-12-22

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US15/207,839 Continuation US20160317126A1 (en) 2014-12-22 2016-07-12 Ultrasound observation apparatus, method for operating ultrasound observation apparatus, and computer-readable recording medium

Publications (1)

Publication Number Publication Date
WO2016103849A1 true WO2016103849A1 (ja) 2016-06-30

Family

ID=56149886

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2015/078548 WO2016103849A1 (ja) 2014-12-22 2015-10-07 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム

Country Status (4)

Country Link
US (1) US20160317126A1 (ja)
EP (1) EP3238632B1 (ja)
CN (1) CN105916448B (ja)
WO (1) WO2016103849A1 (ja)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2010246640A (ja) * 2009-04-13 2010-11-04 Aloka Co Ltd 超音波診断装置
WO2014054469A1 (ja) * 2012-10-01 2014-04-10 オリンパスメディカルシステムズ株式会社 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
WO2014192954A1 (ja) * 2013-05-29 2014-12-04 オリンパスメディカルシステムズ株式会社 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム

Family Cites Families (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040158154A1 (en) * 2003-02-06 2004-08-12 Siemens Medical Solutions Usa, Inc. Portable three dimensional diagnostic ultrasound imaging methods and systems
JP4585326B2 (ja) * 2005-02-08 2010-11-24 富士フイルム株式会社 超音波撮像装置及び超音波撮像方法
WO2007003058A1 (en) * 2005-07-06 2007-01-11 National Research Council Of Canada Method and system for determining material properties using ultrasonic attenuation
US8187191B2 (en) * 2009-01-08 2012-05-29 Volcano Corporation System and method for equalizing received intravascular ultrasound echo signals
US20110288410A1 (en) * 2010-02-22 2011-11-24 Speyer Gavriel A Methods and systems for diagnostic ultrasound based monitoring of high intensity focused ultrasound therapy
JPWO2012063928A1 (ja) * 2010-11-11 2014-05-12 オリンパスメディカルシステムズ株式会社 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
CN102802536B (zh) * 2010-11-11 2015-01-07 奥林巴斯医疗株式会社 超声波诊断装置、超声波诊断装置的动作方法以及超声波诊断装置的动作程序
JP5079177B2 (ja) * 2010-11-11 2012-11-21 オリンパスメディカルシステムズ株式会社 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
WO2013022454A1 (en) * 2011-08-11 2013-02-14 Delphinus Medical Technologies, Inc. Method for imaging a volume of tissue
EP2878271A4 (en) * 2012-07-27 2016-04-13 Olympus Corp ULTRASOUND OBSERVATION APPARATUS, METHOD FOR OPERATING ULTRASOUND OBSERVATION APPARATUS, AND OPERATING PROGRAM FOR ULTRASONIC OBSERVATION APPARATUS
EP2859851B1 (en) * 2013-09-24 2018-10-31 Siemens Medical Solutions USA, Inc. Shear wave estimation from analytic data
US10512450B2 (en) * 2013-09-24 2019-12-24 Siemens Medical Solutions Usa, Inc. Shear wave estimation from analytic data

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2010246640A (ja) * 2009-04-13 2010-11-04 Aloka Co Ltd 超音波診断装置
WO2014054469A1 (ja) * 2012-10-01 2014-04-10 オリンパスメディカルシステムズ株式会社 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
WO2014192954A1 (ja) * 2013-05-29 2014-12-04 オリンパスメディカルシステムズ株式会社 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム

Also Published As

Publication number Publication date
EP3238632B1 (en) 2019-08-07
CN105916448B (zh) 2019-03-01
EP3238632A4 (en) 2018-07-25
US20160317126A1 (en) 2016-11-03
CN105916448A (zh) 2016-08-31
EP3238632A1 (en) 2017-11-01

Similar Documents

Publication Publication Date Title
WO2016151951A1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP5974210B2 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
US10201329B2 (en) Ultrasound observation apparatus, method for operating ultrasound observation apparatus, and computer-readable recording medium
JP6513224B2 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP6289772B2 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP2016202567A (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP6892320B2 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
WO2016103848A1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
WO2016181869A1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP6157790B1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
WO2015198712A1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP2017113145A (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP6138402B2 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP5981072B1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP6253572B2 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP5927367B1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
WO2016103849A1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP5953457B1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
WO2016181856A1 (ja) 超音波診断装置、超音波診断装置の作動方法および超音波診断装置の作動プログラム
JP5932184B1 (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム
JP2017217313A (ja) 超音波観測装置、超音波観測装置の作動方法および超音波観測装置の作動プログラム

Legal Events

Date Code Title Description
ENP Entry into the national phase

Ref document number: 2016505644

Country of ref document: JP

Kind code of ref document: A

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

Ref document number: 15872418

Country of ref document: EP

Kind code of ref document: A1

REEP Request for entry into the european phase

Ref document number: 2015872418

Country of ref document: EP

NENP Non-entry into the national phase

Ref country code: DE