An adaptive high-frequency signal filtering and enhancement method for seabed pipeline SAS imaging
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- HARBIN INST OF TECH (SHENYANG) INTELLIGENT IND TECH CO LTD
- Filing Date
- 2026-05-26
- Publication Date
- 2026-07-21
Smart Images
Figure CN122260330B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of signal processing, specifically to an adaptive high-frequency signal filtering and enhancement method for SAS imaging of subsea pipelines. Background Technology
[0002] In synthetic aperture sonar imaging systems, high-frequency sound waves are emitted and seabed echo signals are received. After processing such as acoustic propagation correction, pulse compression, and azimuth-matched filtering, high-resolution images of seabed targets are reconstructed.
[0003] However, existing synthetic aperture sonar imaging systems still face several technical challenges in processing high-frequency echo signals. Firstly, due to the combined effects of seabed inhomogeneity and distance migration during high-frequency sound propagation, the statistical characteristics corresponding to different regional attributes in the echo signal exhibit complex coupling relationships. Existing systems, relying solely on local statistics, struggle to accurately identify the region to which each imaging unit belongs. Secondly, due to the limited accuracy of region type identification, subsequent filtering processes cannot adaptively match the corresponding processing mechanism based on the actual regional attributes of the imaging unit. This makes it difficult for existing systems to simultaneously achieve effective suppression of speckle noise and accurate preservation of key structural features such as pipe edges in complex seabed environments. Summary of the Invention
[0004] To achieve the above objectives, the present invention provides the following technical solution: an adaptive high-frequency signal filtering and enhancement method for SAS imaging of subsea pipelines, comprising: S1. Perform acoustic propagation compensation on the original high-frequency echo signal. Acoustic propagation compensation includes calculating the ray bending correction amount based on the sound velocity profile data in the marine environment and performing frequency and distance-related acoustic absorption compensation based on the sound absorption coefficient of the seawater medium to obtain the echo signal after acoustic propagation correction. Perform range-direction pulse compression processing on the echo signal after acoustic propagation correction to obtain the first distance data, and perform homomorphic logarithmic transformation on it to obtain logarithmic domain distance data. S2. Traverse each resolution cell in the logarithmic distance data. For the current resolution cell, calculate the local gray mean and local gray variance within its neighborhood window, and construct the speckle noise intensity factor of the resolution cell based on the ratio of the local variance to the local mean. At the same time, extract the distance migration trajectory features in the resolution cell, and perform weighted fusion of the speckle noise intensity factor and the distance migration trajectory feature value of each resolution cell to construct a joint discriminant factor. S3. Adaptively determine a global noise threshold set based on the joint discriminant factor of all resolution units; S4. Compare the joint discriminant factor of each resolution unit with the global noise threshold set, and classify the resolution unit into one of the following categories based on the comparison results: noise-dominated region, edge texture region, or flat region. S5. For logarithmic distance data, use the corresponding differential filtering strategies for different types of region-based resolution units to perform distance filtering; S6. Combine the results of all resolution units after differential filtering to obtain the complete logarithmic domain filtered distance data; S7. Perform an inverse exponential transform on the distance-filtered data to obtain the second distance data; perform matched filtering on the second distance data along the azimuth direction to output the enhanced image of the subsea pipeline.
[0005] As a further technical solution, the specific implementation of sound propagation compensation is as follows: Based on measured or historical sound velocity profile data, the sound ray bending correction amount at different distances and grazing angles is calculated using the ray tracing method. The ray tracing method discretizes the sound velocity profile into several layers based on Snell's law, calculates the sound ray propagation path and time delay layer by layer, obtains the time delay correction value at each distance gate, and performs time delay correction on the echo signal at each distance gate. Simultaneously, based on the relationship between the sound absorption coefficient of seawater and frequency, salinity, temperature, and depth, for example, the Thorp formula or Franc's formula is used. The OIS-Garrison acoustic absorption model calculates the high-frequency acoustic absorption attenuation factor related to propagation distance and performs frequency-dependent compensation on the amplitude of the echo signal to ensure that the echo signal after acoustic propagation correction maintains consistent acoustic scattering characteristics at different distances. During the compensation process, attenuation coefficients are calculated for different frequency components, and the compensated frequency components are recombined to form an amplitude-corrected echo signal. The calculated ray bending correction amount is recorded as a time delay offset matrix according to the resolution unit and transferred to the azimuth processing stage as the input parameter for phase compensation in azimuth matched filtering.
[0006] As a further technical solution, seabed sediment type information is introduced during the construction of the joint discriminant factor. Specifically, based on the seabed sediment classification results, the sediment type of the area where the current resolution unit is located is identified. This sediment classification result can be obtained in advance through multibeam echo intensity inversion, historical sediment database matching, or online clustering based on the echo statistical characteristics of this system. If the sediment type is silt or clay-like soft sediment, the neighborhood window size is expanded to enhance the statistical stability of speckle patterns. The expanded window size is 1.5 to 2 times the original window size, and this range depends on... The rationale is that soft substrates typically exhibit more uniform scattering characteristics, and the speckle field tends to be fully developed, requiring a larger window to obtain stable second-order statistical estimates. Based on the ratio of the resolution cell size to the correlation length in the SAS system, a window enlargement factor of 1.5 to 2 times increases the number of independent samples from approximately 30 to over 70, satisfying the stationarity requirements of speckle statistics. If the substrate type is rock or hard, the neighborhood window size is reduced to preserve detailed structure; the reduced window size is 0.5 to 0.8 times the original window size. This range is based on the fact that hard substrates often... Exhibiting non-uniform, strong scattering characteristics, the structural detail scale is often between 2 and 4 resolution units. An excessively large window size will smooth out key structural features; a reduction of 0.5 to 0.8 times allows the window size to match the structural feature scale, preserving edge information while maintaining the basic reliability of statistical estimation. Based on the preset scattering statistical relationship of the substrate type, the speckle noise intensity factor is corrected by multiplying the speckle noise intensity factor by the scattering modulation coefficient corresponding to the substrate type. The scattering modulation coefficient for soft substrates is less than 1, specifically ranging from 0.6 to 0. The factor of 0.8 is based on the fact that the echo intensity distribution of soft substrates is closer to the Rayleigh distribution, and the theoretical ratio of its variance to mean is about 0.523. Multiplying by a modulation coefficient of 0.6 to 0.8 can make the corrected factor close to this theoretical value. The scattering modulation coefficient corresponding to hard substrates is greater than 1, and the specific value range is 1.2 to 1.5. This is based on the fact that hard substrates have strong scattering points, the echo intensity distribution has a longer tail, and the ratio of variance to mean is usually greater than 1. A modulation coefficient of 1.2 to 1.5 can make the corrected factor consistent with the actual statistical characteristics. The substrate type is also coded as the environmental modulation factor.
[0007] As a further technical solution, the environmental modulation factor is a dimensionless weight value set according to the substrate type. For soft substrates, the environmental modulation factor is between 0.3 and 0.5, based on the fact that structural information is sparse in soft substrate regions, and the environmental modulation factor should be reduced to decrease the contribution of structural features to the joint discriminant factor. A value of 0.3 to 0.5 can reduce the contribution of structural features in the joint discriminant factor to below 20%. For hard substrates, the environmental modulation factor is between 0.8 and 1.2, based on the fact that structural information is abundant in hard substrate regions, and the environmental modulation factor should be increased to enhance the identification weight of structural features. A value of 0.8 to 1.2 can increase the contribution of structural features to over 40%. The environmental modulation factor is then weighted and fused with the speckle noise intensity factor and the distance migration trajectory feature value. The fusion formula is: Joint discriminant factor = α·speckle noise The discriminant is calculated as follows: sound intensity factor + β·distance migration trajectory feature value + γ·environmental modulation factor, where α, β, and γ are preset weighting coefficients. The value of α is based on the fact that the speckle noise intensity factor ranges from 0 to 2, and a value of 0.4 normalizes its contribution to the range of 0 to 0.8. The value of β is based on the fact that the distance migration trajectory feature value is represented by the normalized curvature or offset, and its value ranges from 0 to 1. A value of 0.4 maintains a similar magnitude to the speckle factor. The value of γ is set to 0.2 in the soft substrate region to reduce the contribution of structural features, and to 0.8 in the hard substrate region to increase the contribution of structural features. This differentiated weighting is based on experimental observations that the signal-to-noise ratio of structural features in soft substrates is usually about 10 dB lower than that in hard substrates. Through weight adjustment, the discriminant factor is adaptively balanced under different substrate conditions, forming a joint discriminant factor.
[0008] As a further technical solution, the adaptive determination process of the global noise threshold set is as follows: The joint discriminant factors of all resolution units are constructed into a one-dimensional distribution sequence. Each numerical point is traversed on the numerical axis of the sequence with an initial step size. The initial step size can be set to a value between 0.01 and 0.05 based on the numerical range and resolution of the joint discriminant factors. After weighted fusion, the joint discriminant factors are distributed between 0 and 1.5. Traversing with a step size of 0.01 to 0.05 yields 30 to 150 candidate segmentation points, sufficient to ensure segmentation accuracy while controlling computational complexity. Each numerical point is used as... The internal dispersion of the sequences on both sides of the segmentation boundary is calculated, as well as the difference in dispersion between the sequences on both sides. The internal dispersion is measured by the sum of the variances or the sum of the standard deviations of the sequences on both sides after the segmentation, while the difference in dispersion between the sequences on both sides is measured by the absolute value of the difference between the means of the two sides. The ratio of the difference in dispersion to the sum of the internal dispersion is used as the segmentation evaluation index. The larger the value, the more compact the categories on both sides after the segmentation and the more significant the difference between categories, thus the better the segmentation effect. After traversal, the value point corresponding to the maximum value of the evaluation index is selected as the initial segmentation point.
[0009] The joint discriminant factors are divided into low-value and high-value groups based on the initial segmentation point. The cumulative distribution slopes of all factors in the low-value group and the high-value group are calculated separately. The cumulative distribution slope is calculated as the ratio of the difference between adjacent points on the cumulative distribution curve to the numerical increment at that point. Within the low-value group, the local rate of change of the cumulative distribution slope is calculated point by point from the minimum value towards the initial segmentation point. The local rate of change is the absolute value of the difference between the current point's cumulative distribution slope and the previous point's cumulative distribution slope, divided by the numerical interval between the two points. When the local rate of change first exceeds the preset slope abrupt change threshold, the joint discriminant factor value corresponding to that point is used as the first global noise threshold. This slope abrupt change threshold is set to a value between 0.5 and 1.5 based on experimental data. The range is determined by the fact that when the distribution of the joint discriminant factors transitions from a noise-dominated region to a flat region, the cumulative distribution slope will change significantly. The typical value of this rate of change is between 0.8 and 1.2, and is set between 0.5 and 1. A threshold of 0.5 can robustly capture the abrupt change without being disturbed by random fluctuations. Within the high-value group, the local rate of change of the cumulative distribution slope is calculated point by point from the initial segmentation point towards the maximum value. When the local rate of change first turns from positive to negative and continues to exceed the preset duration, the joint discriminant factor value corresponding to the inflection point is used as the second global noise threshold. The preset duration is set to maintain a negative rate of change for 3 to 5 consecutive points. The range of values is based on the fact that when the distribution of the joint discriminant factor transitions from a flat region to an edge texture region, the rate of change of the cumulative distribution slope will turn from positive to negative, but there may be small local fluctuations. Maintaining a negative value for 3 to 5 consecutive points can effectively filter out false inflections and ensure that the true distribution inflection point is captured. The first global noise threshold and the second global noise threshold constitute a global noise threshold set. The two thresholds correspond to the inflection point position of the transition from the noise-dominated region to the flat region and the inflection point position of the transition from the flat region to the edge texture region in the distribution of the joint discriminant factor, respectively.
[0010] As a further technical solution, the joint discriminant factor of the current resolution unit is compared with the first global noise threshold and the second global noise threshold respectively. If the joint discriminant factor is less than the first global noise threshold, it is temporarily labeled as a noise candidate unit; if it is between the two thresholds, it is temporarily labeled as a flat candidate unit; if it is greater than or equal to the second global noise threshold, it is temporarily labeled as an edge texture candidate unit. The temporary labeling results of all resolution units within the neighborhood window of the current resolution unit are obtained. The size of the neighborhood window is set to 7×7 or 9×9 pixels. In the SAS imaging system, the lateral width of the seabed pipeline structure is about 5 to 10 resolution units. If the neighborhood window is too small, it cannot cover enough structural context information; if it is too large, it will introduce too many irrelevant areas. The window size of 7×7 to 9×9 can achieve a balance between structural feature scale and computational efficiency. For resolution units temporarily labeled as noise candidate units or flat candidate units, if the proportion of edge texture candidate units in the neighborhood window exceeds the preset region consistency threshold, the labeling is corrected to edge texture. For candidate units, the consistency threshold for this region is set between 0.3 and 0.5. This range is based on the fact that when a resolving unit is temporarily labeled as a non-edge texture but 30% to 50% of its surrounding neighboring units are labeled as edge textures, the unit has a high probability of being an extension of an edge texture rather than isolated noise. A value of 0.3 to 0.5 can avoid over-correction or under-correction. For resolving units temporarily labeled as edge texture candidate units, if the proportion of noise candidate units in the neighborhood window exceeds the preset isolated point removal threshold, the unit is corrected to a flat candidate unit. This isolated point removal threshold is set between 0.6 and 0.8. This range is based on the fact that when 60% to 80% of the surrounding neighboring units of a resolving unit temporarily labeled as edge texture are labeled as noise, the unit is very likely a misjudgment caused by strong noise spikes rather than a true edge texture structure. A value of 0.6 to 0.8 can ensure the reliability of removal. The labeling result after neighborhood consistency correction is used as the final classification of the resolving unit.
[0011] As a further technical solution, for the resolving unit identified as the noise-dominant region, a first filtering window is constructed centered on the current resolving unit. The size of the first filtering window is set to a relatively large window between 11×11 and 15×15. The noise-dominant region requires sufficient statistical smoothing to suppress speckle. According to the principle of speckle statistics, the number of independent samples is inversely proportional to the reduction in variance. An 11×11 window contains approximately 121 samples, which can reduce the variance by nearly 10 times to approximately 8% of the original variance. A 15×15 window contains 225 samples, which can reduce the variance to approximately 4% of the original variance. The two sizes are selected based on the actual noise intensity. The logarithmic magnitude of all resolving units within the first filtering window is calculated. The mean value is used as the filter output value to fully suppress speckle noise. For the resolution cell identified as an edge texture area, the distance migration trajectory direction of the current resolution cell is extracted. This direction is determined by the azimuth parameter in the distance migration trajectory feature value. An asymmetric filter window is constructed along this direction. The extension length of the asymmetric filter window in the trajectory direction is set to 9 to 15 resolution cells, and the extension length in the vertical direction is set to 3 to 5 resolution cells. The edge texture area has strong structural continuity along the trajectory direction, requiring a longer extension length to ensure structural integrity. The azimuth continuity length of a typical pipe structure is between 10 and 20 resolution cells. Taking 9 to 15 resolution cells can cover the main structure. The vertical extension length should be... The wall thickness is smaller than the lateral dimension of the structure, typically between 3 and 5 resolution units. A thickness of 3 to 5 avoids cross-structure smoothing. The window's extension length in the trajectory direction is greater than its extension length in the vertical direction. The median of the set of sampling points within the window along the trajectory direction is used as the filtered output value to preserve edge continuity along the trajectory direction. For resolution units identified as flat regions, a second filtering window is constructed centered on the current resolution unit. The size of the second filtering window is set to a medium size between 5×5 and 9×9. Flat regions require moderate smoothing to reduce residual noise, but the window should not be too large to avoid introducing crosstalk from edge texture areas. A 5×5 window containing 25 samples can reduce the variance to the original variance. Approximately 4% of the original variance is reduced by 9×9 windows containing 81 samples, which can reduce the variance to approximately 1.2%, achieving a balance between noise suppression and boundary protection. The weighted mean of the logarithmic magnitudes of all resolving units within the window is calculated as the filtered output value, where the weighting coefficients are inversely proportional to the spatial distance of the sampling point from the current resolving unit. For example, the weighting coefficients adopt the form of a Gaussian function, and its standard deviation is set to one-quarter of the window size, so that the sampling points closer to the current resolving unit have higher weights. The inverse distance weighting effectively avoids excessive deviation from the center position during the smoothing process. The output results of the three filtering strategies are reorganized according to the spatial position of the corresponding resolving units to form the preliminary filtered distance data.
[0012] As a further technical solution, three filtered output values of the current resolution unit after processing by the noise-dominant region filtering strategy, the edge texture region filtering strategy, and the flat region filtering strategy are obtained respectively. These three filtered output values are obtained by pre-calculating and caching. The classification results of all resolution units within the neighborhood window of the current resolution unit are obtained, and the proportion of each of the noise-dominant region, edge texture region, and flat region within the neighborhood window is calculated. This proportion reflects the regional attribute distribution characteristics of the local environment where the current resolution unit is located. Using this proportion as a weighting coefficient, all filtered output values are weighted and summed. The weighted summation formula is: final output value = (noise region proportion × noise filtered output value) + (flat region proportion × flat filtered output value) + (edge texture region proportion × edge texture filtered output value). The summation result is used as the final output value of the current resolution unit. The final output values of all resolution units are arranged according to spatial position to form distance filtered data after logarithmic domain filtering.
[0013] As a further technical solution, second distance data is acquired, and the distance migration trajectory features and ray bending correction values recorded by each resolution unit in the distance processing stage are extracted. The distance migration trajectory features include the trajectory curve parameters of each resolution unit on the distance and azimuth planes, and the ray bending correction value is the value at the corresponding position in the aforementioned time delay offset matrix. A matched filter reference function is constructed along the azimuth direction. This reference function is the complex form of the standard synthetic aperture sonar azimuth matched filter function, specifically: Where j is the imaginary unit, λ is the wavelength of the acoustic wave, v is the velocity of the synthetic aperture sonar platform, t is the azimuth time, and R0 is the shortest slant range; a phase compensation amount corresponding to the range migration trajectory characteristics is introduced into the phase term of the reference function, and the expression for this phase compensation amount is: Where ΔR(τ) is the change in distance migration trajectory with slow time τ, and λ is the wavelength of the sound wave. The theoretical basis for this compensation is that the phase error caused by distance migration is a linear function of twice the wavenumber of the distance change. This compensation is introduced to decouple distance migration correction from azimuth focusing. At the same time, a residual phase error compensation term is calculated based on the ray bending correction. The residual phase error compensation term is a first-order correction term of the ray bending correction to the phase of the matched filter reference function, and its expression is: δt(τ) is the change in time delay deviation caused by ray bending with slow time τ, and c is the reference sound speed. The theoretical basis for this compensation is that ray bending leads to an equivalent slant range change, and the product of the time delay deviation and the sound speed is the equivalent slant range error. A first-order correction is used to effectively compensate for the linear component of this error. The above phase compensation and residual phase error compensation terms are superimposed on the phase term of the matched filter reference function, so that the matched filter process performs point-by-point phase correction based on the actual distance migration trajectory and ray bending effect experienced by each resolution unit.
[0014] As a further technical solution, after azimuth-matched filtering, the azimuth compression values of each resolution unit corresponding to the same range migration trajectory in the filtering result are coherently accumulated. The identification method for the same range migration trajectory is as follows: based on the range migration trajectory characteristics recorded in the range processing stage, resolution units with the same trajectory curve parameters are grouped into the same trajectory set. During accumulation, the classification result of each resolution unit is used as the selection condition. For resolution units that have been identified as edge texture areas, their compression values are included in the coherent accumulation. The theoretical basis is that edge texture areas contain real structural information, and their phases are coherent, so accumulation can enhance the structural signal. For resolution units in noise-dominated areas and flat areas, their compression values are not included in the final imaging accumulation. The theoretical basis is that the echo phases in these areas are random, and coherent accumulation not only cannot enhance the signal but will also introduce coherent speckle noise. However, their range filtering results are still retained for auxiliary analysis, which includes subsequent image quality assessment or seabed classification verification. The accumulation result is used as the final output value corresponding to the range migration trajectory to form the enhanced subsea pipeline image.
[0015] This invention provides an adaptive high-frequency signal filtering and enhancement method for SAS imaging of subsea pipelines. It has the following beneficial effects: 1. This invention constructs a joint discriminant factor that simultaneously includes speckle noise intensity factor and distance migration trajectory characteristics, thereby decoupling the distance migration effect of high-frequency sound propagation in the marine environment from the scattering characteristics of the seabed sediment at the signal level, and achieving accurate differentiation of noise-dominant areas, flat areas and edge texture areas under complex marine conditions.
[0016] 2. This invention adaptively determines the noise threshold set based on the global distribution characteristics of the joint discriminant factors, and after completing the region classification based on this threshold set, it adopts differentiated filtering strategies to achieve targeted suppression of speckle noise and selective preservation of pipeline structural features under different marine conditions.
[0017] 3. This invention achieves high-fidelity focusing imaging of subsea pipelines in environments with strong distance migration and non-uniform sound speed by transferring the distance migration trajectory features and acoustic ray bending correction amount obtained in the distance processing stage to the azimuth matched filtering stage, introducing point-by-point phase compensation in the reference function, and coherently gating and accumulating the compression values on the same trajectory based on the region classification results. Attached Figure Description
[0018] To make the content of this invention easier to understand, the invention will be further described in detail below with reference to specific embodiments and accompanying drawings, wherein: Figure 1 This is a schematic diagram of the process described in this invention. Detailed Implementation
[0019] The technical solutions in the embodiments of the present invention will be clearly and completely described below. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.
[0020] Example 1: Reference Figure 1 This invention provides a specific embodiment of an adaptive high-frequency signal filtering and enhancement method for SAS imaging of subsea pipelines. In this embodiment, the system first receives the original high-frequency echo signal, which is acquired by a synthetic aperture sonar system during subsea pipeline detection. For the original high-frequency echo signal, acoustic propagation compensation processing is performed: based on measured marine environmental sound velocity profile data, a ray tracing method based on Snell's law is used to discretize the sound velocity profile into several layers, and the sound ray propagation path and time delay are calculated layer by layer to obtain the sound ray bending correction amount at different distances and grazing angles. Time delay correction is then performed on the echo signal at each distance gate. Based on the acoustic absorption coefficient of the seawater medium, for example, using the Francois-Ga... The Rirson acoustic absorption model calculates a high-frequency acoustic absorption attenuation factor related to propagation distance. The relationship between this attenuation factor and frequency, salinity, temperature, and depth is determined by a built-in formula in the model. Frequency-dependent compensation is applied to the amplitude of the echo signal. The compensated frequency components are then recombined to form an echo signal corrected for acoustic propagation, ensuring that the signal maintains consistent acoustic scattering characteristics at different distances. Subsequently, the echo signal corrected for acoustic propagation is subjected to range-direction pulse compression to obtain the first distance data, which is then subjected to a homomorphic logarithmic transformation to obtain logarithmic domain distance data. During the acoustic propagation compensation process, the calculated ray bending correction is recorded as a time delay offset matrix according to the resolution unit and transferred to the azimuth processing stage.
[0021] Next, each resolution cell is traversed over the logarithmic distance data. For the current resolution cell, the seabed type is first identified based on the seabed classification results, which are pre-obtained through multibeam echo intensity inversion. If the seabed type is soft sediment, the neighborhood window size is increased from the default 7×7 to 11×11 to enhance speckle statistical stability. If the seabed type is hard rock, the neighborhood window size is reduced to 5×5 to preserve detailed structure. Subsequently, the local gray mean and local gray variance within the neighborhood window are calculated, and the speckle noise intensity factor of the resolution cell is constructed based on the ratio of local variance to local mean. The speckle noise intensity factor is corrected according to the preset scattering statistics of the seabed type: for soft seabed, the speckle noise intensity factor is multiplied by the scattering modulation coefficient 0. 7. To make the corrected factor close to the theoretical variance-to-mean ratio of the Rayleigh distribution, for hard substrates, multiply by a scattering modulation coefficient of 1.3 to make it consistent with actual statistical characteristics. Simultaneously, extract the distance migration trajectory features in the resolution cell, which are represented by normalized curvature. Encode the substrate type as an environmental modulation factor, with the environmental modulation factor for soft substrates set to 0.4 and the environmental modulation factor for hard substrates set to 1.0. Subsequently, the environmental modulation factor is weighted and fused with the corrected speckle noise intensity factor and the distance migration trajectory feature value. The fusion formula is: Joint discriminant factor = 0.4 × speckle noise intensity factor + 0.4 × distance migration trajectory feature value + γ × environmental modulation factor, where γ is set to 0.2 in the soft substrate region and 0.6 in the hard substrate region. Through this differentiated weight allocation, a joint discriminant factor is formed.
[0022] After obtaining the joint discriminant factors of all resolution units, the global noise threshold set is adaptively determined based on their distribution. Specifically, the joint discriminant factors of all resolution units are constructed as a one-dimensional distribution sequence, with values ranging from 0 to 1.5. The sequence is iterated along its numerical axis with an initial step size of 0.02. Using each numerical point as a dividing boundary, the internal dispersion of the sequences on both sides and the difference in dispersion between the sequences are calculated. The internal dispersion is calculated as the sum of the variances of the sequences on both sides, and the difference in dispersion between the sequences is calculated as the absolute value of the difference between their means. The ratio of the difference in dispersion to the sum of the internal dispersion is used as the segmentation evaluation index. After the iteration is complete, the numerical point corresponding to the maximum value of this evaluation index is selected as the initial dividing point. Using the initial dividing point as the boundary, the joint discriminant factors are divided into a low-value group and a high-value group, and the low-value group is calculated separately. The cumulative distribution slope of all factors in the group and the cumulative distribution slope of all factors in the high-value group are calculated. In the low-value group, the local rate of change of the cumulative distribution slope is calculated point by point from the minimum value to the initial segmentation point. When the local rate of change first exceeds the preset slope change threshold of 1.0, the joint discriminant factor value corresponding to that point is taken as the first global noise threshold. This threshold corresponds to the inflection point position in the joint discriminant factor distribution from the noise-dominated area to the flat area. In the high-value group, the local rate of change of the cumulative distribution slope is calculated point by point from the initial segmentation point to the maximum value. When the local rate of change first changes from a positive value to a negative value and continues to exceed the preset duration length of 4 points, the joint discriminant factor value corresponding to the inflection point is taken as the second global noise threshold. This threshold corresponds to the inflection point position from the flat area to the edge texture area. The first global noise threshold and the second global noise threshold constitute the global noise threshold set.
[0023] Next, the joint discriminant factor of each resolution unit is compared with the first global noise threshold and the second global noise threshold. If the joint discriminant factor is less than the first global noise threshold, it is temporarily labeled as a noise candidate unit; if it is between the two thresholds, it is temporarily labeled as a flat candidate unit; if it is greater than or equal to the second global noise threshold, it is temporarily labeled as an edge texture candidate unit. The temporary labeling results of all resolution units within the 7×7 neighborhood window of the current resolution unit are obtained. For a resolution unit temporarily labeled as a noise candidate unit or a flat candidate unit, if the proportion of edge texture candidate units in the neighborhood window exceeds the region consistency threshold of 0.4, its labeling is corrected to edge texture candidate unit. For a resolution unit temporarily labeled as an edge texture candidate unit, if the proportion of noise candidate units in the neighborhood window exceeds the isolated point elimination threshold of 0.7, its labeling is corrected to flat candidate unit. The labeling result after neighborhood consistency correction is used as the final classification of the resolution unit, dividing it into a noise-dominated region, an edge texture region, or a flat region.
[0024] For logarithmic distance data, differentiated filtering strategies are applied to different types of resolution units for range filtering. For resolution units identified as noise-dominant regions, a 13×13 first filtering window is constructed centered on the current resolution unit, and the local mean of the logarithmic magnitude of all resolution units within the first filtering window is calculated as the filtering output value. For resolution units identified as edge texture regions, the distance migration trajectory direction of the current resolution unit is extracted. This direction is determined by the azimuth parameter in the distance migration trajectory feature value. An asymmetric filtering window is constructed along this direction, with an extension length of 12 resolution units in the trajectory direction and an extension length of [missing information] in the vertical direction. The resolution is set to 4 units, ensuring that the window extends longer in the trajectory direction than in the vertical direction. The median of the set of sampling points within the window along the trajectory direction is used as the filtered output value. For a resolution unit determined to be a flat region, a 7×7 second filtering window is constructed centered on the current resolution unit. The weighted mean of the logarithmic magnitudes of all resolution units within the window is calculated as the filtered output value. For example, the weighting coefficients are in the form of a Gaussian function, with a standard deviation set to one-quarter of the window size, i.e., 1.75, so that the weighting coefficients are inversely proportional to the spatial distance of the sampling point from the current resolution unit. The outputs of the above three filtering strategies are then recombined according to the spatial position of the corresponding resolution units.
[0025] To further improve the filtering effect, three filtered output values of the current resolution unit are obtained after processing by the noise-dominated region filtering strategy, the edge texture region filtering strategy, and the flat region filtering strategy. The final classification results of all resolution units within the 7×7 neighborhood window of the current resolution unit are also obtained, and the proportion of each of the noise-dominated region, edge texture region, and flat region within the neighborhood window is calculated. Using this proportion as a weighting coefficient, all filtered output values are weighted and summed. The summation formula is: Final output value = (Noise region proportion × Noise filter output value) + (Flat region proportion × Flat filter output value) + (Edge texture region proportion × Edge texture filter output value). The summation result is used as the final output value of the current resolution unit. The final output values of all resolution units are arranged according to their spatial position to form distance filtered data after logarithmic domain filtering.
[0026] Subsequently, an inverse exponential transform is performed on the range-filtered data to obtain the second range data. After obtaining the second range data, the range migration trajectory features and ray bending correction values recorded by each resolution unit in the range processing stage are extracted. The range migration trajectory features include the trajectory curve parameters of each resolution unit on the range and azimuth planes, and the ray bending correction value is the value at the corresponding position in the aforementioned time delay offset matrix. A matched filter reference function is constructed along the azimuth direction. This reference function is the complex form of the standard synthetic aperture sonar azimuth matched filter function. Where j is the imaginary unit, λ is the wavelength of the acoustic wave, v is the velocity of the synthetic aperture sonar platform, t is the azimuth time, and R0 is the shortest slant range; a phase compensation amount corresponding to the range migration trajectory characteristics is introduced into the phase term of the reference function, and the expression for this phase compensation amount is: Where ΔR(τ) is the change in distance migration trajectory with slow time τ; simultaneously, a residual phase error compensation term is calculated based on the ray bending correction, which is a first-order correction term of the ray bending correction to the phase of the matched filter reference function, and its expression is: , where δt(τ) is the change of time delay deviation caused by ray bending with slow time τ, and c is the reference sound speed; the above phase compensation amount and residual phase error compensation term are superimposed on the phase term of the matched filter reference function, so that the matched filter process performs point-by-point phase correction based on the actual distance migration trajectory and ray bending effect experienced by each resolution unit.
[0027] After azimuth-matched filtering is completed, the azimuth compression values of each resolution unit corresponding to the same range migration trajectory in the filtering results are coherently accumulated. The identification method of the same range migration trajectory is as follows: based on the range migration trajectory characteristics recorded in the range processing stage, resolution units with the same trajectory curve parameters are grouped into the same trajectory set. During accumulation, the classification result of each resolution unit is used as the gating condition. For resolution units that have been identified as edge texture areas, their compression values are included in the coherent accumulation. For resolution units in noise-dominated areas and flat areas, their compression values are not included in the final imaging accumulation, but their range filtering results are still retained for auxiliary image quality assessment. The accumulation result is used as the final output value corresponding to the range migration trajectory to form the enhanced seabed pipeline image.
[0028] Example 2: This example compares the performance of the method of the present invention with that of the prior art using measured data. The experimental data comes from a synthetic aperture sonar detection mission of a submarine pipeline in a certain sea area. The seabed consists of soft sediment and hard rock. The pipeline is a steel pipeline with a diameter of 0.6m. The synthetic aperture sonar system operates at a center frequency of 150kHz, a signal bandwidth of 20kHz, an azimuth resolution of 0.05m, a range resolution of 0.04m, and a platform speed of 2.5kn. The raw echo data is processed by acoustic propagation compensation, range pulse compression, and homomorphic logarithmic transformation to obtain logarithmic domain distance data, which serves as the basic input for processing in this example.
[0029] Two existing technologies are set up for comparison. Method 1 is the traditional synthetic aperture sonar imaging processing, which only performs range pulse compression and standard azimuth matched filtering, without including region identification and differential filtering. Method 2 is a comparative method that uses a single global filter, which performs range processing on the logarithmic domain range data using a median filter with a fixed window size of 9×9, followed by standard azimuth matched filtering. The method of this invention strictly follows all the steps described in the technical solution, including joint discriminant factor construction, adaptive global noise threshold determination, region classification correction, differential filtering, and azimuth point-by-point phase correction and coherent gating accumulation. The key parameter values are the same as in Example 1.
[0030] Three typical regions were selected from the imaging results for quantitative evaluation: the pipe edge region (50 pixel pairs were selected perpendicular to the pipe orientation to calculate the edge preservation index), the flat seabed region near the pipe (uniform 100×100 pixel regions were selected in both soft and hard seabeds to calculate the equivalent number of views), and the pipe body region (20 consecutive resolution units were used to calculate the signal-to-noise ratio and azimuth focusing performance). The evaluation index was calculated using the following formula: Equivalent number of views (ENL) = μ² / σ², where μ is the mean pixel gray level in the flat region, and σ is the standard deviation of pixel gray level in that region; Edge preservation index: Among them I edge I represents the pixel values on both sides of the edge. background The background pixel value is represented by the subscript 'ref', indicating the ideal edge reference value; the signal-to-noise ratio (SNR) is the ratio of the signal power in the target area to the noise power in the background area; the target-background contrast ratio is... ,in The average grayscale value of the target region pixels. The average grayscale value of the background region pixels is used; 5 pixels on each side of the center line of the pipe are selected as the target region in the direction of the pipe cross-section, and 10 pixels on each side of the outside of the pipe are selected as the background region; the azimuth impulse response width is determined by the azimuth profile of the pipe target -3dB main lobe width.
[0031] The comparison results of various evaluation indicators are shown in Table 1: Table 1. Comparison of Imaging Quality Evaluation Indicators for Different Methods
[0032] As shown in Table 1 above, the method of the present invention is significantly superior to the two existing methods in all evaluation indicators. Regarding speckle noise suppression, the equivalent number of views in the soft substrate region increases from 1.2 in Method 1 to 5.6, an improvement of 367%; in the hard substrate region, it increases from 0.9 to 4.7, an improvement of 422%. Regarding edge structure preservation, the edge preservation index increases from 0.42 to 0.89, an improvement of 112%. Regarding the signal-to-noise ratio, it increases from 14.2 dB to 22.8 dB, an improvement of 8.6 dB. Regarding azimuth focusing performance, the impulse response... The width was compressed from 0.11m to 0.06m, approaching the theoretical resolution of 0.05m, and the sidelobe level decreased from -10.2dB to -17.5dB. In the seabed-substrate interface region, the edge texture area of the present invention achieved a correct classification rate of 91.7%, which is significantly better than the comparative scheme that did not incorporate the present invention. Experimental results show that the present invention achieves accurate identification of regional attributes by combining discriminant factors, and effectively suppresses speckle noise and preserves the pipe edge structure by combining adaptive threshold segmentation and differential filtering, thus significantly improving the quality of synthetic aperture sonar imaging in complex seabed environments.
[0033] The above description is merely a specific embodiment of the present invention. Any feature disclosed in this specification, unless specifically stated otherwise, may be replaced by other equivalent or similar alternative features. All disclosed features, or steps in all methods or processes, except for mutually exclusive features and / or steps, may be combined in any way. Any non-essential additions or substitutions made by those skilled in the art based on the technical features of the present invention shall fall within the protection scope of the present invention.
Claims
1. An adaptive high-frequency signal filtering and enhancement method for SAS imaging of subsea pipelines, characterized in that, include: S1. Perform acoustic propagation compensation on the original high-frequency echo signal. Acoustic propagation compensation includes calculating the ray bending correction amount based on the sound velocity profile data in the marine environment and performing frequency and distance-related acoustic absorption compensation based on the sound absorption coefficient of the seawater medium to obtain the echo signal after acoustic propagation correction. Perform range-direction pulse compression processing on the echo signal after acoustic propagation correction to obtain the first distance data, and perform homomorphic logarithmic transformation on it to obtain logarithmic domain distance data. S2. Traverse each resolution cell in the logarithmic distance data. For the current resolution cell, calculate the local gray mean and local gray variance within its neighborhood window, and construct the speckle noise intensity factor of the resolution cell based on the ratio of the local variance to the local mean. At the same time, extract the distance migration trajectory features in the resolution cell, and perform weighted fusion of the speckle noise intensity factor and the distance migration trajectory feature value of each resolution cell to construct a joint discriminant factor. S3. Adaptively determine a global noise threshold set based on the joint discriminant factor of all resolution units; S4. Compare the joint discriminant factor of each resolution unit with the global noise threshold set, and classify the resolution unit into one of the following categories based on the comparison results: noise-dominated region, edge texture region, or flat region. S5. For logarithmic distance data, use the corresponding differential filtering strategies for different types of region-based resolution units to perform distance filtering; S6. Combine the results of all resolution units after differential filtering to obtain the complete logarithmic domain filtered distance data; S7. Perform an inverse exponential transform on the distance-filtered data to obtain the second distance data; perform matched filtering on the second distance data along the azimuth direction to output the enhanced image of the subsea pipeline.
2. The adaptive high-frequency signal filtering and enhancement method for SAS imaging of subsea pipelines according to claim 1, characterized in that: S1 includes: calculating the ray bending correction amount at different distances and grazing angles based on measured or historical sound velocity profile data using ray tracing methods, and performing time delay correction on the echo signal at each distance gate; simultaneously, calculating the high-frequency sound absorption attenuation factor related to the propagation distance based on the relationship between the sound absorption coefficient of the seawater medium and frequency, salinity, temperature and depth, and performing frequency-dependent compensation on the amplitude of the echo signal to ensure that the echo signal after sound propagation correction maintains consistent sound scattering characteristics at different distances; and recording the calculated ray bending correction amount according to the resolution unit and transmitting it to the azimuth processing stage.
3. The adaptive high-frequency signal filtering and enhancement method for SAS imaging of subsea pipelines according to claim 1, characterized in that: S2 includes: Based on the classification results of the seabed sediment, the sediment type of the current resolution unit is identified. If the sediment type is soft sediment such as silt or clay, the neighborhood window size is expanded to enhance the statistical stability of speckle. If the sediment type is rock or hard sediment, the neighborhood window size is reduced to preserve detailed structure. According to the preset scattering statistical relationship of the sediment type, the speckle noise intensity factor is corrected, and the sediment type is encoded as an environmental modulation factor. The environmental modulation factor is weighted and fused with the speckle noise intensity factor and the distance migration trajectory feature value. The weight coefficient of the environmental modulation factor reduces the contribution of structural features in the soft sediment region and increases the contribution of structural features in the hard sediment region, forming a joint discriminant factor.
4. The adaptive high-frequency signal filtering and enhancement method for SAS imaging of subsea pipelines according to claim 1, characterized in that: S3 includes: constructing a one-dimensional distribution sequence of the joint discriminant factors of all resolution units; traversing each numerical point on the numerical axis of the sequence with an initial step size; calculating the internal dispersion of the sequences on both sides and the difference in dispersion between the sequences on both sides after segmentation, using each numerical point as a dividing boundary; using the ratio of the difference in dispersion to the sum of the internal dispersion as a segmentation evaluation index; selecting the numerical point corresponding to the maximum value of the evaluation index as the initial dividing point after traversal; dividing the joint discriminant factors into a low-value group and a high-value group with the initial dividing point as the boundary; calculating the cumulative distribution slope of all factors in the low-value group and the cumulative distribution slope of all factors in the high-value group, respectively; and calculating the cumulative distribution slope of all factors in the low-value group from the minimum value towards the initial dividing point. The local rate of change of the cumulative distribution slope is calculated point by point in the direction. The joint discriminant factor value corresponding to the first time the local rate of change exceeds the preset slope change threshold is taken as the first global noise threshold. In the high value group, the local rate of change of the cumulative distribution slope is calculated point by point from the initial segmentation point to the maximum value. When the local rate of change first turns from positive to negative and continues to exceed the preset duration, the joint discriminant factor value corresponding to the inflection point is taken as the second global noise threshold. The first global noise threshold and the second global noise threshold constitute a global noise threshold set. The two thresholds correspond to the inflection point position of the transition from the noise-dominated area to the flat area in the joint discriminant factor distribution, and the inflection point position of the transition from the flat area to the edge texture area, respectively.
5. The adaptive high-frequency signal filtering and enhancement method for SAS imaging of subsea pipelines according to claim 4, characterized in that: The joint discriminant factor of the current resolution unit is compared with the first global noise threshold and the second global noise threshold. If the joint discriminant factor is less than the first global noise threshold, it is temporarily labeled as a noise candidate unit. If it is between the two thresholds, it is temporarily labeled as a flat candidate unit. If it is greater than or equal to the second global noise threshold, it is temporarily labeled as an edge texture candidate unit. The temporary labeling results of all resolution units in the neighborhood window of the current resolution unit are obtained. For a resolution unit temporarily labeled as a noise candidate unit or a flat candidate unit, if the proportion of edge texture candidate units in the neighborhood window exceeds the preset region consistency threshold, its labeling is corrected to edge texture candidate units. For a resolution unit temporarily labeled as an edge texture candidate unit, if the proportion of noise candidate units in the neighborhood window exceeds the preset isolated point removal threshold, its labeling is corrected to flat candidate units. The labeling result after neighborhood consistency correction is used as the final classification of the resolution unit.
6. The adaptive high-frequency signal filtering and enhancement method for SAS imaging of subsea pipelines according to claim 1, characterized in that: S5 includes: for a resolution unit determined to be a noise-dominant region, constructing a first filtering window centered on the current resolution unit, and calculating the local mean of the logarithmic domain amplitude of all resolution units within the first filtering window as the filtering output value; for a resolution unit determined to be an edge texture region, extracting the distance migration trajectory direction of the current resolution unit, constructing an asymmetric filtering window along this direction, wherein the extension length of the asymmetric filtering window in the trajectory direction is greater than the extension length in the vertical direction, and taking the median of the set of sampling points in the trajectory direction within the window as the filtering output value; for a resolution unit determined to be a flat region, constructing a second filtering window centered on the current resolution unit, wherein the size of the second filtering window is between the trajectory direction extension lengths of the first filtering window and the asymmetric filtering window, and calculating the weighted mean of the logarithmic domain amplitude of all resolution units within the window as the filtering output value, wherein the weighting coefficient is inversely proportional to the spatial distance of the sampling point from the current resolution unit; and recombining the output results of the three filtering strategies according to the corresponding resolution units.
7. The adaptive high-frequency signal filtering and enhancement method for SAS imaging of subsea pipelines according to claim 6, characterized in that: The filtered output values of the current resolution unit after processing by filtering strategies for noise-dominant regions, edge texture regions, and flat regions are obtained respectively. Obtain the classification results of all resolving units within the neighborhood window of the current resolving unit, and count the proportion of noise-dominated areas, edge texture areas, and flat areas within the neighborhood window; Using this proportion as a weighting coefficient, all filtered output values are weighted and summed, and the sum is used as the final output value of the current resolution unit. The final output values of all resolution units are arranged according to their spatial location to form the distance filtered data after logarithmic domain filtering.
8. The adaptive high-frequency signal filtering and enhancement method for SAS imaging of subsea pipelines according to claim 1, characterized in that: S7 includes: acquiring second distance data and extracting the distance migration trajectory features and ray bending correction amount recorded by each resolution unit in the distance processing stage; constructing a matched filter reference function along the azimuth direction, introducing a phase compensation amount corresponding to the distance migration trajectory features and calculating a residual phase error compensation term based on the ray bending correction amount in the phase term of the reference function, wherein the residual phase error compensation term is a first-order correction term of the ray bending correction amount to the phase of the matched filter reference function; and enabling the matched filtering process to perform point-by-point phase correction based on the actual distance migration trajectory and ray bending effect experienced by each resolution unit.
9. The adaptive high-frequency signal filtering and enhancement method for SAS imaging of subsea pipelines according to claim 8, characterized in that: After completing the azimuth-matched filtering, the azimuth compression values of each resolution cell corresponding to the same range migration trajectory in the filtering results are coherently accumulated. The classification results of each resolution cell are used as the gating condition during accumulation. For resolution cells that have been identified as edge texture areas, their compression values are included in the accumulation. For resolution cells in noise-dominated areas and flat areas, their compression values are not included in the final image accumulation, but their range filtering results are still retained for auxiliary analysis. The accumulation result is used as the final output value corresponding to the range migration trajectory to form the enhanced image of the seabed pipeline.