CN109801306B - Tidal ditch extraction method based on high resolution remote sensing image - Google Patents

Tidal ditch extraction method based on high resolution remote sensing image Download PDF

Info

Publication number
CN109801306B
CN109801306B CN201910059401.7A CN201910059401A CN109801306B CN 109801306 B CN109801306 B CN 109801306B CN 201910059401 A CN201910059401 A CN 201910059401A CN 109801306 B CN109801306 B CN 109801306B
Authority
CN
China
Prior art keywords
tidal
ditch
remote sensing
fine
resolution remote
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910059401.7A
Other languages
Chinese (zh)
Other versions
CN109801306A (en
Inventor
宫兆宁
王启为
周德民
张磊
汪星
井然
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Capital Normal University
Original Assignee
Capital Normal University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Capital Normal University filed Critical Capital Normal University
Priority to CN201910059401.7A priority Critical patent/CN109801306B/en
Publication of CN109801306A publication Critical patent/CN109801306A/en
Application granted granted Critical
Publication of CN109801306B publication Critical patent/CN109801306B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Image Analysis (AREA)
  • Image Processing (AREA)

Abstract

A tidal ditch extraction method based on high-resolution remote sensing images is disclosed, and comprises the following steps: and collecting a high-resolution remote sensing image, and eliminating atmospheric interference of the remote sensing image through radiometric calibration and atmospheric correction. The strategy of the method is to extract a wide ditch and a fine ditch separately, wherein the ditch with the width more than or equal to a preset pixel is the wide ditch, and the ditch with the width less than the preset pixel is the fine ditch. And (3) extracting the wide tidal trench by using a normalized water body index (NDWI) and a maximum inter-class variance method (OTSU). And (3) calculating the separability between classes by using a SEATH algorithm and a J-M (Jeffries-Matusita) distance, and selecting a wave band with larger difference between a tidal ditch and a tidal flat to extract a fine tidal ditch. Aiming at the wave band screened by the SEATH algorithm, after the improved fuzzy C-means algorithm is used for homogenizing heterogeneous background, the multi-scale Gaussian matched filtering is used for enhancing the Gaussian-shaped fine tidal creeks. The fine tidal furrow is then extracted using adaptive threshold segmentation based on global mean and standard deviation. The debris plaque in the fine sulcus is then removed. And finally, merging the small tidal ditches and the wide tidal ditches by adopting ditch logic OR operation to form complete tidal ditches.

Description

Tidal ditch extraction method based on high resolution remote sensing image
Technical Field
The invention relates to the technical field of remote sensing, in particular to a tidal ditch extraction method based on high-resolution remote sensing images.
Background
The tidal channel is a typical geomorphic factor on a tidal flat, is a tidal channel which develops in an intertidal zone (particularly a silt intertidal zone) and is formed by ocean power (particularly tidal action), is an important channel for continuously exchanging substances, energy and information with the outside of the tidal flat and the tidal channel system, plays a vital role in distributing tidal water and supplying silt, directly reflects the characteristics of the tidal flat, and is one of important parameters for researching the hydrological connectivity of coastal areas. With the development of construction engineering, mudflat culture and other development activities, the development of the tidal trench is greatly influenced. Therefore, the method has very important theoretical significance and practical application value for accurately extracting the spatial distribution characteristics of the tidal creeks.
Under the image of bad natural conditions, the accessibility of part of coastal areas is poor, in addition, the seawater erosion is serious, the change of the form of a tidal trench is frequent, and the traditional field measurement method cannot be effectively carried out. The tidal channel is characterized in that 1) the width of the tidal channel is different, the width of the tidal channel is different from tens of centimeters at the upper part of a few intertidal zones to thousands of meters at the sea entrance at the lower part of the intertidal zones, and the width difference is obvious. 2) The tidal flat background heterogeneity is strong, the tidal ditches develop from the lower part of the intertidal zone, and flow through three different wetland types of muddy tidal flats, shrub swamps and herbaceous swamps, and the water content of the tidal flats, the vegetation coverage and the sand content of the tidal ditches in the three areas are all greatly different. 3) The tidal channel is strongly anisotropic, and different segmentation orders may exhibit preferential directions that deviate from the random network. These features present a great challenge to the rapid and accurate extraction of the tidal furrow. The existing remote sensing image visual interpretation method has the defects of time and labor consumption and large subjectivity, and is not suitable for large-area repeated observation of a tidal ditch. In consideration of the persistent instability of the tidal creek development process, there is an urgent need to develop a technical method for monitoring the long-term dynamics of the tidal creek evolution process.
In terms of geometric morphology, the tidal canals are distributed in a winding linear shape, and are similar to the land river network, blood vessels and roads, many scholars have conducted intensive research on automatic or semi-automatic extraction of tidal canal-like elements such as roads, blood and land river network, and related extraction methods are generally classified into the following methods: flow model methods, supervised classification methods, object-oriented classification methods, and advanced image processing methods. The flow path model method depends on high-precision DEM data, is easily influenced by structural faults and cracks, can only obtain a tidal channel with a single-pixel width, and ignores the width information of the tidal channel; the supervision classification method is based on data driving, has poor anti-noise interference capability, cannot keep the spatial continuity of geographic elements, and inevitably generates salt and pepper phenomena; the object-oriented classification method is used for dividing an image into homogeneous areas aiming at a high-resolution image, but a large number of tiny tidal ditches developed in the areas cannot be effectively distinguished under the influence of the moisture content of tidal beaches in intertidal zones; the advanced image processing method focuses on the linear distribution characteristics of the tidal channel, is generally divided into three steps of linear enhancement, segmentation and post-processing, and produces good effects in many areas.
For dynamic monitoring of coastal tidal ditches, currently, a method of combining airborne LiDAR (light detection and ranging) DEM (digital elevation model) and aerial photogrammetry data is mostly adopted to extract the tidal ditches, but the high observation cost limits the wide application and popularization of the data in remote sensing observation of coastal zones. The satellite remote sensing image can provide long-time sequential data support for coastal tidal canal monitoring by virtue of the characteristics of high timeliness, large-scale repeated observation and relatively low price. Particularly, with the continuous progress of the space technology, the working wave band of the sensor is continuously expanded, the space resolution is also rapidly improved, and new vitality is injected into the dynamic remote sensing monitoring of the coastal zone.
The above information disclosed in this background section is only for enhancement of understanding of the background of the invention and therefore it may contain information that does not form the prior art that is already known in this country to a person of ordinary skill in the art.
Disclosure of Invention
Problems to be solved by the invention
As described above, the present invention needs to provide a tidal ditch extraction method based on high-resolution remote sensing images, and a tidal ditch extraction method based on digital image processing technology is proposed herein, which uses a multispectral image such as a domestic high-resolution second (GF-2) as a basic data source. Firstly, a wide tidal trench is extracted by utilizing a normalized water body index (NDWI) and a maximum inter-class variance method (OTSU). Secondly, the fine tidal ditches are enhanced on the basis of weakening the heterogeneity of the tidal flat background by using an improved fuzzy C-means algorithm (MFCM) and multi-scale Gaussian matched filtering (MGMF). Then extracting the fine tidal furrows by adaptive threshold segmentation. And finally, combining the fine tide ditches and the wide tide ditches to form a complete tide ditch network. The automatic extraction precision of the invention is obviously improved and the long-time dynamic continuous automatic extraction is realized.
Means for solving the problems
The present inventors have conducted intensive studies to achieve the above object, and more particularly, the present invention provides a tidal creek extraction method based on high resolution remote sensing images, comprising the steps of:
the first step, collecting high resolution remote sensing image containing multispectral data, eliminating atmospheric interference from the remote sensing image through radiometric calibration and atmospheric correction,
and a second step of extracting the wide and fine grooves separately, wherein the grooves with the width greater than or equal to a preset pixel are wide grooves, and the grooves with the width less than the preset pixel are fine grooves. And (3) calculating the separability between classes by using a SEATH algorithm and a J-M (Jeffries-Matusita) distance, and selecting a wave band with larger difference between a tidal ditch and a tidal flat to extract a fine tidal ditch. The SEATH algorithm is as follows:
Figure GDA0002645350900000021
the J-M distance formula is as follows:
J=2(1-e-B) Wherein, in the step (A),
j represents the distance between two classes, mi(i ═ 1,2) and σi(i ═ 1,2) respectively represent the mean and variance of the two classes,
thirdly, extracting the wide tidal trench based on the maximum inter-class variance method, and obtaining the optimal segmentation threshold k when the maximum inter-class variance is met*Wherein, the formula of the maximum inter-class variance method is as follows:
Figure GDA0002645350900000031
Figure GDA0002645350900000032
wherein the content of the first and second substances,
Figure GDA0002645350900000033
is the inter-class variance at threshold k, mGIs the average gray level of the whole image, m (k) is the gray level at threshold k, L is the number of gray levels of the image, P1(k) Is the probability that the pixel is classified into class 1,
a fourth step of homogenizing the mixtureBackground of matter, wherein the clustering centers v are initializediAnd an offset field betajUpdating the membership function muijV cluster centeriAnd an offset field betajWhen | | | vb+1-vb||<When the cluster center calculated in the two previous and subsequent times is smaller than a preset convergence threshold value, stopping calculation to obtain background normalized data, wherein i is a cluster, and j is a pixel and is the preset convergence threshold value; b and b +1 denote the number of iterations, vb+1Is the clustering center of iteration b +1 times, vbIs the cluster center for iteration b times,
and fifthly, enhancing a tiny tidal trench with a Gaussian shape by multi-scale Gaussian matching filtering, rotating the filter to obtain responses of the tidal trenches in different directions, and only keeping the maximum response values of the pixels in multiple directions. The response results of the matched filters of multiple scales are multiplied after normalization,
the sixth step, self-adaptive threshold segmentation is carried out on the fine tidal ditches, segmentation is carried out on the fine tidal ditches on the basis of the global mean value and the standard deviation, and the threshold formula is as follows: t ═ mean + k × std, where mean is the global mean; std is the global standard deviation; t is an optimal threshold value, the range of the parameter k is between 0.01 and 1, and pixels with values larger than T are divided into fine grooves, so that coarse fine groove segmentation results are obtained. Finally, removing the isolated patches with the area smaller than the set threshold value to obtain the final segmentation result of the fine tidal canals.
And a seventh step of combining the fine tide groove of the sixth step with the wide tide groove of the third step to form a complete tide groove.
In the tidal ditch extraction method based on the high-resolution remote sensing image, in the first step, multispectral data is collected by a multispectral camera, the multispectral data comprises a blue light wave band with the wavelength of 0.45-0.52 micrometers, a green light wave band with the wavelength of 0.52-0.59 micrometers, a red light wave band with the wavelength of 0.63-0.69 micrometers and a near infrared wave band with the wavelength of 0.77-0.89 micrometers, and the spatial resolution of the image is 4 m.
In the tidal canal extraction method based on the high-resolution remote sensing image, in the first step, an FLAASH atmospheric correction algorithm based on an MODTRAN radiation transmission model is selected to carry out atmospheric correction on the image.
In the tidal ditch extraction method based on the high-resolution remote sensing image, in the second step, the number of preset pixels is 5. The separability between classes is calculated by using a J-M (Jeffries-Matusita) distance through an SEATH algorithm, a wave band with a large difference between a tidal channel and a tidal flat is selected for fine tidal channel extraction, a GREEN wave band is selected in a mud flat area, an NDWI wave band is selected in a salt marsh area, wherein NDWI is (GREEN-NIR)/(GREEN + NIR)), GREEN refers to a GREEN wave band, and NIR refers to a near infrared wave band.
In the tidal canal extraction method based on the high-resolution remote sensing image, in the fourth step, the observed value r at the pixel j isjIs determined by the true value xjAnd an offset field betajObtained by addition, by subtracting an offset field betajWe can obtain the true reflectivity at pixel j, as follows:
rj=xjj,j∈[1,N]
the global objective function J is such that,
Figure GDA0002645350900000041
wherein C represents the number of cluster centers; n represents the sum of the number of pixels in the image; alpha controls the strength of the neighborhood effect;
Figure GDA0002645350900000042
representing the average value of the pixels in the neighborhood window of the pixel j; mu.sijA membership function representing that the pixel j belongs to the class i; m is a membership factor set to 2; v. ofiA cluster center representing class i; j is a global objective function, and when the global objective function reaches an extreme value, background homogenization data are obtained.
In the tidal trench extraction method based on the high-resolution remote sensing image, in the fifth step, a two-dimensional Gaussian matched filter is used for effectively enhancing the linear tidal trench, and the two-dimensional Gaussian matched filter is defined as follows:
Figure GDA0002645350900000043
where σ represents the distribution of intensity; l is the length of the convolution template along the y-axisThe degree is to smooth noise, x refers to the length of the x-axis of the convolution template; m (x, y) is a two-dimensional Gaussian matched filter; -g "(x) refers to a one-dimensional Gaussian matched filter
Figure GDA0002645350900000044
The derivative is taken negative in both x and y directions.
In the tidal canal extraction method based on the high-resolution remote sensing image, in the fifth step, the rotation interval of the filter is set to 10 degrees, and the filter rotates 18 times to cover all possible directions of the tidal canals. The length of the x-axis of the matched filter is set to | x | ═ 3 σ, and the width of the y-axis of the filter is set to 5 to detect a 5-pixel long tidal channel. The width of the tidal channel to be enhanced is between 1 and 5 pixels,
Figure GDA0002645350900000045
w is half the width of the tidal channel to be enhanced, and σ is the standard deviation of the second derivative of the gaussian function, representing the distribution of intensity. Observing the image to be enhanced, finding that the width of most of the tidal ditches to be enhanced is concentrated on 3 and 5 pixels, the optimal sigma is 0.9 and 1.4, and in order to simplify calculation, the sigma of the mud beach tidal ditch and the salt marsh tidal ditch is concentrated1And σ2Set uniformly to 1 and 1.5.
In the fourth step, after the cluster center selects and screens according to the original image histogram, the cluster center selects at the peak and the trough of the histogram by combining with the selected sample points. The MFCM algorithm is insensitive to α, which is set to 0.1. The cluster center of the mud flat area is [ 0.1601; 0.1846, setting the cluster center of the salt marsh area as [ 0.3198; -0.2344].
In the tidal ditch extraction method based on the high-resolution remote sensing image, in the sixth step of self-adaptive threshold segmentation, the parameter k is 0.2, and the size of the chip plaque removal is set to be 50.
The invention has the beneficial technical effects that:
the complex two-dimensional geometry of the tidal canal itself and the background heterogeneity of the tidal flat are key causes of difficulty in extracting the tidal canal from the optical imagery. Based on the width change characteristics and section Gaussian shape characteristics of the tidal current channel, the invention provides a set of coastal tidal current channel automatic extraction algorithm aiming at the high-resolution No. two 4m multispectral image: dividing a wide tide ditch by combining OTSU with water body index; extracting and segmenting the fine tide ditches by using the MFCM, the MGMF and the self-adaptive threshold; and combining the wide and small tidal ditches to form a complete tidal ditch network. The method overcomes the difficulty that the coastal tidal ditches have obvious scale variation under a heterogeneous background, and completely extracts the tidal ditches with different widths in the yellow river delta area. Quantitative precision evaluation shows that the overall precision of the method is more than 97%, the Kappa coefficient is more than 0.8, the misclassification error is less than 20%, and the misclassification error is less than 11%, so that the method is obviously improved compared with the traditional supervision and classification method, and fine tidal ditches are completely reserved while misclassification is reduced. Most parameters in the algorithm can be obtained through the self characteristics of the image, and the dependence on subjective experience is reduced on the premise of ensuring the automation level of the algorithm.
The above description is only an overview of the technical solutions of the present invention, and in order to make the technical means of the present invention more clearly apparent, and to make the implementation of the content of the description possible for those skilled in the art, and to make the above and other objects, features and advantages of the present invention more obvious, the following description is given by way of example of the specific embodiments of the present invention.
Drawings
Fig. 1 is a schematic step diagram illustrating a tidal ditch extraction method based on high-resolution remote sensing images according to an embodiment of the present invention.
Fig. 2 shows a flowchart of a tidal ditch extraction method based on high-resolution remote sensing images according to an embodiment of the present invention.
Fig. 3 shows a normalized background flow chart of the tidal channel extraction method based on high-resolution remote sensing images according to an embodiment of the present invention.
Fig. 4 shows a background uniformization display diagram of the tidal channel extraction method based on the high-resolution remote sensing image according to an embodiment of the invention.
Fig. 5 shows a cross-sectional view of a multi-scale enhanced tidal ditch based on a tidal ditch extraction method of high-resolution remote sensing images according to an embodiment of the invention.
Fig. 6 shows a multi-scale enhancement display diagram of the tidal channel extraction method based on the high-resolution remote sensing image according to an embodiment of the invention.
Fig. 7 shows an adaptive threshold relationship diagram of a tidal channel extraction method based on high-resolution remote sensing images according to an embodiment of the invention.
Fig. 8 is a schematic diagram illustrating an extraction result of the tidal ditch extraction method based on the high-resolution remote sensing image according to an embodiment of the present invention.
Fig. 9 is a schematic diagram illustrating an extraction result of a tidal ditch extraction method based on high-resolution remote sensing images according to another embodiment of the present invention.
Fig. 10 is a schematic diagram illustrating comparison of extraction results of a tidal ditch extraction method based on high-resolution remote sensing images and different methods according to an embodiment of the invention.
Detailed Description
Specific embodiments of the present invention will be described in more detail below with reference to the accompanying drawings. While specific embodiments of the invention are shown in the drawings, it should be understood that the invention may be embodied in various forms and should not be construed as limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the invention to those skilled in the art.
It should be noted that certain terms are used throughout the description and claims to refer to particular components. As one skilled in the art will appreciate, various names may be used to refer to a component. This specification and claims do not intend to distinguish between components that differ in name but not function. In the following description and in the claims, the terms "include" and "comprise" are used in an open-ended fashion, and thus should be interpreted to mean "include, but not limited to. The description which follows is a preferred embodiment of the invention, but is made for the purpose of illustrating the general principles of the invention and not for the purpose of limiting the scope of the invention. The scope of the present invention is defined by the appended claims.
For the purpose of facilitating an understanding of the embodiments of the present invention, the following description will be made in terms of several specific embodiments with reference to the accompanying drawings, and the drawings are not intended to limit the embodiments of the present invention.
Specifically, as shown in fig. 1, a schematic step diagram of a tidal ditch extraction method based on high-resolution remote sensing images is shown. A tidal ditch extraction method based on high-resolution remote sensing images comprises the following steps:
the first step S1, collecting a high resolution remote sensing image including multispectral data, the remote sensing image eliminating atmospheric interference through radiometric calibration and atmospheric correction,
in the second step S2, the sensing image is calibrated by radiation and corrected by atmosphere to eliminate atmospheric interference,
and a second step (S2) of separately extracting wide and fine grooves, wherein the grooves having a width equal to or greater than a predetermined number of pixels are wide grooves and the grooves having a width less than the predetermined number of pixels are fine grooves. And (3) calculating the separability between classes by using a SEATH algorithm and a J-M (Jeffries-Matusita) distance, and selecting a wave band with larger difference between a tidal ditch and a tidal flat to extract a fine tidal ditch. The SEATH algorithm is as follows:
Figure GDA0002645350900000061
the J-M distance formula is as follows:
J=2(1-e-B) Wherein, in the step (A),
j represents the distance between two classes, mi(i ═ 1,2) and σi(i ═ 1,2) respectively represent the mean and variance of the two classes,
the third step S3, extracting the wide ditch based on the maximum inter-class variance method, and obtaining the optimal segmentation threshold k when the inter-class variance is maximum*Wherein, the formula of the maximum inter-class variance method is as follows:
Figure GDA0002645350900000071
Figure GDA0002645350900000072
wherein the content of the first and second substances,
Figure GDA0002645350900000073
is the inter-class variance at threshold k, mGIs the average gray level of the whole image, m (k) is the gray level at threshold k, L is the number of gray levels of the image, P1(k) Is the probability that the pixel is classified into class 1,
a fourth step S4 of normalizing the heterogeneous background, wherein a clustering center v is initializediAnd an offset field betajUpdating the membership function muijV cluster centeriAnd an offset field betajWhen | | | vb+1-vb||<When the cluster center calculated in the two previous and subsequent times is smaller than a preset convergence threshold value, stopping calculation to obtain background normalized data, wherein i is a cluster, and j is a pixel and is the preset convergence threshold value; b and b +1 denote the number of iterations, vb+1Is the clustering center of iteration b +1 times, vbIs the cluster center for iteration b times,
in a fifth step S5, the multi-scale gaussian matched filter enhances the fine grooves with gaussian shape, the filter rotates to obtain the response of the grooves with different directions, and for each pixel, only the maximum response values of its multiple directions are retained. The response results of the matched filters of multiple scales are multiplied after normalization,
sixth step S6, adaptively thresholding the fine tidal trench, segmenting the fine tidal trench based on the global mean and standard deviation, the thresholding formula: t ═ mean + k × std, where mean is the global mean; std is the global standard deviation; t is an optimal threshold value, the range of the parameter k is between 0.01 and 1, and pixels with values larger than T are divided into fine grooves, so that coarse fine groove segmentation results are obtained. Finally, removing the isolated patches with the area smaller than the set threshold value to obtain the final segmentation result of the fine tidal canals.
In the seventh step S7, the fine tide trench of the sixth step S6 is combined with the logical or operation of the wide tide trench of the third step S3 to form a complete tide trench.
The tide ditch extraction method based on the high-resolution remote sensing image has excellent expansibility by combining with a tide ditch extraction algorithm of background homogenization and linear enhancement, and a clustering center can be selected at the wave crest and the wave trough of a histogram by combining with sample points after rough screening; parameters such as the size, the intensity distribution and the like of the Gaussian matched filter can be selected according to the width and the length of the tidal channel to be enhanced; the optimal adaptive threshold is calculated by a method of gradually increasing the step length, so that the applicability of the method to different sensor images and different areas is improved.
In a preferred embodiment of the method for extracting a tidal canal based on a high-resolution remote sensing image, in the first step S1, the multispectral data is collected by a multispectral camera, the multispectral data includes a blue light band with a wavelength of 0.45-0.52 micrometers, a green light band with a wavelength of 0.52-0.59 micrometers, a red light band with a wavelength of 0.63-0.69 micrometers, and a near infrared band with a wavelength of 0.77-0.89 micrometers, and the spatial resolution of the image is 4 m.
In the preferred embodiment of the tidal canal extraction method based on the high-resolution remote sensing image, in the first step S1, a FLAASH atmospheric correction algorithm based on a MODTRAN radiation transmission model is selected to perform atmospheric correction on the image.
In a preferred embodiment of the tidal canal extraction method based on the high-resolution remote sensing image, in the second step S2, the number of predetermined pixels is 5. The separability between classes is calculated by using a J-M (Jeffries-Matusita) distance through an SEATH algorithm, a wave band with a large difference between a tidal channel and a tidal flat is selected for fine tidal channel extraction, a GREEN wave band is selected in a mud flat area, an NDWI wave band is selected in a salt marsh area, wherein NDWI is (GREEN-NIR)/(GREEN + NIR)), GREEN refers to a GREEN wave band, and NIR refers to a near infrared wave band.
In a preferred embodiment of the tidal canal extraction method based on high-resolution remote sensing images, in the fourth step S4, the observed value r at the pixel j isjIs determined by the true value xjAnd an offset field betajObtained by addition, by subtracting an offset field betajWe can obtain the true reflectivity at pixel j, as follows:
rj=xjj,j∈[1,N]
the global objective function J is such that,
Figure GDA0002645350900000081
wherein C represents the number of cluster centers; n represents the sum of the number of pixels in the image; alpha controls the strength of the neighborhood effect;
Figure GDA0002645350900000082
representing the average value of the pixels in the neighborhood window of the pixel j; mu.sijA membership function representing that the pixel j belongs to the class i; m is a membership factor set to 2; v. ofiA cluster center representing class i; j is a global objective function, and when the global objective function reaches an extreme value, background homogenization data are obtained.
In a preferred embodiment of the tidal current extraction method based on high-resolution remote sensing images, in the fifth step S5, a two-dimensional gaussian matched filter is used to effectively enhance the linear tidal current, and the two-dimensional gaussian matched filter is defined as follows:
Figure GDA0002645350900000083
where σ represents the distribution of intensity; l is the length of the convolution template along the y-axis to smooth out noise, x refers to the length of the convolution template on the x-axis; m (x, y) is a two-dimensional Gaussian matched filter; -g' (x) refers to a one-dimensional Gaussian matched filter
Figure GDA0002645350900000084
The derivative is taken negative in both x and y directions.
In the preferred embodiment of the tidal current extraction method based on high-resolution remote sensing images, in the fifth step S5, the filter rotation interval is set to 10 °, and the filter is rotated 18 times to cover all possible directions of the tidal current. The length of the x-axis of the matched filter is set to | x | ═ 3 σ, and the width of the y-axis of the filter is set to 5 to detect a 5-pixel long tidal channel. The width of the tidal channel to be enhanced is between 1 and 5 pixels,
Figure GDA0002645350900000085
w is half the width of the tidal channel to be enhanced, and σ is the standard deviation of the second derivative of the gaussian function, representing the distribution of intensity. By observing the image to be enhanced, the majority are foundThe width of the tidal channel to be enhanced is concentrated into 3 and 5 pixels, the optimal sigma is 0.9 and 1.4, and the sigma of the tidal channel of the mud beach and the tidal channel of the salt marsh is calculated for simplification1And σ2Set uniformly to 1 and 1.5.
In a preferred embodiment of the tidal ditch extraction method based on the high-resolution remote sensing image, in the fourth step S4, after selecting and screening the cluster center according to the original image histogram, the cluster center is selected at the peak and the trough of the histogram by combining with the selected sample points. The MFCM algorithm is insensitive to α, which is set to 0.1. The cluster center of the mud flat area is [ 0.1601; 0.1846, setting the cluster center of the salt marsh area as [ 0.3198; -0.2344].
In a preferred embodiment of the tidal furrow extraction method based on the high-resolution remote sensing image, in the sixth step S6, in the adaptive threshold segmentation, the parameter k is 0.2, and the size of the debris plaque removal is set to 50.
To further understand the present invention, in an embodiment, high score second (GF-2) multispectral data is selected as research data, the imaging time is 2016, 8, 13 and the period is 8, 13 days, the tide level is low and the saline marsh vegetation is developed vigorously, and a reference can be provided for the influence of the subsequent saline marsh vegetation development on the tide furrow. The sensor using the image at this time is PMS2, and the wave band information and the spatial resolution of the sensor are shown in table 1. The preprocessing step includes radiometric calibration of the atmosphere to eliminate atmospheric disturbances and restore true surface reflectivity. The absolute calibration coefficient used in the preprocessing is obtained from a China resource satellite application center official network, the FLASSH atmospheric correction based on the MODTRAN model is selected as the atmospheric correction model, and the preprocessing step is completed in ENVI 5.3.
TABLE 1 high-resolution No. two image parameter information
Figure GDA0002645350900000091
Fig. 2 shows a flowchart of a tidal ditch extraction method based on high-resolution remote sensing images according to an embodiment of the present invention. The flow chart of the extraction of the wide and small grooves from the high-resolution second image is shown in FIG. 2. According to field investigation, the tidal furrows having a width of less than 20m are densely distributed and difficult to extract, and thus the tidal furrows having a width of less than 5 pixels are defined herein as fine tidal furrows. The fine tide ditch extraction steps are as follows: (1) the complex background is uniform, and the salt marsh tidal ditch area can better distinguish the tidal ditch from the tidal flat only by using the normalized water body index (NDWI, NDWI ═ GREEN-NIR)/(GREEN + NIR)). The method is characterized in that a SEaTH algorithm is adopted in the tidal bank area of the mud beach, the separability between classes is calculated by using J-M (Jeffries-Matusita) distance, and a green light wave band with large difference between the tidal bank and the tidal bank is selected. On the basis, the difference contrast between the target and the background caused by the background heterogeneity of the tidal flat is restrained by using an improved fuzzy C mean value algorithm (MFCM) aiming at the selected wave bands. (2) And multi-scale linear enhancement, namely enhancing the small tidal trench with a Gaussian shape by adopting multi-scale Gaussian matched filtering (MGMF) so as to solve the problem of width change of the tidal trench. (3) The filter is rotated to cover all possible directions, attenuating the effects of the strong aeolian anisotropy. (4) Fine tidal furrows were segmented using adaptive thresholds. (5) And removing the debris and the plaque. The wide tidal trench is filled with seawater and can be divided by NDWI and the maximum inter-class variance method (OTSU). And finally, combining the results of the small and wide tidal ditches by adopting logic OR operation to form a complete tidal ditch.
When water body extraction is carried out, the traditional hard clustering method is influenced by the uncertainty of the remote sensing image, so that the classification result is ambiguous. In fact, in the complex background of tidal beaches, the boundary between the water body and the non-water body is fuzzy, so that a fuzzy clustering method can be used to solve the ambiguity. The standard fuzzy C-means (FCM) algorithm is sensitive to noise and does not take into account the spatial distribution of pixels, which prevents it from being segmented in remote sensing images.
FIG. 3 shows a flow chart of a normalized background of a tidal channel extraction method based on high-resolution remote sensing images, which is to say, firstly, an observed value r at a pixel j is assumedjIs determined by the true value xjAnd an offset field betajObtained by addition, by subtracting an offset field betajWe can obtain the true reflectivity at pixel j, as follows:
rj=xjj,j∈[1,N]
then, the adjacent spatial information is integrated into the standard FCM algorithm to reduce the influence of isolated noise.
Figure GDA0002645350900000101
Wherein C represents the number of cluster centers; n represents the sum of the pixel numbers in the image; alpha controls the strength of the neighborhood effect;
Figure GDA0002645350900000102
representing the average value of the pixels in the neighborhood window of the pixel j; mu.sijA membership function representing that the pixel j belongs to the class i; m is a membership factor, typically set to 2; v. ofiA cluster center representing class i; j is the global objective function. And when the objective function reaches an extreme value, the image to be classified obtains the optimal classification scheme. For this purpose, a membership function mu is introducedijThe constraint of (2):
Figure GDA0002645350900000103
each parameter is solved using the lagrange multiplier method.
Figure GDA0002645350900000104
Fig. 4 is a background normalized display diagram of a tidal channel extraction method based on high-resolution remote sensing images according to an embodiment of the present invention, specifically referring to fig. 4, wherein (a) green light band original images; (b) green light band histogram; (c) a background homogenization result; (d) background normalized histogram.
The contrast between the fine tidal current and the image background is poor due to the influence of the longitudinal change of the tidal current and the mixed spectrum effect of the surrounding background pixels, the fine tidal current cannot be effectively identified only by using the OTSU, and the threshold value is reduced, so that more small tidal currents can be obtained, but more noise can be caused. Through observing the cross section of the tidal trench, the cross section of the tidal trench is mostly in an inverted U shape or an inverted V shape and can be approximated by a Gaussian curve, so that the linear tidal trench can be effectively enhanced by using a two-dimensional Gaussian matched filter. Fig. 5 is a cross-sectional view of a multi-scale enhanced tidal ditch of the tidal ditch extraction method based on high-resolution remote sensing images, which is shown in fig. 5.
The two-dimensional gaussian matched filter is defined as follows:
Figure GDA0002645350900000105
where σ represents a distribution of intensity whose value is closely related to the width of the target tidal channel; l is the length of the convolution template along the y-axis to smooth out noise.
In view of the strong anisotropy of the tidal canal system, it is necessary to rotate the filter by the possible angle to obtain the response of the tidal canal in different directions, and for each pixel, only the maximum response value of its multiple directions is retained. The rotation interval is set here to 10 °, the filter is rotated 18 times to cover all possible directions of the tidal channel. Considering that the width span of the tidal channel is very large, all tiny tidal channels cannot be enhanced by using a single-scale matched filter, and many tiny tidal channels can be submerged by Gaussian noise. To this end, a multi-scale enhancement method is adopted to solve this problem, and the response results of the matched filters of multiple scales are multiplied after normalization, so that the linear feature can be further enhanced while suppressing noise. Fig. 6 shows a multi-scale enhancement display diagram of the tidal channel extraction method based on the high-resolution remote sensing image according to an embodiment of the invention. Referring specifically to fig. 6, wherein (a) scale 1 enhances the results; (b) scale 2 enhancement results; (c) and (5) multi-scale enhancing the result.
In the fine tide trench segmentation process, the optimal threshold calculated by the OTSU is higher, so that a plurality of fuzzy fine tide trenches are ignored, for this reason, the fine tide trench is segmented by adopting an adaptive threshold combining a global mean and a standard deviation, and the threshold formula is as follows:
T=mean+k*std
where mean is the global mean; std is the global standard deviation; t is an optimal threshold; and dividing the pixels with the value larger than T into the tidal ditches.
In one embodiment, the parameters of the tidal ditch extraction algorithm are divided into three parts (1): a cluster center of an improved fuzzy C-means algorithm; (2) length, width and sigma of the multi-scale Gaussian matched filter; (3) threshold k in adaptive threshold segmentation.
And after selecting and screening the clustering center according to the original image histogram, selecting the clustering center at the peak and the trough of the histogram by combining with the selected sample points. The cluster center for the beach tidal channel herein is [ 0.1601; 0.1846, setting the cluster center of the salt marsh tidal current ditch as [ 0.3198; -0.2344 ]; the MFCM algorithm is insensitive to α, which is uniformly set to 0.1.
Because more than 99% of the area under the Gaussian curve is positioned in [ -3 sigma, 3 sigma [ -3 sigma [ ]]Thus, the length of the x-axis of the matched filter is set to be | x | ═ 3 σ, the width of the y-axis of the filter is set to be 5 to detect the 5-pixel long tidal channel, and the width of the 4 test zones to be enhanced to the tidal channel is between 1-5 pixels. In the multi-scale enhancement,
Figure GDA0002645350900000111
the optimum sigma corresponding to the tiny tidal ditches is between 0.3 and 1.4, the width of most of the tidal ditches to be enhanced is found to be concentrated into 3 and 5 pixels by observing images to be enhanced, the optimum sigma is between 0.9 and 1.4, and the sigma of the mud beach tidal ditches and the salt marsh tidal ditches is calculated for simplifying the calculation1And σ2Set uniformly to 1 and 1.5.
In the self-adaptive threshold segmentation, the setting of a parameter k has a large influence on the extraction precision, and how to keep balance between wrong segmentation and missed extraction is the key for selecting the optimal threshold, so that the range of the parameter k is set to be 0.01-1, 0.01 is used as a step length, and when a Kappa coefficient is the highest, the global optimal threshold is obtained. As shown in fig. 10, the maximum values of Kappa coefficients are all around 0.2 at 0.2, 0.25, 0.23 and 0.22 in the four test regions, and therefore k is set to 0.2. The dimensions of the debris plaque removal were set at 50. Fig. 7 shows an adaptive threshold relationship diagram of a tidal channel extraction method based on high-resolution remote sensing images, in which the relationship between the threshold k and the Kappa coefficient is shown in fig. 7.
In order to test the reliability of the method, a tidal channel extraction method combining complex background homogenization and multi-scale linear characteristic enhancement is compared with a Maximum Likelihood Method (ML) and a Support Vector Machine (SVM), and when a training sample is selected, the mud flat tidal channel is divided into four types, namely an intertidal zone, a fine tidal channel and a wide tidal channel; the salt marsh tidal ditches are divided into five categories, namely clouds, tidal beaches, vegetation, fine tidal ditches and wide tidal ditches. And (3) using a Radial Basis Function (RBF) as a kernel Function of the SVM, and selecting the highest Kappa coefficient from various parameter settings as a final result of the two comparison methods.
A Confusion Matrix (Consion Matrix) is calculated by a Consion Matrix-Using group Truth Image module of ENVI5.3, and the precision of the extraction result is evaluated. At this time, Overall Accuracy (OA), Kappa Coefficient (KC), leakage error (omision, EO) and error score error (Commission, EC) are selected as evaluation indexes.
The method is realized in a MATLAB 2017a platform, the complete extraction results of a mud beach tidal channel and a salt marsh tidal channel are shown in fig. 8 and 9, the quantitative evaluation of the extraction results of four methods in a local area is shown in fig. 10, and the quantitative evaluation of the extraction results of different methods in the local area is shown in table 2.
As can be seen by visual interpretation comparison, most of the tidal ditches are extracted more completely by the method, particularly the tiny tidal ditches with the width between 1 and 5 pixels, the extraction effect is better, and the integrity and the continuity of the tidal ditches are better maintained.
In the four test areas, ML and SVM are not ideal for fine tidal channel extraction from a visual interpretation perspective only, and all four test areas produce severe scoring artifacts. Especially in the area of mud beach gullies, a large number of beaches are misclassified as gullies.
Quantitative accuracy evaluation shows that the tidal ditch extraction method combining complex background homogenization and multi-scale linear characteristic enhancement provided by the invention is superior to other three methods in the aspects of Kappa coefficient, overall accuracy and leakage error, the Kappa coefficient is more than 0.8, and the matching degree of the extraction result of the method and ground actual measurement data is better.
In the tidal canal areas of the two mud beaches, the water content in the fine tidal canal is low, the difference between the fine tidal canal and the tidal beach is small in terms of spectrum, the phenomenon of foreign matters in the same spectrum is generated, and the extraction of linear ground objects and fine ground objects is seriously influenced. When a sample is selected, the degree of separation among classes is not ideal, the tiny tidal ditches cannot be effectively identified by supervision and classification, the tidal flat and the tiny tidal ditches cannot be effectively separated, and various evaluation indexes are generally low. From the geometrical aspect, the tidal flat is distributed in a block shape, the tidal ditch shows obvious linear characteristics, and after background interference is reduced by using the MFCM, the response of non-linear elements can be greatly inhibited by using linear characteristic enhancement, and errors are reduced.
TABLE 2 evaluation of the accuracy of tidal ditch extraction
Figure GDA0002645350900000121
Industrial applicability
The tidal ditch extraction method based on the high-resolution remote sensing image can be used in the field of remote sensing detection.
Although the embodiments of the present invention have been described above with reference to the accompanying drawings, the present invention is not limited to the above-described embodiments and application fields, and the above-described embodiments are illustrative, instructive, and not restrictive. Those skilled in the art, having the benefit of this disclosure, may effect numerous modifications thereto without departing from the scope of the invention as defined by the appended claims.

Claims (8)

1. A tidal ditch extraction method based on high-resolution remote sensing images comprises the following steps:
a first step (S1) of acquiring a high resolution remote sensing image including multispectral data, the remote sensing image being subjected to radiometric calibration and atmospheric correction to eliminate atmospheric interference,
a second step (S2) of separately extracting wide and fine grooves, wherein the grooves with the width greater than or equal to a predetermined number of pixels are wide grooves, the grooves with the width less than the predetermined number of pixels are fine grooves, calculating inter-class separability by using J-M distance by using SEATH algorithm, selecting wave bands to extract the fine grooves, and the SEATH algorithm formula is as follows:
Figure FDA0002645350890000011
the J-M distance formula is as follows:
J=2(1-e-B) Wherein, in the step (A),
j represents the distance between two classes, m1、m2Respectively representing the mean values, σ, of two classes1、σ2The variances of the two classes are represented separately,
a third step (S3) of extracting a wide ditch based on the maximum inter-class variance method, and obtaining an optimal segmentation threshold k when the maximum inter-class variance is satisfied*Wherein, the formula of the maximum inter-class variance method is as follows:
Figure FDA0002645350890000012
Figure FDA0002645350890000013
wherein the content of the first and second substances,
Figure FDA0002645350890000014
is the inter-class variance at threshold k, mGIs the average gray level of the whole image, m (k) is the gray level when the threshold value is k, L is the number of gray levels of the image, P1(k) Is the probability that the pixel is classified into class 1,
a fourth step (S4) of normalizing the heterogeneous background, wherein the clustering center v is initializediAnd an offset field betajUpdating the membership function muijV cluster centeriAnd an offset field betajWhen | | | vb+1-vb||<When the cluster center calculated in the two previous and subsequent times is smaller than a preset convergence threshold value, stopping calculation to obtain background normalized data, wherein i is a cluster, and j is a pixel and is the preset convergence threshold value; b and b +1 denote the number of iterations, vb+1Is the clustering center of iteration b +1 times, vbIs the cluster center for iteration b times,
a fifth step (S5) of multi-scale Gaussian matched filtering enhancing the tiny gullies with Gaussian shape, filter rotation obtaining the response of the gullies with different directions, for each pixel, only keeping the maximum response value of the multiple directions, multiplying the response results of the matched filters with multiple scales after normalization,
a sixth step (S6) of adaptively thresholding the fine tidal trench, segmenting the fine tidal trench based on the global mean and the standard deviation, the thresholding formula: t ═ mean + k × std, where mean is the global mean; std is the global standard deviation; t is an optimal threshold value, the range of the parameter k is between 0.01 and 1, pixels with the value larger than T are divided into fine grooves, coarse fine groove segmentation results are obtained, finally isolated patches with the area smaller than the set threshold value are removed to obtain final fine groove segmentation results,
and a seventh step (S7) of merging the final segmentation result of the fine tide ditch of the sixth step (S6) with the logic OR operation of the wide tide ditch of the third step (S3) to form a complete tide ditch.
2. The tidal canal extraction method based on the high-resolution remote sensing image according to claim 1, wherein in the first step (S1), the multispectral data is collected by a multispectral camera, the multispectral data comprises a blue light band with a wavelength of 0.45-0.52 micrometers, a green light band with a wavelength of 0.52-0.59 micrometers, a red light band with a wavelength of 0.63-0.69 micrometers and a near infrared band with a wavelength of 0.77-0.89 micrometers, and the spatial resolution of the image is 4 m.
3. The tidal canal extraction method based on high-resolution remote sensing images as claimed in claim 1, wherein in the first step (S1), a FLAASH atmospheric correction algorithm based on a MODTRAN radiation transmission model is selected to perform atmospheric correction on the images.
4. The tidal channel extraction method based on the high-resolution remote sensing image according to claim 1, wherein in the second step (S2), the predetermined number of pixels is 5, a SEaTH algorithm is adopted, the separability between classes is calculated by using a J-M distance, the band is selected for fine tidal channel extraction, the GREEN band is selected for a mud beach area, and the NDWI band is selected for a salt marsh area, wherein NDWI is (GREEN-NIR)/(GREEN + NIR), GREEN refers to the GREEN band, and NIR refers to the near infrared band.
5. The tidal channel extraction method based on high-resolution remote sensing images according to claim 1, wherein in the fourth step (S4),
observed value r at pixel jjIs determined by the true value xjAnd an offset field betajObtained by addition, by subtracting an offset field betajWe obtain the true reflectivity at pixel j, as follows:
rj=xjj,j∈[1,N]
the global objective function J is such that,
Figure FDA0002645350890000021
wherein C represents the number of cluster centers; n represents the sum of the number of pixels in the image; alpha controls the strength of the neighborhood effect;
Figure FDA0002645350890000022
representing the average value of the pixels in the neighborhood window of the pixel j; mu.sijA membership function representing that the pixel j belongs to the class i; m is a membership factor set to 2; v. ofiA cluster center representing class i; j is a global objective function, and when the global objective function reaches an extreme value, background homogenization data are obtained.
6. The tidal channel extraction method based on high-resolution remote sensing images according to claim 1, wherein in the fifth step (S5), the rotation interval of the filter is set to 10 degrees, the filter is rotated 18 times to cover the tidal channels in all possible directions, the length of the x-axis of the matched filter is set to | x | -3 σ, the width of the y-axis of the filter is set to 5 to detect the tidal channel with the length of 5 pixels, the width of the tidal channel to be enhanced is between 1-5 pixels,
Figure FDA0002645350890000031
w is the tide to be enhancedOne half of the trench width, σ, is the standard deviation of the second derivative of the gaussian function, representing the distribution of intensity.
7. The tidal channel extraction method based on high-resolution remote sensing images as claimed in claim 1, wherein in the fourth step (S4), after the cluster center performs selection screening according to the original image histogram, the cluster center is selected at the peaks and valleys of the histogram by combining with the selected sample points.
8. The tidal channel extraction method based on high-resolution remote sensing images according to claim 1, wherein in the sixth step (S6), the parameter k is 0.2 in adaptive threshold segmentation.
CN201910059401.7A 2019-01-22 2019-01-22 Tidal ditch extraction method based on high resolution remote sensing image Active CN109801306B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910059401.7A CN109801306B (en) 2019-01-22 2019-01-22 Tidal ditch extraction method based on high resolution remote sensing image

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910059401.7A CN109801306B (en) 2019-01-22 2019-01-22 Tidal ditch extraction method based on high resolution remote sensing image

Publications (2)

Publication Number Publication Date
CN109801306A CN109801306A (en) 2019-05-24
CN109801306B true CN109801306B (en) 2020-10-27

Family

ID=66559980

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910059401.7A Active CN109801306B (en) 2019-01-22 2019-01-22 Tidal ditch extraction method based on high resolution remote sensing image

Country Status (1)

Country Link
CN (1) CN109801306B (en)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110689549B (en) * 2019-09-30 2020-12-18 北京航天宏图信息技术股份有限公司 Object extraction method, device and equipment
CN110991019B (en) * 2019-11-22 2023-05-05 首都师范大学 Automatic grading algorithm for coastal tidal ditches
CN117346744B (en) * 2023-12-04 2024-03-19 山东科技大学 Method for inverting measured water depth based on satellite-borne active and passive remote sensing information during rising and falling tide
CN118072192B (en) * 2024-04-25 2024-06-21 烟台哈尔滨工程大学研究院 GEE-based large-sized algae harmful algal bloom long-time sequence automatic extraction method

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101634705B (en) * 2009-08-19 2011-12-07 西安电子科技大学 Method for detecting target changes of SAR images based on direction information measure
KR101924278B1 (en) * 2017-02-16 2019-02-22 목포대학교산학협력단 Apparatus and method for extraction of tidal creeks in tidal flat zones using drone

Also Published As

Publication number Publication date
CN109801306A (en) 2019-05-24

Similar Documents

Publication Publication Date Title
CN109801306B (en) Tidal ditch extraction method based on high resolution remote sensing image
CN109190538B (en) Sediment-laden river delta coastal zone evolution analysis method based on remote sensing technology
Witharana et al. Understanding the synergies of deep learning and data fusion of multispectral and panchromatic high resolution commercial satellite imagery for automated ice-wedge polygon detection
CN111881816B (en) Long-time-sequence river and lake ridge culture area monitoring method
CN110991393A (en) Method and device for remote sensing monitoring and analysis of coastline transition
CN111027446B (en) Coastline automatic extraction method of high-resolution image
Li et al. Accurate water extraction using remote sensing imagery based on normalized difference water index and unsupervised deep learning
CN111402169B (en) Method for repairing remote sensing vegetation index time sequence under influence of coastal tide
CN114724049A (en) Inland culture pond water surface identification method based on high-resolution remote sensing image data
CN111310681B (en) Mangrove forest distribution remote sensing extraction method integrated with geoscience knowledge
Jin et al. River body extraction from sentinel-2A/B MSI images based on an adaptive multi-scale region growth method
CN109409265B (en) Floating raft culture area extraction method based on land resource satellite images
CN111597930A (en) Coastline extraction method based on remote sensing cloud platform
Kaiser Environmental changes, remote sensing, and infrastructure development: The case of Egypt's East Port Said harbour
CN108764132A (en) A kind of lake and marshland remote sensing images error detection method
Gong et al. Extracting tidal creek features in a heterogeneous background using Sentinel-2 imagery: a case study in the Yellow River Delta, China
CN112509134A (en) Tidal flat digital elevation model construction method and system
Zhao et al. Stability evaluation of tidal flats based on time-series satellite images: A case study of the Jiangsu central coast, China
CN112037244A (en) Landsat-8 image culture pond extraction method combining index and contour indicator SLIC
CN114519824A (en) Rapid detection method for SAR image flood inundation area
Teodoro et al. Extraction of Cabedelo sand spit area (Douro estuary) from satellite images through image processing techniques
Yu et al. Coastline detection using optical and synthetic aperture radar images
Wang et al. Deriving natural coastlines using multiple satellite remote sensing images
CN104596448A (en) Eutrophic water aquatic vegetation remote sensing extraction method based on alga index frequency method
CN110033460B (en) Satellite image offshore culture area extraction method based on scale space transformation

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant