Sonar image filtering method and system based on dual-tree discrete orthogonal S transformation

Abstract

The invention discloses a sonar image filtering method and system based on dual-tree discrete orthogonal S transformation, and belongs to the field of sonar image processing, and the method comprises the following steps: obtaining a sonar image, and carrying out the DOST transformation; the frequency spectrum image obtained after DOST transformation is decomposed into a real part frequency spectrum image and an imaginary part frequency spectrum image through the orthogonality of a frequency domain; respectively acquiring optimal segmentation thresholds of the real part spectrum image and the imaginary part spectrum image by using an improved OTSU algorithm; and obtaining a filter function by using the optimal segmentation threshold and completing the filtering of the sonar image. According to the method, optimization can be carried out more flexibly and finely for different types of signal components in the frequency domain filtering process, better image denoising and information retention effects are achieved, and image information can be better processed.

Description

Technical Field

[0001] The present invention relates to the technical field of sonar image preprocessing, and more specifically, to a sonar image filtering method and system based on a dual-tree discrete orthogonal S transform. Background Art

[0002] Underwater imaging detection technologies include optical imaging and sonar imaging. The operating range of optical imaging is relatively short, generally between ten meters and dozens of meters, and it basically fails in turbid water areas. Sonar images have the advantages of long operating range and strong penetration ability, and are particularly suitable for turbid water areas. However, due to the complex and variable characteristics of the water medium and its boundaries in the aquatic channel, combined with the projection characteristics of sound waves themselves, the images collected by imaging sonars often have characteristics such as strong noise and low resolution, seriously affecting underwater detection and operations. Therefore, image filtering, as a way of image preprocessing, plays an important role in subsequent image segmentation and edge detection, etc.

[0003] Sonar image filtering is similar to traditional optical image filtering, both being planar distribution maps of energy. There are two filtering methods for sonar images. One is spatial domain filtering (adjusting the gray values of image pixels and directly processing the image), including linear filtering and non-linear filtering. The other is transform domain filtering (converting the signal to another space, indirectly filtering, and finally performing an inverse transformation to restore it), mainly involving the mutual conversion between the time domain and the frequency domain. Fourier transform, wavelet transform, and numerous time-frequency conversion methods improved based on them have many applications in the field of image processing. Among them, the S transform proposed by R.G. Stockwell in 1996 has great application prospects in signal processing. The S transform combines the advantages of the Short Time Fourier Transform (STFT) and wavelet transform, and is a multi-resolution analysis method. The S transform improves the disadvantage of the low time-frequency resolution of STFT, and can be regarded as an extension or special case of wavelet transform theory under the multi-resolution structure. The S transform inherits and develops the concept of locality in the time-frequency domain of the short-time Fourier transform and wavelet transform, has good time-frequency focusing and contains a phase factor, and is suitable for time-frequency analysis of non-stationary signals. In 1997, R.G. Stockwell also proposed a method for processing two-dimensional data using the S transform, which has strong applicability in the field of image processing. Aiming at the disadvantage that the S transform coefficient matrix is highly redundant and not convenient for practical applications, R.G. Stockwell proposed a non-redundant discrete orthogonal S transform in 2006. In 2009, Wang Y. and Orchard J. proposed the Fast Discrete Orthonormal Stockwell Transform by borrowing the idea of improving the operation efficiency in the calculation process of the Fast Fourier Transform, and applied it to the field of image compression.

[0004] Through experiments, the above methods cannot meet the filtering requirements of sonar images in the sonar image field. The filtering requirements of sonar images are: it must be able to remove speckle noise and should try to maintain the original spatial resolution, which means that the result of image filtering should not only smooth the image but also keep the local information of the image from being destroyed as much as possible. After a large number of experiments, it is proved that using a single spatial domain processing method or transform domain processing method to filter sonar images has an unsatisfactory effect. Although it can make the image smoother, it also filters out a large amount of useful information, making the image blurred and not conducive to subsequent operations such as image segmentation. The reason is that most traditional filtering methods are aimed at additive noises such as Gaussian white noise, while the noise in sonar images is mainly multiplicative speckle noise, and its components are relatively complex. Summary of the Invention

[0005] In view of this, the present invention provides a sonar image filtering method and system based on the dual-tree discrete orthogonal S-transform, which is used to at least partially solve the technical problems mentioned in the background art and improve the filtering quality of sonar images.

[0006] In order to achieve the above object, the present invention adopts the following technical solutions:

[0007] A sonar image filtering method based on the dual-tree discrete orthogonal S-transform includes the following steps:

[0008] Obtain a sonar image and perform a DOST transform;

[0009] Utilize the orthogonality in the frequency domain to decompose the spectral image obtained after the DOST transform into a real-part spectral image and an imaginary-part spectral image;

[0010] Utilize an improved OTSU algorithm to respectively obtain the optimal segmentation thresholds of the real-part spectral image and the imaginary-part spectral image;

[0011] Utilize the optimal segmentation thresholds to obtain a filter function and complete the filtering of the sonar image.

[0012] Specifically, in the step of obtaining a sonar image and performing a DOST transform in the present invention, the sonar image is obtained through the following steps:

[0013] Use a transducer to emit multiple ultrasonic waves towards the seabed detection area at a certain pitch angle and fan-shaped;

[0014] Receive the scattered echo signals generated by each ultrasonic wave;

[0015] Determine the pixel brightness of the target sonar image according to the intensity of the scattered echo signals;

[0016] Perform spatial position analysis on the scattered echo signals to obtain the position coordinates of the target sonar image;

[0017] Generate the target sonar image based on the pixel brightness and position coordinates of the target sonar image.

[0018] Specifically, performing a DOST transform on the obtained sonar image specifically includes the following steps:

[0019] S11. For the obtained sonar image I(x, y) with a size of M×N, perform a two-dimensional discrete Fourier transform to obtain the spectrum after the two-dimensional discrete Fourier transform; this process can be expressed as:

[0020]

[0021] Among them, x and y respectively represent the abscissa and ordinate in the sonar image I(x, y), H[m, n] represents the spectrum after two-dimensional discrete Fourier transform, and m and n represent the abscissa and ordinate in the frequency domain after two-dimensional discrete Fourier transform.

[0022] S12. Perform frequency-domain partitioning on the spectrum after two-dimensional discrete Fourier transform to obtain a number of frequency band blocks;

[0023] During the frequency-domain partitioning process, it is divided into orthogonal frequency band blocks according to the rules of DOST basis functions. The frequency band partitioning needs to satisfy: partitioning along the radial direction (frequency scale) and angular direction (direction), usually partitioning the frequency according to binary or logarithmic scale; each frequency band block corresponds to a specific frequency range and direction to ensure the orthogonality of the basis functions.

[0024] S13. Perform inverse two-dimensional discrete Fourier transform on each of the several frequency band blocks to obtain the time-frequency representation result of each frequency band block;

[0025] In a specific embodiment, the processes of steps S12 and S13 can be expressed as:

[0026]

[0027] Among them, the sampling frequency p in the x direction during the frequency-domain partitioning process x and the sampling frequency p in the y direction y are respectively p x = 2, …, log2(M) - 1, p y = 2, …, log2(N) - 1; ν x , ν y respectively represent the horizontal and vertical sound frequencies, and the sound frequencies are defined by the following formula:

[0028]

[0029] Among them, p represents the sampling frequency p x or p y , and is an integer greater than or equal to 0.

[0030] S14. Stitch and combine all the time-frequency representation results to obtain the DOST transform spectrum image of the sonar image, and this process is expressed as:

[0031]

[0032] Specifically, in a specific embodiment of the present invention, the spectrum image obtained after DOST transform is decomposed into a real part spectrum image and an imaginary part spectrum image by using the orthogonality in the frequency domain, including the following formula:

[0033] DOST(x′, y′) = DOST Re (x′, y′) + jDOSTIm (x′, y′)

[0034] Among them, DOST(x′, y′) represents the spectral image obtained after the DOST transformation; DOST Re The real part spectral image after the decomposition of (x′, y′); DOST Im The imaginary part spectral image after the decomposition of (x′, y′), where x′ and y′ are the horizontal and vertical coordinates of the spectral image, and j represents the imaginary part marker.

[0035] In a specific embodiment of the present invention, in the step of using the improved OTSU algorithm to respectively obtain the optimal segmentation thresholds of the real part spectral image and the imaginary part spectral image, the steps of obtaining the optimal segmentation threshold of the real part spectral image include the following steps:

[0036] Divide the real part spectral image values into N parts, and the size of each part is (A max -A min ) / N, where A max is the maximum value of the amplitude of the real part spectral image, and A min is the minimum value of the amplitude of the real part spectral image. Denote the midpoint of each part as α n , and each part is the neighborhood of the middle α n ;

