A method for calculating sea ice resistance of an ice crushing channel for first icebreaking of layer ice

Abstract

The application discloses a kind of ice channel sea ice resistance calculation methods for layer ice first ice breaking, including collecting and obtaining the ice distribution video image in ice channel, to obtain the ice area feature map;And according to the ice area feature map, obtain the fitting probability density function for obtaining the ice size probability density distribution histogram, the ice shape feature distribution and the ice density of the ice distribution in the ice channel, and then the ice channel sea ice resistance calculation optimization model about the ice channel sea ice parameter can be constructed;Based on the ice channel sea ice resistance calculation optimization model, realize the ice breaking pilot guidance of ship in ice channel;Solve the problem that existing method cannot accurately obtain the ice channel sea ice parameter of layer ice first ice breaking, and then cannot effectively guide the ice breaking pilot of ship in ice channel.

Description

Technical Field

[0001] The present invention relates to the field of polar shipping technology, and in particular to a method for calculating sea ice resistance in an ice-breaking channel for first icebreaking in layered ice. Background Art

[0002] As the global shipping industry continues to expand into polar ice regions, icebreaking pilotage is crucial as a navigation strategy to ensure smooth access to polar ports and waterways. Sea ice parameters such as size, shape, and density are key factors influencing ice resistance. This study, through actual ship tests and image analysis, quantified the dynamic variations in ice size, shape, and density in icebreaking channels. A correlation model was established between icebreaking speed and the dynamic variation in ice density, providing a basis for following ship safety, ice resistance simulation, and pilotage optimization.

[0003] In the study of ice resistance in shipbuilding and ocean engineering, ice pool tests are currently used to explore the effects of different sea ice parameters on ice resistance. These tests usually use the control variable method to set different parameters such as ice density, size distribution and shape to quantify their effect on ship navigation resistance. However, this type of research has obvious limitations. The ice conditions in the test are artificially set static states, which cannot reflect the dynamic evolution of broken ice in the actual icebreaking channel. The real-time changing characteristics of broken ice in the channel after icebreaking (such as retraction rate, secondary distribution, etc.) are ignored; the shape and size of the broken ice are artificially set, and the actual shape of the broken ice after icebreaking is not considered; no dynamic correlation model between icebreaking speed, ship maneuvering and broken ice parameters (size, density) has been established.

[0004] In the field of polar science, there have been many studies on the movement of natural ice fragments, mainly focusing on the drift and aggregation behavior of ice fragments under the influence of environmental forces such as wind and current. Although such studies have revealed the macroscopic mechanism of ice-sea-air interaction, they have the following shortcomings: (1) The evolution process of natural ice fragments is significantly different from the ice fragments in the channel after active fragmentation by icebreakers; (2) The scale is not matched. Natural studies mostly focus on large-scale ice areas (km level), while channel contraction involves smaller-scale ice fragment dynamics (10-100m level); (3) The engineering applicability is weak, and it is not coupled with the navigation parameters of ships, making it difficult to directly guide icebreaking pilotage practices. Summary of the Invention

[0005] The present invention provides a method for calculating sea ice resistance in an ice-breaking channel for first icebreaking in layer ice, so as to overcome the above technical problems.

[0006] In order to achieve the above object, the technical solution of the present invention is:

[0007] A method for calculating sea ice resistance in a broken ice channel for first icebreaking in layer ice comprises the following steps:

[0008] S1: Using a pre-set aerial monitoring platform, collect and obtain video images of the ice distribution within the icebreaking channel from a vertical perspective;

[0009] S2: Obtaining a characteristic map of the crushed ice area based on the crushed ice distribution video image;

[0010] S3: Obtain the power law distribution of ice size to quantify the degree of ice damage based on the ice crushing area characteristic map;

[0011] S4: Obtaining the area-weighted probability density function of ice crushing size according to the power law distribution of ice crushing size;

[0012] and obtaining a fitting probability density function for obtaining a probability density distribution histogram of the crushed ice size based on the probability density function of the crushed ice size;

[0013] S5: Obtaining the distribution of crushed ice shape characteristics according to the crushed ice area characteristic map;

[0014] S6: Obtain the crushed ice density of the crushed ice distribution in the icebreaking channel according to the crushed ice area characteristic map;

[0015] Based on the density of broken ice, a dynamic relationship model between the ship's icebreaking speed and the density of broken ice is established;

[0016] S7: Based on the dynamic relationship model, the density of broken ice is obtained in real time when the ship breaks through the ice layer for the first time at different speeds, and the probability density function is fitted. By combining the distribution of broken ice shape characteristics and broken ice size, an optimization model for calculating the sea ice resistance of the broken ice channel is obtained based on the sea ice parameters of the broken ice channel. Based on the optimization model for calculating the sea ice resistance of the broken ice channel, icebreaking pilotage guidance for ships in the broken ice channel is achieved.

[0017] The sea ice parameters of the broken ice channel include the size of the broken ice, the shape characteristics of the broken ice, and the density of the broken ice.

[0018] Furthermore, the S2 specifically includes the following steps:

[0019] S21: converting the crushed ice distribution video image into a plurality of single-frame images, and preprocessing each single-frame image to obtain a preprocessed image;

[0020] The preprocessing operation includes performing median filtering denoising on the single-frame image to obtain a filtered denoised image; performing contrast-limited adaptive histogram equalization on the filtered denoised image to obtain an equalized image; and performing adaptive sharpening on the equalized image to obtain a preprocessed image.

[0021] S22: performing unmarked sea ice instance segmentation on the pre-processed image based on an image segmentation model (SAM) to identify and obtain a sea ice feature map of sea ice edge contour features;

[0022] S23: Based on the Phansalkar adaptive threshold segmentation technology, the broken ice area in the sea ice feature map is extracted and a binary image of the broken ice area is obtained, and the binary image of the broken ice area is used as the broken ice area feature map.

[0023] Furthermore, the S3 specifically includes the following steps:

[0024] S31: Obtaining the crushed ice diameter in the crushed ice area characteristic map;

[0025] And the expression of crushed ice diameter is

[0026]

[0027] Where: S represents the area of ​​broken ice obtained by counting the number of foreground pixels of the sea ice edge contour; D f Indicates the diameter of crushed ice;

[0028] S32: Obtain a power function distribution model that satisfies the crushed ice size distribution according to the crushed ice diameter, and the expression of the power function distribution model is:

[0029] N(D f )=βD f -α

[0030] Where: N(D f ) represents the cumulative number distribution function of broken ice, that is, the diameter of broken ice per unit area is less than D f The cumulative number of broken ice; α represents the exponent of the power function distribution, and its value can reflect the size of the broken ice; β represents the parameter reflecting the total number of broken ice or the overall scale of its distribution;

[0031] S33: quantifying the distribution of ice fragment size using a power function distribution model to obtain a power law distribution of ice fragment size that quantifies the degree of ice fragment damage;

[0032] Methods for quantifying ice size distribution include:

[0033] The power function distribution model is logarithmized to convert it into a linear model, and the expression of the linear model is:

[0034] log N(D f )=logβ-αlog D f

[0035] Based on the least squares method, log N(D f ) and log D f A linear fit is performed to obtain the fitting slope -α, which is used as a parameter of the power-law distribution of ice fragment size to quantify the degree of ice fragment damage.

[0036] Furthermore, the S4 specifically includes the following steps:

[0037] S41: Based on the power law distribution of ice crush size, define a set d containing multiple ice crush size categories i , and i = 1...M, where M represents the number of size boxes and M = N 0.5 , N represents the total number of ice fragments identified from the power-law distribution of ice fragment sizes;

[0038] S42: Setting a constant interval width Δd, wherein the constant interval width Δd is the size range of the crushed ice divided by the number of size boxes;

[0039] The probability density distribution of crushed ice size is obtained according to the constant interval width Δd, and the probability density distribution function of crushed ice size is obtained as follows:

[0040]

[0041] Where: A represents the total area of ​​the broken ice pieces; a j represents the area of ​​the jth piece of crushed ice belonging to the i-th size box; p A (d i ) represents the probability density distribution function value of ice crush size;

[0042] S43: Using the least squares method, the probability density distribution of the crushed ice size is fitted with a Weibull function according to the elliptical contour characteristics of the crushed ice, and a fitting probability density function for obtaining a histogram of the probability density distribution of the crushed ice size is obtained, which is expressed as follows:

[0043]

[0044] Where: k represents the crushed ice shape parameter; λ represents the crushed ice scale parameter; min(D f ) represents the minimum diameter of broken ice; N(>L) represents the number of broken ice with a size greater than a certain size threshold L; d i represents the size set of crushed ice; N0 represents the total number of crushed ice.

[0045] Furthermore, in S5 , the crushed ice shape characteristic distribution is obtained according to the crushed ice area characteristic map;

[0046] The distribution of the shape characteristics of the crushed ice includes any one of a roundness characteristic, an elliptical contour characteristic, and a fitted polygonal characteristic of each crushed ice;

[0047] The formula for obtaining the roundness of the crushed ice is:

[0048]

[0049] Where: R represents the roundness of the crushed ice; D P P represents the diameter of the circle around the crushed ice;f Indicates the circumference of the crushed ice;

[0050] Obtain and count the elliptical contour features of each piece of crushed ice, that is, the size features of the crushed ice;

[0051] The formula for obtaining the elliptical contour feature of the crushed ice is:

[0052]

[0053]

[0054] Where: d max Indicates the maximum caliper diameter of crushed ice; d min represents the minimum caliper diameter of the crushed ice; (x, y) represents the coordinates of the pixel position of the crushed ice in the image, represents the coordinates of the centroid of the broken ice obtained by taking the average of all the pixel positions of the broken ice; μ 20 ,μ 02 ,μ 11 represents an intermediate variable and μ 00 Indicates the crushed ice area of ​​each;

[0055] The size of the crushed ice is reflected by the caliper diameter ratio of the crushed ice to obtain the shape characteristics of the crushed ice;

[0056] Method for obtaining the fitting polygon features:

[0057] Call the cv2.approx Poly DP() tool in Open CV to fit the polygon;

[0058] Based on the Douglas-Peucker algorithm, the crushed ice contour point set is simplified according to the crushed ice area feature map to fit and obtain the crushed ice fitting polygon features.

[0059] Furthermore, the step S6 specifically includes the following steps:

[0060] S61: Obtaining crushed ice density data within the icebreaking channel according to the crushed ice area characteristic map;

[0061] And the formula for obtaining crushed ice density data is:

[0062]

[0063] Where: C(t) represents the density of crushed ice at time t; A ice (t) represents the total area of ​​the crushed ice region in the crushed ice region characteristic map at time t; A total represents the total area of ​​the entire broken ice area feature map; t represents the length of time after the stern leaves the shooting area;

[0064] S62: Obtain the time series data of the ice floe concentration data at different times after t > 0 according to the ice floe concentration data, and smooth the curve by fitting with a third-order B-spline. Construct the ice floe concentration dynamic function according to the time series data;

[0065] And the expression of the ice floe concentration dynamic function is

[0066] When t ≤ T1, the ice floe concentration is C1(t) = Kt + C O ;

[0067] When T1 < t ≤ T2, the ice floe concentration is

[0068] When t > T2, the ice floe concentration is C3(t) = C F ;

[0069] where, T1 and T2 represent time thresholds; C1(t), C2(t), and C3(t) represent the values of the ice floe concentration C(t) at corresponding different times; C O represents the initial concentration constant; K represents the concentration change rate control coefficient; C F represents the peak height of the ice floe concentration; σ represents the concentration change rate control parameter;

[0070] S63: Obtain the calculation parameters of the sea ice concentration C(t) at different icebreaking speeds according to the ice floe shape feature distribution combined with the ice floe concentration dynamic function, which include K, C O , t1, T2, σ, and C F ;

[0071] Set the calculation parameters of the sea ice concentration C(t) as the parameters of a quadratic function of the ship speed v;

[0072] Then establish a quadratic function of the ship speed v and the calculation parameters of the sea ice concentration C(t) based on quadratic polynomial regression, and use it as the dynamic relationship model between the ship icebreaking speed and the icebreaking parameters.

[0073] Further, the S7 specifically includes the following steps:

[0074] S71: According to the dynamic relationship model, obtain the ice floe concentration in real time when the ship breaks ice for the first time facing stratified ice at different speeds, and obtain the initial sea ice resistance calculation model by combining the ice floe shape feature distribution and the ice floe concentration;

[0075] And the expression of the initial sea ice resistance calculation model is

[0076] R ice=f(d,v,G,C(t)...)=f1(d,v,C(t)...)×G

[0077]

[0078] Where: R ice represents the sea ice resistance; d represents the average size of broken ice; v represents the ship speed; G represents the distribution of broken ice shape characteristics, i.e., any one of the roundness characteristics, elliptical contour characteristics, and fitted polygon characteristics of the broken ice; C(t) represents the concentration of broken ice at time t; f(d, v, G, C(t)...) represents the broken ice resistance function considering the sea ice parameters of the broken ice channel; f1(d, v, C(t)...) represents the common sea ice resistance calculation function; p1, p2, p3, and p4 represent the area proportions of quadrilaterals, pentagons, polygons, and triangles in the broken ice fitting polygon, respectively.

[0079] S72: By introducing the fitting probability density function and combining it with the initial sea ice resistance calculation model, an optimization model for calculating the sea ice resistance of the broken ice channel is obtained based on the sea ice parameters of the broken ice channel;

[0080] The expression of the optimization model for calculating sea ice resistance in broken ice channel is:

[0081] R ice =∑ i ∑ j f(d i ,v,G,C(t j )...)p A (d i )C(t j )Δd

[0082] Where: C(t j ) represents t j The density of crushed ice at the moment; A (d i ) represents the probability density distribution function value of ice crush size;

[0083] The sea ice resistance of the broken ice channel is obtained based on the optimization model of sea ice resistance calculation in the broken ice channel, so as to provide icebreaking pilotage guidance for ships in the broken ice channel.

[0084] Beneficial effects: The present invention provides a method for calculating the sea ice resistance of an ice-breaking channel for the first icebreaking of layer ice. The method collects and obtains a video image of the distribution of ice fragments in the ice-breaking channel to obtain a characteristic map of the ice fragment area. According to the characteristic map of the ice fragment area, a fitting probability density function for obtaining a probability density distribution histogram of the size of ice fragments, a characteristic distribution of ice fragment shapes, and a density of ice fragments of the ice fragment distribution in the ice-breaking channel are obtained respectively, so as to construct an optimization model for calculating the sea ice resistance of an ice-breaking channel based on the sea ice parameters of the ice-breaking channel. The present invention uses image technology to quickly count the size and shape of ice blocks, and can effectively predict the contraction of the channel after the icebreaker breaks ice for the first time based on the optimization model for calculating the sea ice resistance of the ice-breaking channel, so as to provide a basis for icebreaking pilotage guidance and following safety of following ships, ice resistance simulation, and pilotage optimization in the ice-breaking channel. BRIEF DESCRIPTION OF THE DRAWINGS

[0085] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following is a brief introduction to the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative labor.

[0086] Figure 1 This is a flow chart of the method for calculating sea ice resistance in a broken ice channel for first icebreaking in layer ice according to the present invention;

[0087] Figure 2 This is a schematic diagram of recording the ice-breaking test video in this embodiment;

[0088] Figure 3 Schematic diagram of image preprocessing in this embodiment;

[0089] Figure 4 This is the SAM segmentation effect diagram in this embodiment;

[0090] Figure 5 This is the Phansalkar adaptive threshold segmentation map in this embodiment;

[0091] Figure 6 This is a power law distribution diagram of crushed ice size in this embodiment;

[0092] Figure 7 : This is the probability density distribution diagram of the crushed ice size in this embodiment;

[0093] Figure 8 This is a fitting diagram of the crushed ice roundness distribution in this embodiment;

[0094] Figure 9 This is a distribution fitting diagram of the crushed ice caliper diameter ratio in this embodiment;

[0095] Figure 10 This is a schematic diagram of the classification of crushed ice shape fitting in this embodiment;

[0096] Figure 11 is the probability histogram of the crushed ice shape distribution in this embodiment;

[0097] Figure 12 This is a time series diagram of crushed ice density in this embodiment;

[0098] Figure 13 This is a fitting diagram of crushed ice density in this embodiment;

[0099] Figure 14 This is a resistance result diagram of the introduction of the dynamic function of crushed ice density in this embodiment;

[0100] Figure 15 This is a resistance result diagram for introducing the probability density function of crushed ice size in this embodiment. DETAILED DESCRIPTION

[0101] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without making any creative efforts shall fall within the scope of protection of the present invention.

[0102] This embodiment provides a method for calculating sea ice resistance in a broken ice channel for the first icebreaking of layer ice. Figure 1 As shown, the following steps are included:

[0103] S1: Using a pre-set aerial monitoring platform, collect and obtain video images of the ice distribution within the icebreaking channel from a vertical perspective;

[0104] Specifically, this embodiment conducted a multi-condition icebreaking test on the ice layer. During the test, an icebreaking tugboat and a drone were used to collect data in collaboration. The test area was located in the Bohai Bay (approximately 40 degrees north latitude and 121 degrees east longitude). The environmental conditions met the target research requirements and provided an ideal test site for studying the characteristics of ice in the channel after the first breaking of the ice layer. Figure 2 As shown, this embodiment uses a multi-rotor drone equipped with a Zenmuse X7 camera to record the ice-breaking process at a resolution of 1080P and a frame rate of 30 frames per second, ensuring high-precision capture of the size and shape characteristics of the broken ice. The gimbal is calibrated to a vertical downward position, and the hovering height is fixed at 27m (error ±0.5m). A single frame image covers an area of ​​approximately 20×10m. 2 , the resolution is 1280×720px, record the whole process of icebreaking video and save it in mp4 format;

[0105] S2: Obtaining a characteristic map of the crushed ice area based on the crushed ice distribution video image;

[0106] like Figure 3 The specific steps shown include:

[0107] S21: In order to effectively improve the quality of the original image, the crushed ice distribution video image is converted into a number of single-frame images, and each single-frame image is preprocessed to obtain a preprocessed image;

[0108] The preprocessing operation includes performing median filtering denoising on the single-frame image to obtain a filtered denoised image; performing contrast-limited adaptive histogram equalization on the filtered denoised image to obtain an equalized image; and performing adaptive sharpening on the equalized image to obtain a preprocessed image.

[0109] S22: Based on the image segmentation model SAM (SegmentAnything Model), the pre-processed image is segmented into unmarked sea ice instances to identify and obtain a sea ice feature map of sea ice edge contour features; Figure 4 As shown, the identified sea ice edge outline is saved for subsequent calculation of ice break size and shape;

[0110] S23: Based on the Phansalkar adaptive threshold segmentation technology, the broken ice area in the sea ice feature map is extracted and a binary image of the broken ice area is obtained, and the binary image of the broken ice area is used as the broken ice area feature map, such as Figure 5 As shown;

[0111] S3: Obtain the power law distribution of ice size to quantify the degree of ice damage based on the ice crushing area characteristic map;

[0112] The specific steps include:

[0113] S31: Obtaining the crushed ice diameter in the crushed ice region characteristic map; wherein, in this embodiment, the crushed ice diameter is described by the diameter of a circle having the same area as the crushed ice;

[0114] The expression of the crushed ice diameter is:

[0115]

[0116] Where: S represents the area of ​​broken ice obtained by counting the number of foreground pixels of the sea ice edge contour; D f Indicates the diameter of crushed ice;

[0117] S32: Obtain a power function distribution model that satisfies the crushed ice size distribution according to the crushed ice diameter, and the expression of the power function distribution model is:

[0118] N(D f )=βDf -α

[0119] Where: N(D f ) represents the cumulative number distribution function of broken ice, that is, the diameter of broken ice per unit area is less than D f The cumulative number of crushed ice; α represents the exponent of the power function distribution, and its value can reflect the size of the crushed ice; β represents a parameter reflecting the total number of crushed ice or the overall scale of the distribution; this embodiment can calculate the number of crushed ice in each interval by grouping the original data by size, and make a cumulative number of crushed ice curve, such as Figure 6 As shown;

[0120] S33: quantifying the distribution of ice fragment size using a power function distribution model to obtain a power law distribution of ice fragment size that quantifies the degree of ice fragment damage;

[0121] Methods for quantifying ice size distribution include:

[0122] The power function distribution model is logarithmized to convert it into a linear model, and the expression of the linear model is:

[0123] log N(D f )=log β-αlog D f

[0124] Based on the least squares method, log N(D f ) and log D f Performing a linear fit to obtain a fitting slope -α, and using the fitting slope -α as a parameter of the power law distribution of ice fragment size to quantify the degree of ice fragment damage;

[0125] S4: obtaining an area-weighted probability density function of crushed ice size according to a power-law distribution of crushed ice size; and obtaining a fitting probability density function for obtaining a probability density distribution histogram of crushed ice size based on the probability density function of the crushed ice size;

[0126] In order to more clearly reflect the characteristics of the crushed ice size distribution, this embodiment considers the distribution of crushed ice with a large area and a small number. The probability density function is used to describe the crushed ice size distribution, and a probability density distribution histogram of the crushed ice size is established according to the standardized process. Figure 7 As shown;

[0127] The specific steps include:

[0128] S41: Based on the power law distribution of ice crush size, define a set d containing multiple ice crush size categories i , and i = 1...M, where M represents the number of size boxes and M = N 0.5N represents the total number of crushed ice identified from the power-law distribution of crushed ice size. The number of size bins is described as follows: size categories are size intervals divided according to the size of floating ice. The number of size bins is determined by the number of bins. In this embodiment, the number of size bins is determined based on the total amount of floating ice, so that the number of ice cubes in each bin remains consistent.

[0129] S42: Setting a constant interval width Δd, where the constant interval width Δd is the size range of the crushed ice divided by the number of size bins, thereby ensuring that the width of each interval is the same and avoiding interference with the results due to differences in partition settings;

[0130] The probability density distribution of crushed ice size is obtained according to the constant interval width Δd, and the probability density distribution function of crushed ice size is obtained as follows:

[0131]

[0132] Where: A represents the total area of ​​the broken ice pieces; a j represents the area of ​​the jth piece of crushed ice belonging to the i-th size box; p A (d i ) represents the probability density distribution function value of ice crush size;

[0133] S43: Using the least squares method, the probability density distribution of the crushed ice size is fitted with a Weibull function according to the elliptical contour characteristics of the crushed ice, and a fitting probability density function for obtaining a histogram of the probability density distribution of the crushed ice size is obtained, which is expressed as follows:

[0134]

[0135] Where: k represents the crushed ice shape parameter; λ represents the crushed ice scale parameter; min(D f ) represents the minimum diameter of broken ice; N(>L) represents the number of broken ice with a size greater than a certain size threshold L; d i represents the size set of crushed ice; N0 represents the total number of crushed ice;

[0136] In this embodiment, by establishing an area-weighted probability density distribution histogram of crushed ice size after Weibull function fitting, fitting parameters k and λ can be obtained, and parameters k and λ can reflect the distribution of crushed ice size; and the probability density function parameters of the Weibull function fitting are shown in Table 1 below;

[0137] Table 1. Probability density function parameters

[0138]

[0139] S5: Obtaining the distribution of crushed ice shape characteristics according to the crushed ice area characteristic map;

[0140] In a specific embodiment, in S5, a characteristic distribution of crushed ice shape is obtained according to the characteristic map of the crushed ice region;

[0141] The distribution of the shape characteristics of the crushed ice includes any one of a roundness characteristic, an elliptical contour characteristic, and a fitted polygonal characteristic of each crushed ice;

[0142] The formula for obtaining the roundness of the crushed ice is:

[0143]

[0144] Where: R represents the roundness of the crushed ice; D P P represents the diameter of the circle around the crushed ice; f Indicates the circumference of the crushed ice;

[0145] Obtain and count the elliptical contour features of each crushed ice, that is, the crushed ice size features, such as Figure 8 As shown; the formula for obtaining the elliptical contour feature of the crushed ice is

[0146]

[0147] Where: d max Indicates the maximum caliper diameter of crushed ice; d min represents the minimum caliper diameter of the crushed ice; (x, y) represents the coordinates of the pixel position of the crushed ice in the image, represents the coordinates of the centroid of the broken ice obtained by taking the average of all the pixel positions of the broken ice; μ 20 ,μ 02 ,μ 11 represents an intermediate variable and μ 00 Indicates the crushed ice area of ​​each;

[0148] like Figure 9 As shown, the size of the crushed ice shape is reflected by the caliper diameter ratio of the crushed ice to obtain the crushed ice shape characteristics;

[0149] In this embodiment, although the shapes of the crushed ice vary, their shape information such as perimeter and area have certain statistical characteristics. In this embodiment, the roundness R and the ratio of the maximum and minimum caliper diameters d are used to calculate the roundness R and the maximum and minimum caliper diameters d. max / d min To describe the shape of crushed ice, as shown in Table 2;

[0150] Table 2. Roundness, maximum and minimum caliper diameter ratio

[0151]

[0152]

[0153] Method for obtaining the fitting polygon features:

[0154] Call the cv2.approxPolyDP() tool in OpenCV for fitting polygons;

[0155] And based on the Douglas - Peucker algorithm, simplify the ice floe contour point set according to the ice floe area feature map to fit and obtain the ice floe fitting polygon features; in this embodiment, by fitting the ice floe shape with different geometric models (such as triangles, quadrilaterals, pentagons, polygons), as Figure 10 shown, obtain the shape distribution probability histogram, as Figure 11 shown;

[0156] S6: Obtain the ice floe density of the ice floe distribution in the icebreaking channel according to the ice floe area feature map;

[0157] To establish a dynamic relationship model between the ship's icebreaking speed and the ice floe density according to the ice floe density;

[0158] Specifically, it includes the following steps:

[0159] S61: Obtain the ice floe density data in the icebreaking channel according to the ice floe area feature map;

[0160] And the acquisition formula of the ice floe density data is

[0161]

[0162] In the formula: C(t) represents the ice floe density at time t; A ice (t) represents the total area of the ice floe area in the ice floe area feature map at time t; A total represents the total area of the entire ice floe area feature map; t represents the time length after the ship's stern leaves the shooting area;

[0163] S62: According to the ice floe density data, obtain the time - series data of the ice floe density data at different times after t > 0, and smooth the curve by fitting with a third - order B - spline, and construct a dynamic function of the ice floe density according to the time - series data;

[0164] Specifically, as Figure 12 shown, in this embodiment, the ice floe density data at different times are sorted into a time series, and the curve is smoothed by fitting with a third - order B - spline. The derivative function of the smoothed curve is the change rate of the ice floe density; by analyzing the characteristics of the change rate curve, the curve is divided into three piece - wise functions. The first piece - wise function is a constant, and the time value range is t ≤ T1; the second piece - wise function is a Gaussian function, and the time value range is T1 < t ≤ T2; the value of the third piece - wise function is zero, and the time value range is t > T2. And the second piece - wise function is fitted into a generalized Gaussian function Integrate the generalized Gaussian function to obtain the original function Based on this, it is used to establish the dynamic function of ice floe concentration; where α represents the ice floe shape parameter; β represents the change rate parameter; C F represents the peak height of the concentration; σ represents the change rate control parameter; μ represents the position of the concentration peak;

[0165] and the expression of the ice floe concentration dynamic function is

[0166] When t ≤ T1, the concentration increases linearly, and the ice floe concentration is C1(t) = Kt + C O ;

[0167] When T1 < t ≤ T2, the concentration follows a Gaussian distribution, and the ice floe concentration is

[0168]

[0169] When t > T2, the concentration remains stable, and the ice floe concentration is C3(t) = C F ;

[0170] where T1 and T2 represent time thresholds; C1(t), C2(t), and C3(t) represent the values of the ice floe concentration C(t) at different times; C O represents the initial concentration constant; K represents the concentration change rate control coefficient; C F represents the peak height of the ice floe concentration; σ represents the concentration change rate control parameter;

[0171] S63: As<00​​​​​​​​​​​​​​​​​​​​

[0177]

[0178] S7: Based on the dynamic relationship model, the density of broken ice is obtained in real time when the ship breaks through the ice layer for the first time at different speeds. The ice resistance calculation optimization model for the broken ice channel is obtained based on the sea ice parameters of the broken ice channel by combining the distribution of broken ice shape characteristics and broken ice size according to the fitted probability density function. Based on the sea ice resistance calculation optimization model for the broken ice channel, icebreaking pilotage guidance for the ship in the broken ice channel is achieved. The sea ice parameters of the broken ice channel include broken ice size, broken ice shape characteristics, and broken ice density.

[0179] The steps in this embodiment are intended to convert the probability density distribution p based on area weighting into A (d i ), the distribution of crushed ice shape characteristics and crushed ice density are introduced into ice resistance calculation to make it more consistent with actual working conditions and reduce errors caused by uneven ice size, different shapes, and density changes over time;

[0180] The specific steps include:

[0181] S71: Based on the dynamic relationship model, the density of broken ice when the ship first breaks through the ice layer at different speeds is obtained in real time, and the initial sea ice resistance calculation model is obtained by combining the distribution of broken ice shape characteristics and the density of broken ice;

[0182] And the expression of the initial sea ice resistance calculation model is

[0183] R ice =f(d,v,G,C(t)...)=f1(d,v,C(t)...)×G

[0184]

[0185] Where: R ice represents the sea ice resistance; d represents the average size of the broken ice; v represents the ship speed; G represents the distribution of the broken ice shape characteristics, i.e., any one of the roundness characteristics, elliptical contour characteristics, and fitted polygon characteristics of the broken ice; C(t) represents the concentration of broken ice at time t; f(d, v, G, C(t)...) represents the broken ice resistance function considering the sea ice parameters of the broken ice channel; f1(d, v, C(t)...) represents a common sea ice resistance calculation function, "..." represents other parameters in the common sea ice resistance function, such as those derived from ice water tank tests or numerical calculation results; p1, p2, p3, and p4 represent the area proportions of quadrilaterals, pentagons, polygons, and triangles in the broken ice fitting polygon, respectively; and the method for obtaining the common sea ice resistance calculation function is an existing well-known technical means and will not be elaborated on here;

[0186] S72: By introducing the fitting probability density function and combining it with the initial sea ice resistance calculation model, an optimization model for calculating the sea ice resistance of the broken ice channel is obtained based on the sea ice parameters of the broken ice channel;

[0187] The expression of the optimization model for calculating sea ice resistance in broken ice channel is:

[0188] R ice =∑ i ∑ j f(d i ,v,G,C(t j )...)p A (d i )C(t j )Δd

[0189] Where: C(t j ) represents t j The density of crushed ice at the moment; A (d i ) represents the probability density distribution function value of ice crush size;

[0190] The sea ice resistance of the broken ice channel is obtained based on the optimization model of sea ice resistance calculation in the broken ice channel, so as to provide icebreaking pilotage guidance for ships in the broken ice channel.

[0191] The resistance calculation results in this implementation are as follows Figure 14 and 15 As shown, Figure 14 Other parameters used in the common sea ice resistance function include ship width 8m, channel width 10m, ice density 900kg / m 3 , ice thickness 0.3m, ship speed 3m / s; Figure 15 Other parameters used in the common sea ice resistance function include ship width 8m, channel width 10m, ice density 900kg / m 3 , ice thickness 0.3m, ship speed 3m / s. The size and shape of the broken ice are experimental data, and the ice density is the convergent value of the ice density after icebreaking.