[0037] Take the value of the middle α n of each part of the data as the candidate threshold α k . Based on the candidate threshold α k , divide the real part spectral image values into the class sets C1[α k+1 , α k+1 , α k+2 , … α N-1 and C2[α0, α2, α3, … α k ;

[0038] Calculate the between-class variance of all candidate thresholds according to the class sets C1 and C2, and take the candidate threshold corresponding to the maximum between-class variance as the optimal segmentation threshold of the real part spectral image.

[0039] Specifically, in a specific embodiment of the present invention, to calculate the between-class variance of the candidate threshold, the following formula is specifically adopted:

[0040]

[0041] Among them, P1(k) is the probability that the neighborhood of each element of the class set C1 occurs, and P2(k) is the probability that the neighborhood of each element of the class set C2 occurs; m1(k) and m2(k) are the average values in the sets C1 and C2 respectively, and m G is the global mean.

[0042] Specifically, in the specific embodiments of the present invention, obtaining the filter function by using the optimal segmentation threshold specifically includes the following steps:

[0043] Obtaining the filter function by using the optimal segmentation threshold specifically includes the following steps:

[0044] Denote the optimal threshold of the real part spectrum as T1 and the optimal threshold of the imaginary part spectrum as T2, and obtain the following filter function:

[0045]

[0046] where F(x′, y′) represents the filter function, DOST Re (x′, y′) represents the decomposed real part spectrum image, DOST Im (x′, y′) is the decomposed imaginary part spectrum image.

[0047] Specifically, in the specific embodiments of the present invention, filtering the sonar image specifically includes the following steps:

[0048] Multiply the spectrum image after the DOST transformation of the sonar image by convolution with the obtained filter function to obtain the filtered spectrum image; it can be expressed as:

[0049]

[0050] where represents the obtained filtered spectrum image, F(x′, y′) represents the filter function, and × represents the convolution operation.

[0051] Perform DOST reconstruction on the filtered spectrum image by using the inverse Fourier transform to obtain the filtered sonar image. DOST reconstruction is the process of changing the spectrum image back to the sonar image. The specific operation is to apply the forward 2D-FFT to the "sound" of each region to reverse the spectrum segmentation and reconstruct the spectrum of the image. The reconstruction formula is as follows:

[0052]

[0053] H[m', n', p x , p y is the Fourier spectrum after frequency shift of H[m, n].

[0054]

[0055] where

[0056] m = m - v

[0057] n = n - v

[0058]

[0059] Finally, the image is restored by inverse FT:

[0060]

[0061] On the other hand, the present invention also discloses a computer system, including a computer program, which can implement the sonar image filtering method based on the dual-tree discrete orthogonal S transform according to any one of the present invention when executed.

[0062] Through the above technical solutions, compared with the prior art, the present invention discloses a sonar image filtering method and system based on the dual-tree discrete orthogonal S transform, which has the following beneficial effects:

[0063] The present invention decomposes the frequency spectrum after discrete orthogonal S transform into a binary tree structure, one part is the real part and the other part is the imaginary part. Filter the two parts separately and then recombine them into a new frequency spectrum. And the method of dividing the threshold is improved, and the OTSU algorithm based on the gray value probability applicable to binarization processing is improved to make it applicable to frequency domain filtering. After quality evaluation, the method proposed in this paper has obvious improvement in filtering quality compared with the commonly used transform domain filtering methods. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required to be used in the description of the embodiments or the prior art. Obviously, the drawings in the following description are only the embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained according to the provided drawings without creative efforts.

[0065] Figure 1 It is a schematic diagram of the overall flow of the method provided by the present invention.

[0066] Figure 2 It is a schematic diagram of the division of the real part spectrum of DOST when N = 16 provided by the embodiment of the present invention.

[0067] Figure 3 It is the original image of the sonar image and the spectrum image after DOST transform provided by the embodiment of the present invention.

[0068] Figure 4 It is the decomposition and reconstruction of two-stage DTDOST provided by the embodiment of the present invention.

[0069] Figure 5 It is the original image, the noise image and their spectral images provided by the embodiment of the present invention.

[0070] Figure 6 It is the original image and the gray histogram provided by the embodiment of the present invention.

[0071] Figure 7 The sonar images and filtered images with uneven histograms provided by the embodiments of the present invention.

[0072] Figure 8 Another group of original images and grayscale histograms provided by the embodiments of the present invention.

[0073] Figure 9 The sonar images and filtered images with uniform histograms provided by the embodiments of the present invention.

[0074] Figure 10 The third group of original images and grayscale histograms provided by the embodiments of the present invention.

[0075] Figure 11 The sonar images and filtered images with complex image information provided by the embodiments of the present invention. Detailed implementation manners

[0076] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0077] Embodiment 1

[0078] Embodiment 1 discloses a sonar image filtering method based on the dual-tree discrete orthogonal S transform, as Figure 1 shown, including the following steps:

[0079] Obtain a sonar image and perform a DOST transform;

[0080] Decompose the spectral image obtained after the DOST transform into a real-part spectral image and an imaginary-part spectral image by using the orthogonality in the frequency domain;

[0081] Use the improved OTSU algorithm to obtain the optimal segmentation thresholds of the real-part spectral image and the imaginary-part spectral image respectively;

[0082] Obtain a filter function by using the optimal segmentation threshold and complete the filtering of the sonar image.

[0083] Specifically, in the step of obtaining a sonar image and performing a DOST transform in the present invention, the sonar image is obtained through the following steps:

[0084] Use a transducer to emit multiple ultrasonic waves to the seabed detection area at a certain pitch angle and fan-shaped;

[0085] Receive the scattered echo signals generated by each ultrasonic wave;

[0086] Determine the pixel brightness of the target sonar image according to the intensity of the scattered echo signal;

[0087] Perform spatial position analysis on the scattered echo signal to obtain the position coordinates of the target sonar image;

[0088] Generate a target sonar image based on the pixel brightness and position coordinates of the target sonar image.

[0089] Specifically, perform a DOST transform on the acquired sonar image, which specifically includes the following steps:

[0090] S11. For the acquired sonar image I(x, y) with a size of M×N, perform a two-dimensional discrete Fourier transform to obtain the spectrum after the two-dimensional discrete Fourier transform; this process can be expressed as:

[0091]

[0092] Among them, x and y respectively represent the abscissa and ordinate in the sonar image I(x, y), H[m, n] represents the spectrum after the two-dimensional discrete Fourier transform, and m and n represent the abscissa and ordinate in the frequency domain after the two-dimensional discrete Fourier transform.

[0093] S12. Partition the spectrum after the two-dimensional discrete Fourier transform in the frequency domain to obtain several frequency band blocks;

[0094] During the frequency domain partitioning process, it is divided into orthogonal frequency band blocks according to the rules of the DOST basis function. The frequency band partitioning needs to satisfy: partition along the radial direction (frequency scale) and angular direction (direction), usually divide the frequency by binary or logarithmic scale; each frequency band block corresponds to a specific frequency range and direction to ensure the orthogonality of the basis function.

[0095] S13. Perform an inverse two-dimensional discrete Fourier transform on each of the several frequency band blocks to obtain the time-frequency representation result of each frequency band block;

[0096] In a specific embodiment, the processes of steps S12 and S13 can be expressed as:

[0097]

[0098] Among them, the sampling frequency p in the x direction during the frequency domain partitioning process x and the sampling frequency p in the y direction y are respectively p x = 2,..., log2(M) - 1, p y = 2,..., log2(N) - 1; ν x , ν y respectively represent the horizontal and vertical sound frequencies, and the sound frequencies are defined by the following formula:

[0099] When p = 0, v x = 0, v y = 0;

[0100] When p = 1, v x = 1, v y = 1;

[0101] When p > 1,

[0102] where p represents the sampling frequency p x or p y , and is an integer greater than or equal to 0.

[0103] S14. Combine all the time-frequency representation results to obtain the DOST transform frequency spectrum image of the sonar image, and this process is expressed as:

[0104]

[0105] Specifically, in the specific embodiment of the present invention, the frequency spectrum image obtained after the DOST transform is decomposed into a real part frequency spectrum image and an imaginary part frequency spectrum image by using the orthogonality in the frequency domain, including the following formula:

[0106] DOST(x′, y′) = DOST Re (x′, y′) + jDOST Im (x′, y′)

[0107] where DOST(x′, y′) represents the frequency spectrum image obtained after the DOST transform; DOST Re (x′, y′) is the real part frequency spectrum image after decomposition; DOST Im (x′, y′) is the imaginary part frequency spectrum image after decomposition, x′, y′ are the horizontal and vertical coordinates of the frequency spectrum image, and j represents the imaginary part marker.

[0108] In the specific embodiment of the present invention, in the step of using the improved OTSU algorithm to respectively obtain the optimal segmentation thresholds of the real part frequency spectrum image and the imaginary part frequency spectrum image, the steps of obtaining the optimal segmentation threshold of the real part frequency spectrum image include the following steps:

[0109] Divide the real part frequency spectrum image values into N parts, and the size of each part is (A max - A min ) / N, where A max is the maximum value of the real part frequency spectrum image amplitude, A min is the minimum value of the real part frequency spectrum image amplitude, and take the midpoint of each part as α n , and each part is the neighborhood of the middle α n ;

[0110] For the middle α of each part of the datan The value is used as the candidate threshold α k , based on the candidate threshold α k , the real part spectral image values are divided into class sets C1[α k+1 , α k+1 , α k+2 , … α N-1 and C2[α0, α2, α3, … α k ;

[0111] Calculate the between-class variance of all candidate thresholds according to the class sets C1 and C2, and take the candidate threshold corresponding to the maximum between-class variance as the optimal segmentation threshold of the real part spectral image.

[0112] Specifically, to calculate the between-class variance of the candidate threshold, the following formula is specifically adopted:

[0113]

[0114] Among them, P1(k) is the probability of the neighborhood occurrence of each element in the class set C1, P1(k)P2(k) is the probability of the neighborhood occurrence of each element in the class set C2; m1(k) and m2(k) are the average values in the sets C1 and C2 respectively, and m G is the global mean.

[0115] The traditional OTSU algorithm can only take integers between 0 and 255 as the threshold, but there are decimals and negative numbers in the real part of the DOST spectrum, so the traditional OTSU cannot be used for calculation. Therefore, in this paper, the spectral image of the real part is divided into N parts, and the middle value α n of each part is taken for calculation, from which the probability can be calculated. When m(k) is the largest, the value of α n is the threshold T to be taken.

[0116] Specifically, in the specific embodiment of the present invention, to obtain the filter function using the optimal segmentation threshold, the following steps are specifically included:

[0117] To obtain the filter function using the optimal segmentation threshold, the following steps are specifically included:

[0118] Record the obtained optimal threshold of the real part spectrum as T1 and the optimal threshold of the imaginary part spectrum as T2, and obtain the following filter function:

[0119]

[0120] Among them, F(x′, y′) represents the filter function, and DOST Re (x′, y′) represents the decomposed real part spectral image, and DOST Im (x′, y′) is the decomposed imaginary part spectral image.

[0121] Specifically, in the specific embodiments of the present invention, the filtering of the sonar image is completed, specifically including the following steps:

[0122] The filtered frequency spectrum image is obtained by convolving and multiplying the frequency spectrum image after the DOST transformation of the sonar image with the obtained filtering function; it can be expressed as:

[0123]

[0124] Wherein, represents the obtained filtered frequency spectrum image, F(x′,y′) represents the filter function, and × represents the convolution operation.

[0125] The filtered frequency spectrum image is subjected to DOST reconstruction using the inverse Fourier transform to obtain the filtered sonar image. DOST reconstruction is the process of changing the frequency spectrum image back into the sonar image. The specific operation is to apply the forward 2D-FFT to the "sound" of each region to reverse the frequency spectrum segmentation and reconstruct the frequency spectrum of the image. The reconstruction formula is as follows:

[0126]

[0127] H[m',n',p x ,p y is the Fourier spectrum after frequency shift of H[m,n].

[0128]

[0129] Wherein,

[0130] m = m - v

[0131] n = n - v

[0132]

[0133]

[0134] Finally, the image is restored through the inverse Fourier transform FT:

[0135]

[0136] The technical principle in the present invention is further described below.

[0137] The S transform is a time-frequency transform method that can provide the frequency and phase information of a signal. Given a one-dimensional continuous signal h(t), the definition of its corresponding S transform is

[0138]

[0139] Where f is the frequency of the signal, and t and τ are the time and space variables of the signal. The basic wavelet in the transform is defined as

[0140]

[0141] It can be seen from the definition formula of the one-dimensional S transform that the Gaussian window in the S transform is inversely proportional to the frequency, and higher time and frequency resolutions are obtained through the narrower time window at high frequencies and the wider time window at low frequencies. This enables the S transform to describe the signal characteristics more precisely.

[0142] The discrete orthogonal s transform solves the redundancy of the original S transform. It is an algorithm that uses less time-domain decomposition and correspondingly fewer frequency bands according to the sampling theory. DOST has a longer period at low frequencies, corresponding to a lower sampling interval; similarly, it has a higher sampling interval at high frequencies. DOST divides the time-frequency domain into N regions by constructing N orthogonal basis vectors. The k-th basis vector can be defined as:

[0143]

[0144] where k = 0,..., N - 1, υ represents the center of the frequency band, β represents the bandwidth, and τ represents the time. After summation:

[0145]

[0146] where The selection of the parameters υ, β, and τ affects the orthogonality of the basis function. Stockwell gave the variable p, which represents the number of octaves of the signal's multiple-frequency sampling. The basis vector parameters υ, β, and τ of the positive frequency of the DOST algorithm are:

[0147] When p = 0, v = 0; β = 1; τ = 0;

[0148] When p = 1, v = 1; β = 1; τ = 0;

[0149] When p = 2, 3,..., log2N - 1,

[0150] Based on the above orthogonal basis vectors and the relationship between the S transform and the Fourier transform, the formula of DOST can be deduced as follows:

[0151]

[0152] where h(t) is the initial image signal. The inverse transform of the S transform is the inverse Fourier transform, that is

[0153]

[0154] The embodiments of the present invention give the relevant proof of the inverse transform of the S transform. From the convolution integral property

[0155] x(t)*h(t) = y(t)

[0156] y (-1) (t) = x (-1) (t)*h(t) = x(t)*h (-1) (t)

[0157] It can be obtained that

[0158]

[0159] By the Gaussian integral formula Therefore

[0160] Thus

[0161]

[0162] Through the above proof, the inverse transform of the S transform can be transformed into the original image by the inverse Fourier transform.

[0163] For sonar images, the brightness value of each pixel in the sonar image can be extracted, forming a two-dimensional matrix, and expanded to obtain one-dimensional data as the one-dimensional signal of the DOST transform.

[0164] The frequency spectrum image after the DOST transform satisfies the orthogonality in the frequency domain. The frequency spectrum image is decomposed into two parts, the real part and the imaginary part, which are respectively denoted as DOST Re and DOST Im . Then the frequency spectrum image after the DOST transform is

[0165] DOST(τ, f) = DOST Re (τ, f) + jDOST Im (τ, f)

[0166] Expand the DOST formula using Euler's formula. Then the real part DOST Re The specific formula is

[0167]

[0168] where k = (0, 1, 2…N - 1). Taking N = 16 as an example, such as Figure 2Shown is the order of converting two-dimensional DOST coefficients into one-dimensional vectors and the division of the time-frequency domain. The bold blue part in the middle is the DC component, that is, the part where f = 0. In the positive frequency part, it is divided in sequence according to the basis vector parameters of the positive frequency of the DOST algorithm. The selected values of the parameters υ, β, τ determine the frequency band center, bandwidth, and time of the image division, which are the main factors for the spectrogram to present multiple rectangular blocks. And the selected values of the parameters υ, β, τ do not interfere with the results of the dual-tree decomposition. The operation of the imaginary part is basically similar to that of the real part, and the formula for the imaginary part is

[0169]

[0170] Since DOST is an orthogonal transform and the basis functions are orthogonal bases, the real part and the imaginary part after decomposition are independent of each other. To a certain extent, it solves the problem of information redundancy of the S transform. And the spectral information existing in the real part and the imaginary part is different, and filtering them separately helps in the analysis of image information.

[0171] The present invention uses an improved OTSU algorithm to obtain the optimal threshold for each part and performs threshold division on the two parts separately. The principle is as follows. Taking the real part as an example, the imaginary part has commonality with the real part. Let the maximum value of the amplitude of the real part spectral image be A max , and the minimum value be A min . The total number of frequency points is N. The real part image values are divided into N parts, and the size of each part is (A max - A min ) / N. Take the midpoint of each part and denote it as α n , that is, each part is the neighborhood of α n .

[0172] Assume that a threshold T has been selected. C1 is a set of values of [α k+1 , α k+1 , α k+2 , … α N-1 , and C2 is a set of values of [α0, α2, α3, … α k . The OTSU algorithm is to select the optimal threshold. In a sense, it selects the threshold T such that its maximum between-class variance is

[0173]

[0174] In the formula, P1(k) is the probability that the value belongs to the neighborhood of each element in the set C1, and P2(k) is the probability that the value belongs to the neighborhood of each element in the set C2

[0175]

[0176] m1(k) and m2(k) are the average values in the sets C1 and C2 respectively. m Gis the global mean.

[0177]

[0178] Expand the expression of, the between-class variance can be written as

[0179]

[0180] where m(k) is the average value from α0 to α k The larger the variance, the closer it is to the threshold for segmenting the image. For sonar images, the image gray value is close to the noise value. The ordinary threshold division method will misjudge some parts of the image and regard some parts of the image as noise and remove them. The method of selecting the threshold through variance combines the image characteristics and avoids the errors caused by randomness.

[0181] After obtaining the optimal thresholds of the real part and the imaginary part, the filter function is as follows

[0182]

[0183] Finally, the filtered sonar image is obtained through convolution and quasi-transformation.

[0184] Embodiment 2

[0185] Embodiment 2 discloses the practical application of the method of the present invention, including the following steps:

[0186] 1) Acquisition of sonar images

[0187] In this embodiment, the sonar image is first generated by the transducer emitting ultrasonic waves obliquely downward at a certain pitch angle and fan-shaped towards the detection area on the seabed. When the emitted ultrasonic waves encounter small seabed targets or underwater media, scattering will occur, and all the scattered echo signals return to the transducer along the original path to form the image. The echo received by the transducer is often a group of pulse trains with different amplitudes and time lengths. Among them, the echo generated by the seabed area relatively close to the receiving array is collected first, and the echo generated by the seabed area farther from the receiving array is collected later. The echo of hard and convex substances on the seabed is stronger, the echo of soft and concave substances is weaker, and there is no echo in other areas blocked by objects. The closer to the sonar system, the stronger the echo. These echoes are combined into a position array through spatial position analysis to represent the position of the emission point. The data in each row is composed of the data of several beams in a certain fan-shaped surface at the same distance, and each column represents the corresponding different azimuths.

[0188] Figure 3(a) is a range-azimuth two-dimensional image with a size of 330 (rows) × 224 (columns). In the figure, the larger bright area corresponds to the area with stronger echoes and is called the target bright area. The darker part around the target bright area is the target dark area, which is caused by small underwater targets or acoustic shadows of various underwater objects or the inability of sound waves to penetrate small targets resulting in no echoes. Except for the target bright area and the target dark area, the remaining part is the seabed reverberation area. The reverberation signal formed in this part is a non-stationary random signal that decays over time according to a certain law. Through the study of its phase and amplitude, its statistical characteristics also follow a certain law of distribution. Assuming that various scattering of the sonar system is evenly distributed, the sound wave emission and reception are at the same position, and the reverberation noise is mainly caused by the scattering of various substances on the seabed, it can be obtained that the statistical model of speckle noise in the sonar image is similar to the statistical models of common ultrasonic images and SAR images in medicine, that is

[0189] h(i,j) = v(i,j)g(i,j)

[0190] where i and j are the rows and columns of the sonar image respectively, v(i,j) represents the noiseless sonar image, and g(i,j) represents the speckle noise of the image. The amplitude of the noise g(i,j) follows a Rayleigh distribution. Since the multiplicative speckle noise is the main reason for the degradation of the sonar image quality, removing the speckle noise in the sonar image is crucial for the analysis and discrimination of the sonar image.

[0191] 2) DTDOST filtering process

[0192] The DTDOST filtering method adopted in this embodiment uses a two-channel filter bank with a binary tree structure for signal decomposition and reconstruction. The first tree generates the real part, and the second tree generates the imaginary part. The low-pass filters of the real and imaginary part trees are reasonably designed. The sampling frequencies of the filters of the two trees are the same and independent of each other. In the decomposition process, the improved OTSU algorithm is used for threshold decomposition each time. After multiple decompositions, the filtered real and imaginary parts can be obtained. The two are recombined to obtain the amplitude matrix at the corresponding frequency, and then the inverse transform is used to restore it to the filtered image. Figure 4 It is the decomposition and reconstruction process of 2-level DTDOST.

[0193] Since the sonar image noise is speckle noise, and speckle noise is multiplicative noise. The noise processed by the S transform can only be additive noise, so a logarithmic transform is needed to convert the multiplicative noise into additive noise. And to satisfy the logarithmic transform, 1 is added to the pixel points with a gray value of 0 during the logarithmic transform. In the picture, R1 and I1 are the real and imaginary parts respectively after the spectral image is decomposed. The improved OTSU first-level decomposition is performed on both of them. After decomposition, R21 and R22 are the low-frequency and high-frequency parts of the real part first-level decomposition respectively, and I21 and I22 are the high-frequency and low-frequency parts of the imaginary part first-level decomposition respectively. The low-frequency part on the spectral image is the target area of the original image, so a low-pass filter is designed to remove the noise. After secondary decomposition, the low-frequency parts R31 and I31 of the real and imaginary parts are obtained. The meaning of the elements in the spectral image matrix is the amplitude at the corresponding frequency. The image reconstructed from R31 and I31 The elements in it are the 1 / 2 power of the sum of the squares of each corresponding element in R31 and I31.

[0194] Due to reasons such as the acquisition and imaging factors of different devices, the gray information of the same tissue in the image is inconsistent. Image normalization is an image conversion method that reduces or even eliminates the gray inconsistency in the image while retaining the diagnostically valuable gray differences. Image normalization is more conducive to computer automatic analysis and processing. Therefore, image normalization is to convert the original image to be processed into the corresponding standard form through a series of transformations. The normalization of the gray value of the filtered image can be defined as

[0195]

[0196] where i, j = (1, 2, 3…N - 1).

[0197] To verify the feasibility of the invented filtering method, first, a classic optical image is used as the target for filtering. The optical image has similarities with the sonar image, but compared with the sonar image, the optical image has a larger contrast, a higher resolution, and weaker interference.

[0198] Such as Figure 5 (a), the experimental object is an optical image with a size of 256×256 pixels. (b) is its spectral image after the DOST transform. The bright area in the figure is the low-frequency area of the image, the dark area is the high-frequency area of the image, and the middle cross is the area with a frequency of 0. Multiple rectangular blocks appear on the spectrogram, which is the result of the influence of the DOST orthogonal basis. Its meaning is that when processing an N×N two-dimensional image, the traditional S transform will generate N 4 coefficients, while the DOST transform generates a total of N 2For a coefficient, DOST significantly reduces the number of coefficients generated after transformation, reducing the computational complexity. (c) is the grayscale image of the original image. (d) is the salt-and-pepper noise image with a noise density of 0.2 added to the grayscale image. The values of salt-and-pepper noise are 0 or 255, appearing as black and white speckles on the image. (e) is the spectrogram of the noisy image. Since the noise is high-frequency and has a long bandwidth (β is large), the dark area of the spectrogram becomes larger and is concentrated at the corners. (f) is the image after filtering by the method in this paper. Comparing the noisy image with the original image, the noise in the filtered image has been reduced. This also shows the feasibility of the method in this paper for filtering. In addition, further experiments on sonar images need to be carried out to test whether the method in this paper can effectively filter sonar images.

[0199] Experimental Results and Quality Assessment

[0200] In this paper, filtering experiments were carried out on three real sonar images. Comparisons were also made with classical filters, including classical low-pass filters such as DOST filter, exponential filter, Gaussian filter, and Butterworth filter. Since the resolution of the human eye is limited, and there are biases in subjective human observation, and human eye observation often depends on the experience and subjective judgment of the evaluator, several commonly used performance indicators were also selected in this paper to judge the quality of the filtering effect. These include peak signal-to-noise ratio (PSNR), structural similarity (SSIM), equivalent number of looks (ENL), and edge preservation index (EPI).

[0201] In this paper, MSE and PSNR are used to evaluate the distortion degree of image filtering. MSE represents the average of the differences in pixel values at each position of two images. The larger the value, the lower the similarity between the two pictures. PSNR is one of the most commonly used parameters to measure the filtering effect. It improves MSE. The larger the value, the larger the proportion of useful information, and the better the filtering effect.

[0202]

[0203] Among them

[0204]

[0205] where k i,j is the grayscale value of the filtered image, f i,j is the grayscale value of the original image, and L MAX is the maximum grayscale value.

[0206] Another tool is structural similarity (SSIM). SSIM is an index to measure the similarity between two pictures. Like PSNR, SSIM is also often used to evaluate the distortion degree of images. The input of SSIM is two images. Assuming the two images we input are x and y respectively, SSIM can be defined as

[0207]

[0208] Among them, μ x and μ y represent the average values of x and y respectively, and σ x and σ y represent the standard deviations of x and y respectively. σ xy represents the covariance of x and y. And C1, C2, and C3 are constants respectively to avoid system errors caused by a denominator of 0. SSIM is a number between 0 and 1. The larger it is, the smaller the difference between the output image and the distortion-free image, that is, the better the image quality. When the two images are identical, SSIM = 1.

[0209] Equivalent number of looks (ENL), which is an index to measure the smoothness of a homogeneous region. ENL is a standard parameter widely used to evaluate the filtering results of homogeneous regions. ENL is defined as

[0210]

[0211] where μ z and σ z are the estimated average value and standard deviation of the original or filtered sonar image. The larger the ENL, the stronger the noise suppression ability.

[0212] While removing noise from the image by image filtering, it is also necessary to ensure that the image edges are not lost. Therefore, the edge preservation index (EPI) is selected in this paper to measure the edge preservation ability of the image. EPI is defined as

[0213]

[0214] where, where I f and I o represent the filtered sonar image and the original sonar image respectively, and i and j are the coordinates of the pixels. The closer the EPI value is to 1, the stronger the edge preservation ability.

[0215] Experiment 1 - Sonar image with non-uniform gray histogram: In this experiment, the tested image is a sonar image with a non-uniform gray histogram of size 224×330 pixels. Figure 6 The gray histogram shown has the largest number of gray values in the range of 80 - 120, with a relatively concentrated gray value range and uneven distribution.

[0216] Figure 7 are the sonar image with non-uniform histogram and the filtered images. Among them, (a) is the original sonar image. (b) is the image filtered by the DOST filter. (c) is the image filtered by the exponential filter. (d) is the image filtered by the Gaussian filter. (e) is the image filtered by the Butterworth filter. (f) is the image filtered by the DTDOST filter.

[0217] Table 1 gives the quantitative measures of the distortion degree of sonar images with uneven gray histograms after filtering by the DTDOST filtering method and four other filtering methods. Among them, the filtering method in this paper performs the best in PSNR, with the largest value and the smallest distortion degree. The DOST filtering method differs from the method in this paper by 1%. Compared with the other three methods, the performance of the filtering method in this paper is improved by about 2%, and the distortion degrees are not much different. Moreover, the SSIM value of the filtering method in this paper is closest to 1, and the distortion degree is the smallest. Among the other four methods, the DOST filtering method performs the best, with an improvement of about 4%. In SSIM, the filtering method in this paper and the DOST filtering method have relatively obvious advantages.

[0218] Table 1 Quantitative measures of the distortion degree of different filtering methods for sonar images with uneven gray histograms

[0219]

[0220] Experiment 2 - Sonar images with uniform gray histograms: The images tested in this experiment are sonar images with uniform gray histograms of size 256×256 pixels. Since the range of gray values of gray images is small, Figure 8 the number of pixel points with gray values below 120 in the displayed gray histograms is very close, so this sonar image can be called relatively uniform.

[0221] Figure 9 are sonar images with uniform histograms and filtered images. (a) Original sonar image. (b) Image filtered using the DOST filter. (c) Image filtered using the exponential filter (d) Image filtered using the Gaussian filter. (e) Image filtered using the Butterworth filter. (f) Image filtered using the DTDOST filter.

[0222] Table 2 gives the quantitative measures of the distortion degree and smoothness of sonar images with uniform gray histograms after filtering by the DTDOST filtering method and four other filtering methods. Among them, the filtering method in this paper performs the best in PSNR, with the largest value and the smallest distortion degree. The DOST filtering method is better than the other three methods, with a performance improvement of about 1%, and the distortion degrees are not much different. Moreover, the SSIM value of the filtering method in this paper is closest to 1, and the distortion degree is the smallest. Among the other four methods, the DOST filtering method performs the best, with an improvement of about 4%. In SSIM, the filtering method in this paper and the DOST filtering method have relatively obvious advantages. Comparing the ENL values of uniform images, the exponential filtering method performs the best and is not much different from the Gaussian filtering and Butterworth filtering. Compared with the method in this paper, the ENL performance of other methods is better than that of the method in this paper. Because the method in this paper not only effectively removes noise during the filtering process but also retains the original image without loss, there is a slight deficiency in smoothness compared to other methods.

[0223] Table II Quantitative measurement of the distortion degree of sonar images with uniform gray - level histograms using different filtering methods

[0224]

[0225]

[0226] Experiment Three - Sonar Images with Complex Image Information: In this experiment, the sonar images with complex image information tested have a size of 172×164 pixels. Figure 10 (b) shows that this sonar image is also an image with non - uniform gray - level values. However, compared with the images used in Experiment One and Experiment Two, the images used in Experiment Three have more complex image information, the shapes of the target bright areas and target dark areas are irregular, and the noise distribution in the seabed reverberation area is more chaotic.

[0227] Figure 11 Sonar images with complex image information and filtered images. (a) Original sonar image. (b) Image filtered using the DOST filter. (c) Image filtered using the exponential filter. (d) Image filtered using the Gaussian filter. (e) Image filtered using the Butterworth filter. (f) Image filtered using the DTDOST filter.

[0228] Table III gives the quantitative measurements of the distortion degree, smoothness degree, and edge - preserving ability of the DTDOST filtering method and four other filtering methods for sonar images with complex image information after filtering. For PSNR, SSIM, and ENL, the quantitative measurements of the images are consistent with the conclusions obtained in Experiment One and Experiment Two. For EPI, the numerical value of the filtering method in this paper is closest to 1, and the edge - preserving ability is the strongest. The EPI of the DOST filtering method is 15% lower than that of this paper's method and about 30% higher than the other three methods.

[0229] Table III Quantitative measurement of the distortion degree of sonar images with complex image information using different filtering methods

[0230]

[0231] After experiments, in the method of this paper, the distortion degree, smoothness degree, and detail - preserving ability are also contradictory like most methods. However, the method of this paper performs better than several common low - pass filters. Compared with the original DOST filtering method, the method of this paper only has a slight improvement in several evaluation indicators.

[0232] The sonar image filtering method of DTDOST proposed by the present invention is a transform domain filtering method. It is improved on the basis of the original DOST method. By leveraging the characteristics of the original orthogonal basis, it performs binary tree decomposition on it, and processes the real part and the imaginary part separately, enabling more flexible and refined optimization for different types of signal components during the frequency domain filtering process, achieving better image denoising and information retention effects, and being able to better process image information. And the threshold division method is improved, enabling the OTSU threshold division method originally applicable to grayscale image binarization to be used in the frequency domain. It avoids the error of artificially dividing the threshold, improves the stability and robustness of the processing results, and performs particularly well in dealing with complex noise and detail retention.

[0233] Experiments have proved that the method of the present invention has good effects in both noise suppression and detail preservation, weakens the phenomenon of image distortion during filtering, and also retains the image edges. Generally, the present invention can obtain a good filtering effect on sonar images.

[0234] In this specification, the various embodiments are described in a progressive manner. The key point of each embodiment is to illustrate the differences from other embodiments. The same or similar parts among the various embodiments can be referred to each other. For the devices disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the description is relatively simple, and the relevant parts can be referred to the description of the method part.

[0235] The above description of the disclosed embodiments enables those skilled in the art to implement or use the present invention. Various modifications to these embodiments will be obvious to those skilled in the art. The general principles defined herein can be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to these embodiments shown herein, but rather will conform to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A sonar image filtering method based on dual-tree discrete orthogonal S transform, characterized in that: The following steps are involved: Acquire sonar images and perform DOST transformation; The spectrum image obtained after DOST transformation is decomposed into real spectrum image and imaginary spectrum image by using the orthogonality of frequency domain; The improved OTSU algorithm is used to obtain the optimal segmentation thresholds of the real spectrum image and the imaginary spectrum image respectively; The optimal segmentation threshold is used to obtain a filter function and complete the filtering of the sonar image.

2. The sonar image filtering method based on dual-tree discrete orthogonal S transform according to claim 1 is characterized in that: In the step of acquiring a sonar image and performing a DOST transformation, the sonar image is acquired by the following steps: The transducer is used to send multiple ultrasonic beams to the seabed detection area at a certain pitch angle and fan direction; Receive the scattered echo signal generated by each ultrasonic beam; Determining the pixel brightness of the target sonar image according to the intensity of the scattered echo signal; Performing spatial position analysis on the scattered echo signal to obtain the position coordinates of the target sonar image; Generate a target sonar image based on the pixel brightness and position coordinates of the target sonar image.

3. The sonar image filtering method based on dual-tree discrete orthogonal S transform according to claim 1 is characterized in that: The DOST transformation of the acquired sonar image includes the following steps: S11. For the acquired sonar image I(x, y) of size M×N, perform a two-dimensional discrete Fourier transform to obtain a spectrum after the two-dimensional discrete Fourier transform; S12. partitioning the spectrum after the two-dimensional discrete Fourier transform in the frequency domain to obtain a number of frequency band blocks; S13. Perform an inverse two-dimensional discrete Fourier transform on each of the several frequency band blocks to obtain a time-frequency representation result of each frequency band block; S14. All the time-frequency representation results are combined and spliced ​​to obtain the DOST transformed spectrum image of the sonar image.

4. The sonar image filtering method based on dual-tree discrete orthogonal S transform according to claim 1 is characterized in that: The orthogonality of the frequency domain is used to decompose the spectrum image obtained after the DOST transform into a real spectrum image and an imaginary spectrum image, which specifically includes the following formulas: DOST(x′,y′)=DOST Re (x′,y′)+jDOST Im (x′,y′) Where DOST(x′,y′) represents the spectrum image obtained after DOST transformation; DOST Re (x′, y′) Real part spectrum image after decomposition; DOST Im (x′, y′) is the imaginary part spectrum image after decomposition, x′, y′ are the horizontal and vertical coordinates of the spectrum image, and j represents the imaginary part label.

5. The sonar image filtering method based on dual-tree discrete orthogonal S transform according to claim 1 is characterized in that: In the step of respectively obtaining the optimal segmentation threshold of the real spectrum image and the imaginary spectrum image by using the improved OTSU algorithm, obtaining the optimal segmentation threshold of the real spectrum image specifically includes: Divide the real spectrum image value into N parts, the size of each part is (A max -A min ) / N, where A max is the maximum amplitude of the real spectrum image, A min is the minimum amplitude of the real spectrum image, and the midpoint of each part is recorded as α n , each of which is the middle α n Neighborhood of; The middle α of each data n The value of k , based on the candidate threshold α k , divide the real spectrum image numerically into the category set C1[α k+1 ,α k+1 ,α k+2 ,…α N-1 ] and C2[α0,α2,α3,…α k ]; The inter-class variance of all candidate thresholds is calculated according to the category sets C1 and C2, and the candidate threshold corresponding to the maximum inter-class variance is taken as the optimal segmentation threshold of the real spectrum image.

6. The sonar image filtering method based on dual-tree discrete orthogonal S transform according to claim 5 is characterized in that: Calculate the inter-class variance of the candidate threshold using the following formula: Among them, P1(k) is the probability of occurrence of the neighborhood of each element of the category set C1, and P2(k) is the probability of occurrence of the neighborhood of each element of the category set C2; m1(k) and m2(k) are the average values ​​in the sets C1 and C2 respectively, and m G is the global mean.

7. The sonar image filtering method based on dual-tree discrete orthogonal S transform according to claim 1 is characterized in that: The optimal segmentation threshold is used to obtain a filter function, which specifically includes the following steps: The optimal threshold of the real spectrum is recorded as T1, and the optimal threshold of the imaginary spectrum is recorded as T2, and the following filter function is obtained: Among them, F(x′,y′) represents the filter function, DOST Re (x′, y′) represents the real part spectrum image after decomposition, DOST Im The imaginary spectrum image after (x′,y′) decomposition.

8. The sonar image filtering method based on dual-tree discrete orthogonal S transform according to claim 6 is characterized in that: The filtering of the sonar image is completed, which specifically includes the following steps: The spectrum image after DOST transformation of the sonar image is convolved with the obtained filter function to obtain a filtered spectrum image; The filtered spectrum image is reconstructed by inverse Fourier transform to obtain the filtered sonar image.

9. A computer system, characterized in that: It includes a computer program, which, when executed, can implement the sonar image filtering method based on dual-tree discrete orthogonal S transform as described in any one of claims 1 to 8.

Images