[0192] Beneficial effect: The method described in this embodiment collects and obtains video images of the distribution of broken ice in the icebreaking channel to obtain a characteristic map of the broken ice area; and according to the characteristic map of the broken ice area, obtains the fitting probability density function for obtaining the probability density distribution histogram of the broken ice size, the characteristic distribution of the broken ice shape and the crushed ice density of the broken ice distribution in the icebreaking channel, thereby being able to construct an optimization model for calculating the sea ice resistance of the icebreaking channel based on the sea ice parameters of the icebreaking channel; the present invention uses image technology to quickly count the size and shape of ice blocks, and can effectively predict the contraction of the channel after the icebreaker breaks the ice for the first time based on the optimization model for calculating the sea ice resistance of the icebreaking channel, thereby providing a basis for icebreaking pilotage guidance and following ship following safety, ice resistance simulation and pilotage optimization in the icebreaking channel.

[0193] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the above embodiments, or replace some or all of the technical features therein with equivalents. However, these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for calculating sea ice resistance in a broken ice channel for first icebreaking, characterized in that: The specific steps include: S1: Using a pre-set aerial monitoring platform, collect and obtain video images of the ice distribution within the icebreaking channel from a vertical perspective; S2: Obtaining a characteristic map of the crushed ice area based on the crushed ice distribution video image; S3: Obtain the power law distribution of ice size to quantify the degree of ice damage based on the ice crushing area characteristic map; S4: Obtaining the area-weighted probability density function of ice crushing size according to the power law distribution of ice crushing size; and obtaining a fitting probability density function for obtaining a probability density distribution histogram of the crushed ice size based on the probability density function of the crushed ice size; S5: Obtaining the distribution of crushed ice shape characteristics according to the crushed ice area characteristic map; S6: Obtain the crushed ice density of the crushed ice distribution in the icebreaking channel according to the crushed ice area characteristic map; Based on the density of broken ice, a dynamic relationship model between the ship's icebreaking speed and the density of broken ice is established; S7: Based on the dynamic relationship model, the density of broken ice is obtained in real time when the ship breaks through the ice layer for the first time at different speeds. Based on the fitted probability density function, the distribution of broken ice shape characteristics and broken ice size are combined to obtain the optimization model for calculating the sea ice resistance of the broken ice channel based on the sea ice parameters of the broken ice channel. Based on the optimization model for calculating the sea ice resistance of the broken ice channel, icebreaking pilotage guidance for ships in the broken ice channel is achieved. The sea ice parameters of the broken ice channel include the size of the broken ice, the shape characteristics of the broken ice and the density of the broken ice; The specific steps include: S71: Based on the dynamic relationship model, the density of broken ice when the ship first breaks through the ice layer at different speeds is obtained in real time, and the initial sea ice resistance calculation model is obtained by combining the distribution of broken ice shape characteristics and the density of broken ice; And the expression of the initial sea ice resistance calculation model is In the formula: In the formula: represents the sea ice resistance; Indicates the average size of crushed ice; Indicates the ship's speed; Represents the distribution of crushed ice shape features, i.e., any one of the roundness feature, elliptical contour feature, and fitted polygon feature of the crushed ice; express The density of crushed ice at each moment; represents the ice resistance function considering the sea ice parameters of the ice-breaking channel; represents the sea ice resistance calculation function; 、 、 、 They represent the area proportions of quadrilaterals, pentagons, polygons, and triangles in the crushed ice fitting polygon respectively; S72: By introducing the fitting probability density function and combining it with the initial sea ice resistance calculation model, an optimization model for calculating the sea ice resistance of the broken ice channel is obtained based on the sea ice parameters of the broken ice channel; The expression of the optimization model for calculating sea ice resistance in broken ice channel is: Where: express The density of crushed ice at each moment; represents the probability density distribution function value of ice crush size; The sea ice resistance of the broken ice channel is obtained based on the optimization model of sea ice resistance calculation in the broken ice channel, so as to provide icebreaking pilotage guidance for ships in the broken ice channel.

2. The method for calculating sea ice resistance in a crushed ice channel for first icebreaking in layer ice according to claim 1 is characterized in that: The S2 specifically includes the following steps: S21: converting the crushed ice distribution video image into a plurality of single-frame images, and preprocessing each single-frame image to obtain a preprocessed image; The preprocessing operation includes performing median filtering denoising on the single-frame image to obtain a filtered denoised image; performing contrast-limited adaptive histogram equalization on the filtered denoised image to obtain an equalized image; and performing adaptive sharpening on the equalized image to obtain a preprocessed image. S22: performing unmarked sea ice instance segmentation on the pre-processed image based on an image segmentation model (SAM) to identify and obtain a sea ice feature map of sea ice edge contour features; S23: Based on Adaptive threshold segmentation technology is used to extract the broken ice area in the sea ice feature map and obtain a binary image of the broken ice area, and the binary image of the broken ice area is used as the broken ice area feature map.

3. The method for calculating sea ice resistance in a crushed ice channel for first icebreaking in layer ice according to claim 2 is characterized in that: The S3 specifically includes the following steps: S31: Obtaining the crushed ice diameter in the crushed ice area characteristic map; And the expression of crushed ice diameter is Where: represents the area of ​​broken ice obtained by counting the number of foreground pixels of the sea ice edge outline; Indicates the diameter of crushed ice; S32: Obtain a power function distribution model that satisfies the crushed ice size distribution according to the crushed ice diameter, and the expression of the power function distribution model is: Where: It represents the cumulative number distribution function of broken ice, that is, the diameter of broken ice per unit area is less than Cumulative amount of crushed ice; The exponent representing the power function distribution, the value of which can reflect the size of the ice fragments; represents a parameter reflecting the overall size of the total amount or distribution of crushed ice; S33: quantifying the distribution of ice fragment size using a power function distribution model to obtain a power law distribution of ice fragment size that quantifies the degree of ice fragment damage; Methods for quantifying ice size distribution include: The power function distribution model is logarithmized to convert it into a linear model, and the expression of the linear model is: Based on the least squares method, the linear model and Perform linear fitting and obtain the fitting slope , and the fitted slope Parameters of the power-law distribution of ice fragment size that quantify the extent of ice fragment damage.

4. The method for calculating sea ice resistance in a crushed ice channel for first icebreaking in layer ice according to claim 3 is characterized in that: The S4 specifically includes the following steps: S41: Based on the power law distribution of ice crush size, define a set containing multiple ice crush size categories ,and ,in Indicates the number of size boxes and , N represents the total number of ice fragments identified from the power-law distribution of ice fragment sizes; S42: Set constant interval width , and the constant interval width This is the size range of crushed ice divided by the number of size bins; And according to the constant interval width Obtain the probability density distribution of crushed ice size, and obtain the probability density distribution function of crushed ice size as follows: Where: A represents the total area of ​​the crushed ice fragments; Indicates that it belongs to i Size box j The area of ​​crushed ice; represents the probability density distribution function value of ice crush size; S43: Using the least squares method, the probability density distribution of the crushed ice size is fitted with a Weibull function according to the elliptical contour characteristics of the crushed ice, and a fitting probability density function for obtaining a histogram of the probability density distribution of the crushed ice size is obtained, which is expressed as follows: ， Where: k represents the shape parameters of crushed ice; λ represents the crushed ice scale parameter; Indicates the minimum crushed ice diameter; Indicates that the size is greater than a certain size threshold L The amount of crushed ice; A set of sizes representing crushed ice; Indicates the total amount of crushed ice.

5. The method for calculating sea ice resistance in a crushed ice channel for first icebreaking in layer ice according to claim 4 is characterized in that: In S5, the characteristic distribution of the crushed ice shape is obtained according to the characteristic map of the crushed ice area; The distribution of the shape characteristics of the crushed ice includes any one of a roundness characteristic, an elliptical contour characteristic, and a fitted polygonal characteristic of each crushed ice; The formula for obtaining the roundness of the crushed ice is: , = Where: Indicates the roundness of crushed ice; The diameter of the circle representing the perimeter of the crushed ice; Indicates the circumference of the crushed ice; Obtain and count the elliptical contour features of each piece of crushed ice, that is, the size features of the crushed ice; The formula for obtaining the elliptical contour feature of the crushed ice is: ， Where: Indicates the maximum caliper diameter of crushed ice; Indicates the minimum caliper diameter of crushed ice; represents the coordinates of the ice-crushed pixel locations in the image, represents the coordinates of the centroid of the broken ice obtained by taking the average of the coordinates of all broken ice pixels; , , represents an intermediate variable and , , ; Indicates the crushed ice area of ​​each; The size of the crushed ice is reflected by the caliper diameter ratio of the crushed ice to obtain the shape characteristics of the crushed ice; Method for obtaining the fitting polygon features: Call the cv2.approx Poly DP() tool in Open CV to fit the polygon; Based on the Douglas-Peucker algorithm, the crushed ice contour point set is simplified according to the crushed ice area feature map to fit and obtain the crushed ice fitting polygon features.

6. The method for calculating sea ice resistance in a crushed ice channel for first icebreaking in layer ice according to claim 4 is characterized in that: The S6 specifically includes the following steps: S61: Obtaining crushed ice density data within the icebreaking channel according to the crushed ice area characteristic map; And the formula for obtaining crushed ice density data is: Where: express The density of crushed ice at each moment; Indicates time The total area of ​​the crushed ice area in the crushed ice area characteristic map at each moment; The total area of ​​the region representing the entire broken ice area feature map; Indicates the length of time since the stern left the shooting area; S62: According to the crushed ice density data, obtain The time series data of crushed ice density at different moments are collected, and a smooth curve is fitted by third-order B-spline to construct a dynamic function of crushed ice density based on the time series data. And the expression of the dynamic function of crushed ice density is when , the crushed ice density is ; when , the crushed ice density is ; when , the crushed ice density is ; in, , Indicates the time threshold; Indicates the density of crushed ice at different times The value of represents the initial density constant; represents the density change rate control coefficient; Indicates the peak height of crushed ice density; represents the density change rate control parameter; S63: Based on the distribution of broken ice shape characteristics and the dynamic function of broken ice density, the sea ice density at different icebreaking speeds is obtained The calculation parameters include 、 、 、 、 as well as ; Sea ice concentration The calculation parameter is set as the ship speed v The number of parameters of the quadratic function; The ship speed is established based on the quadratic polynomial regression v and sea ice density The quadratic function of the calculation parameters is used as the dynamic relationship model between the ship's icebreaking speed and ice-breaking parameters.

Images