A method for rapid change detection in multi-source satellite remote sensing images

Abstract

The present invention provides a method for rapid change detection in multi-source satellite remote sensing images. By using a pixel-oriented projection transformation method, relative radiation correction and change area pre-detection of multi-source sensors are simultaneously achieved. Areas that may have changed are preliminarily screened and only the areas that may have changed are processed, avoiding global operations for multi-temporal image classification and segmentation, which can greatly reduce the amount of image segmentation calculations. Four types of ground objects, namely buildings, water bodies, vegetation, and bare land, are then identified and segmented in local areas that may have changed. This avoids the problem that pixel-level change detection methods are affected by illumination and noise, while also avoiding the problem that object segmentation errors caused by a small target range on the global image, which leads to inaccurate change detection.

Description

Technical Field

[0001] The invention belongs to the technical field of rapid remote sensing image processing and computer vision, and in particular relates to a method for detecting rapid changes in remote sensing images. Background Art

[0002] Change detection is the process of detecting changes on the Earth's surface by comparing two or more remote sensing images acquired over the same geographic area at different points in time. It has been widely used in areas such as land use change analysis, disaster assessment, and resource and environmental monitoring. Traditional change detection methods are sensitive to differences in radiometric intensity and noise in multi-temporal remote sensing imagery. Differences in shooting angles and times can cause geometric errors and radiometric differences between remote sensing images captured by different sensors at different times, affecting the accuracy of subsequent change detection.

[0003] Existing change detection methods are mainly divided into pixel-based and object-based methods. Pixel-based change detection methods are computationally simple and fast, but the resulting results are fragmented and prone to the appearance of false change information similar to noise. Object-based detection methods, on the other hand, require separate classification and segmentation of two remote sensing images, processing large amounts of data, and their accuracy is affected by the segmentation algorithm. Summary of the Invention

[0004] The purpose of the present invention is to provide a method for rapid change detection in multi-source satellite remote sensing images. The IR-MAD transformation is used to simultaneously perform relative radiation correction of multi-source sensors and pre-detection of change areas on multi-temporal images. The trained YOLO8-seg segmentation network is then used to segment objects in areas where changes may occur. Finally, a differential calculation is performed to obtain change detection results. This method can reduce the amount of computational complexity in multi-temporal image classification and segmentation while improving its accuracy and efficiency.

[0005] The present invention is implemented by the following technical solutions:

[0006] A method for detecting rapid changes in satellite remote sensing images comprises the following steps:

[0007] Step 1: Geometric registration: Taking the first moment image Img1 as the reference image, a feature extraction algorithm is used to extract geometric features from the first moment image Img1 and the second moment image Img2, and feature matching is performed. The geometric mapping relationship between the second moment image Img2 and the first moment image Img1 is obtained by establishing an affine transformation model, and the second moment image IMg2 is geometrically transformed according to the geometric mapping relationship to obtain the registration image IMg'2 at the second moment compared to the first moment, so that the two images are in the same spatial coordinate system.

[0008] Step 2, IR-MAD transformation: Perform an iterative weighted multivariate detection algorithm IR-MAD transformation on the first moment image Img1 and the registered image Img'2 obtained in step 1 to find the unchanged pixel sample PIF and the changed pixel sample CIF. The unchanged pixel sample PIF and the changed pixel sample CIF record the position information of the unchanged pixel in the image and the position information of the changed pixel in the image respectively;

[0009] Step 3: Relative radiation correction: Based on the unchanged pixel sample set PIF obtained in step 2, a relative radiation correction model is constructed, and the relative radiation correction model is solved to obtain radiation correction model parameters. The relative radiation correction model parameters are used to perform radiation correction on the registration image Img'2 obtained in step 1 to obtain the corrected image Img at the second moment. 2_new ;

[0010] Step 4: Target segmentation: For the first moment image Img1 and the second moment corrected image Img obtained in step 3 2_new , set the slice window size Size, slice into blocks, and obtain the slice sequence Tile1 of the first moment image Img1 and the corrected image Img at the second moment 2_new Slice sequence Tile2; traverse the slices, when the first moment image Img1 and the second moment corrected image Img 2_new When slices C1 and C2 at the same position contain the changed pixel sample CIF obtained in step 2, target segmentation is performed on slices C1 and C2 respectively, otherwise it is skipped; the trained target segmentation network model is used to segment slices C1 and C2, and the target type class1 and target shape vector shape1 in slice C1, and the target type class2 and target shape vector shape2 in slice C2 are extracted, and stored in the segmentation results Tar1 and Tar2 respectively;

[0011] Step 5: Differential detection: perform differential calculation on the target shape vectors shape1 and shape2 in the segmentation results Tar1 and Tar2 obtained in step 4, and perform threshold segmentation to obtain the changed target shape vector shape change , the target type class1 corresponding to Tar1 and the class2 corresponding to Tar2 are respectively used as the target type class before and after the change before ,class after , store the change result C = [shape change ,class before ,class after ].

[0012] Furthermore, step one specifically includes: taking the first image Img1 as the reference image, using the feature extraction algorithm SIFT operator to extract geometric features of the first image Img1 and the second moment image Img2 respectively, and generating a 128-dimensional SIFT feature descriptor; setting the neighborhood radius to calculate the distance between each pair of key points of the SIFT feature descriptor, and returning the point corresponding to the minimum distance as the best matching point; according to the coordinates (x1, y1) of the best matching point on the first moment image Img1 and the coordinates (x2, y2) of the second moment image Img2, substituting them into the affine transformation model to obtain the affine transformation parameters; according to the affine transformation parameters, performing a geometric transformation on the second moment image Img2 to obtain the registered image Img'2 of the second moment compared to the first moment, so that the two images are in the same spatial coordinate system.

[0013] Furthermore, step 2 specifically includes: performing linear transformations U and V respectively according to the spectral vector values ​​of the two images Img1 and Img'2, and calculating the MAD variable according to the difference between the two transformation vectors; calculating the invariant probability P of each pixel according to the MAD no_change As the iteration weight, P no_change Repeat the iteration as the next iteration weight until the two typical correlation coefficients ρ N-i+1 The difference is less than the threshold λ or the number of iterations reaches the threshold T. no_change The set of pixel positions greater than the threshold p is used as PIF, and the remaining pixels are used as changed pixel samples CIF.

[0014] Furthermore, step three specifically includes: according to the unchanged pixel sample set PIF obtained in step two, obtaining the pixel reflectance set I1 of the first moment image Img1 corresponding to PIF and the pixel reflectance set I2 of the registration image Img'2 obtained in step one, constructing a relative radiation correction model based on I1 and I2, and using the least squares method to calculate the radiation correction model parameters; using the least squares method to solve the linear regression parameters; according to the solved linear regression parameters, performing radiation correction on the registration image Img'2 obtained in step one to obtain the corrected image Img at the second moment 2_new .

[0015] Furthermore, in step 4, the slice window size is 512 or 1024.

[0016] Furthermore, the target segmentation network model trained in step 4 is a YOLOv8m-seg network model.

[0017] Furthermore, the target shape vector in step 4 is used to record the XY set of vertices of each target polygon.

[0018] Furthermore, the target types in steps 4 and 5 include buildings, water bodies, vegetation, and bare land.

[0019] The present invention uses a pixel-oriented projection transformation method to simultaneously achieve relative radiation correction of multi-source sensors and pre-detection of changed areas, conduct preliminary screening of areas where changes may occur, and only process areas where changes may occur, avoiding global operations for multi-temporal image classification and segmentation, which can greatly reduce the amount of image segmentation calculations; then, four types of ground objects, namely buildings, water bodies, vegetation, and bare land, are identified and segmented in local areas where changes may occur. On the one hand, this avoids the problem that pixel-level change detection methods are affected by light and noise, and at the same time, it also avoids the problem that object segmentation errors caused by too small a target range on the global image, resulting in inaccurate change detection. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] Figure 1 This is a flow chart of a method for detecting rapid changes in remote sensing images according to one embodiment of the present invention;

[0021] Figure 2 It is a schematic diagram of the detection results of one embodiment of the present invention. DETAILED DESCRIPTION

[0022] 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 in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts shall fall within the scope of protection of the present invention.

[0023] See also Figure 1 The embodiment of the present invention provides a method for detecting rapid changes in remote sensing images, comprising the following steps:

[0024] Step 1: Geometric registration. Using the first moment image Img1 as the reference image, a feature extraction algorithm is used to extract geometric features from the first moment image Img1 and the second moment image Img2, respectively, and feature matching is performed. An affine transformation model is established to obtain the geometric mapping relationship between the second moment image Img2 and the first moment image Img1. Based on this relationship, a geometric transformation is performed on the second moment image Img2 to obtain the registered image Img'2 at the second moment compared to the first moment Img1, so that the two images are in the same spatial coordinate system.

[0025] The step 1 specifically includes:

[0026] Using the first image Img1 as the reference image, the SIFT operator is used to extract geometric features from both the first image Img1 and the second image Img2, generating 128-dimensional SIFT feature descriptors. A neighborhood radius is set to calculate the distance between each pair of key points in the SIFT feature descriptors, and the point with the minimum distance is returned as the best matching point. The coordinates (x1, y1) of the best matching point on the first image Img1 and the coordinates (x2, y2) of the second image Img2 are substituted into the affine transformation model to obtain the affine transformation parameters. Based on the affine transformation parameters, a geometric transformation is performed on the second image Img2 to obtain the registered image Img'2 at the second moment relative to the first image Img1, aligning the two images to the same spatial coordinate system.

[0027] Step 2: IR-MAD transformation. Perform an iterative weighted multivariate detection algorithm IR-MAD transformation on the first moment image Img1 and the registered image Img'2 obtained in step 1 to find the unchanged pixel sample PIF and the changed pixel sample CIF. The unchanged pixel sample PIF and the changed pixel sample CIF respectively record the position information of the unchanged pixel in the image and the position information of the changed pixel in the image.

[0028] The second step specifically includes:

[0029] According to the spectral vector values ​​of the two images Img1 and Img'2, linear transformations U and V are performed respectively, and the MAD variable is calculated based on the difference between the two transformation vectors; the invariant probability P of each pixel is calculated based on the MAD no_change As the iteration weight, P no_change Repeat the iteration as the next iteration weight until the two typical correlation coefficients ρ N-i+1 The difference is less than the threshold λ or the number of iterations reaches the threshold T. no_change The set of pixel positions greater than the threshold p is used as PIF, and the remaining pixels are used as changed pixel samples CIF.

[0030] Step 3: Relative radiation correction: Based on the unchanged pixel sample set PIF obtained in step 2, a relative radiation correction model is constructed, and the relative radiation correction model is solved to obtain radiation correction model parameters. The relative radiation correction model parameters are used to perform radiation correction on the registration image Img'2 obtained in step 1 to obtain the corrected image Img at the second moment. 2_new ;

[0031] The step three specifically includes:

[0032] Based on the unchanged pixel sample set PIF obtained in step 2, the pixel reflectance set I1 of the first moment image Img1 corresponding to PIF and the pixel reflectance set I2 of the registration image Img'2 obtained in step 1 are obtained. Based on I1 and I2, a relative radiation correction model is constructed, and the radiation correction model parameters are calculated using the least squares method. The linear regression parameters are solved using the least squares method; based on the solved linear regression parameters, the registration image Img'2 obtained in step 1 is subjected to radiation correction to obtain the corrected image Img at the second moment. 2_new ;

[0033] Step 4: Target segmentation. For the first moment image Img1 and the second moment corrected image Img obtained in step 3 2_new , manually set the slice window size Size (Size is usually 512 or 1024), and obtain the slice sequence Tile1 of the first moment image Img1 and the corrected image Img at the second moment 2_new Slice sequence Tile2; traverse all slices, when Img1 and Img 2_new When the slices C1 and C2 at the same position contain the changed pixel samples CIF obtained in step 2, the slice is segmented, otherwise the slice is skipped. Based on the pre-trained YOLOv8m-seg network model target segmentation network, the pixel sets C1 and C2 are segmented respectively, and the target type class1 (target types include buildings, water bodies, vegetation, bare land) and target shape vector shape1 (recording the XY set of each target polygon vertex) in C1, and the target type class2 and target shape vector shape2 (recording the XY set of each target polygon vertex) in C2 are extracted, and stored in the target sets Tar1 and Tar2 respectively, that is, [class1, shape1]∈Tar1, [class2, shape2]∈Tar2;

[0034] Step 5: Differential detection: perform differential calculation on the target shape vector shape2 and shape1 in the segmentation results Tar1 and Tar2 obtained in step 4, and perform threshold segmentation to obtain the changed target shape vector shape change , the target type class1 corresponding to Tar1 and the class2 corresponding to Tar2 are respectively used as the target type class before and after the change before ,class after , store the change result c = [shape change ,class before ,class after ].

[0035] The following is a specific example to illustrate:

[0036] Step 1: Geometric registration: Perform geometric registration on the two images Img1 and Img2 at different times to obtain the registered images Img1 and Img'2.

[0037] (1) Taking the first moment image Img1 as the reference image, extract the SIFT features of Img1 and Img2, and generate 128-dimensional SIFT feature descriptors S1 and S2.

[0038] (2) Feature matching. Set the neighborhood radius R, calculate the distance between the 128-dimensional SIFT feature descriptors S1 and S2 between each pair of key points within R, and return the minimum distance among the k best matches for each key point (k = 5 in this example). The point with the minimum distance is the matching point of the key point.

[0039] (3) Geometric transformation: Based on the feature matching results, the least squares method is used to solve the affine transformation parameters.

[0040] Assume that the matching point of feature point Y on Img1 on Img2 is X, then:

[0041] x1=ax2+cy2+m (1)

[0042] y1=bx2+dy2+l (2)

[0043] Right now

[0044]

[0045] Among them, a, b, c, d, m, and l are affine transformation parameters.

[0046] make: Then Y=AX, find the affine transformation parameter A, then

[0047] A=(X T X) -1 X T Y (4)

[0048] (4) Use formula 3 to perform affine transformation on Img2 to obtain Img'2.

[0049] Step 2: IR-MAD transformation. Perform an iterative weighted multivariate detection algorithm IR-MAD transformation on the first moment image Img1 and the registered image Img'2 obtained in step 1 to find the unchanged pixel sample PIF and the changed pixel sample CIF. The unchanged pixel sample PIF and the changed pixel sample CIF respectively record the position information of the unchanged pixel in the image and the position information of the changed pixel in the image.

[0050] (1) Set the initial weight of all pixels to 1,

[0051] (2) Let variable X = [X1…X N ] T ,Y=[Y1…Y N ] T Corresponding to the spectral vector values ​​of the two images Img1 and Img'2, N is the number of bands. Perform linear transformation on X and Y respectively:

[0052] U=a T X=a1X1+…a i X i +a N X N

[0053] V=b T X=b1Y1+…b i Y i +b N Y N

[0054] Let the MAD variable be the difference between U and V, and obtain the linear transformation parameters of the maximum correlation coefficient between vectors U and V through CCA:

[0055] MAD i =U N-i+1 -V N-i+1 , i=1…N

[0056] make MAD i The standard deviation of is the variance, ρ N-i+1 is the typical correlation coefficient, then

[0057]

[0058] but:

[0059]

[0060] N MAD is the normalized MAD variable, and the MAD mean for unchanged pixels is approximately 0.

[0061] (3) Calculate the unchanged probability of each pixel through the chi-square test:

[0062] P no_change =P{χ 2 (N)>N MAC}

[0063] (4)P no_change As the weight for the next iteration;

[0064] (5) Repeat steps (2) to (4) until ρ is twice N-i+1The difference is less than the threshold λ or the number of iterations reaches the threshold T.

[0065] (6)P no_change The set of pixel positions greater than the threshold p is used as PIF, and the remaining pixels are used as changed pixel samples CIF.

[0066] Step 3: Relative radiation correction. Based on the unchanged pixel sample set PIF obtained in step 2, a relative radiation correction model is constructed, and the relative radiation correction model is solved to obtain the radiation correction model parameters. The relative radiation correction model parameters are used to perform radiation correction on the registration image Img'2 obtained in step 1 to obtain the corrected image Img at the second moment. 2_new .

[0067] (1) Use the least squares method to find the linear regression parameters of the radiation correction model. Let the spectral values ​​of the pixels Img1 and Img'2 at the corresponding positions in PIF be PX i ,PY i ,but:

[0068] PX=A·PY+B

[0069] A and B are linear regression parameters. Least squares method is used to solve for parameters A and B.

[0070] (2) According to the solved A and B values, perform projection transformation on Img'2:

[0071] Img 2_new =A·Img′2+B

[0072] Step 4: Target segmentation. For Img1 and Img 2_new ∈CIF, perform target segmentation based on the target segmentation network, extract the target type within the window W, and obtain Tar1 and Tar2.

[0073] (1) Manually set the slice window size Size = 512, and compare the first moment image Img1 and the second moment corrected image Img obtained in step 3 2_new The whole image is sliced ​​into blocks of size Size, and the slice sequence Tile1 of the first moment image Img1 and the corrected image Img at the second moment are obtained. 2_new Slice sequence Tile2;

[0074] (2) When the image Img1 at the first moment and the image Img after correction at the second moment 2_newIf the slices C1 and C2 at the same position on the image contain the changed pixel samples CIF obtained in step 2, target segmentation is performed on C1 and C2 respectively; otherwise, the process is skipped. The trained target segmentation network model is used to segment C1 and C2. The target type class1 (target types include buildings, water bodies, vegetation, and bare land) and target shape vector shape1 (recording the XY set of each target polygon vertex) in C1, and the target type class2 and target shape vector shape2 in C2 are extracted and stored in the segmentation results Tar1 and Tar2, respectively.

[0075] Step 5: Differential detection: Perform differential calculation on Tar1 and Tar2 obtained in step 4 to obtain the change result C.

[0076] (1) For the position (x, y) ∈ CIF, obtain the segmentation results of Tar1 and Tar2, which are [class1, shape1] and [class2, shape2] respectively, and determine whether it is a change

[0077]

[0078] (2) Mark the changed Tar2 and obtain the change-unchange binary classification result. change , the target type class1 and class2 of the target are respectively used as the target type class before and after the change before ,class after , store the change result c = [shape change ,class before ,class after ].

[0079] Figure 2 This is the detection result of one embodiment of the present invention. Satellite images of a reservoir area from 2017 and 2022 were registered, transformed using IR-MAD, and corrected for relative radiometry to identify candidate change regions (CIFs). The 2017 and registered and corrected 2022 images were then divided into blocks, and blocks containing CIFs of changed pixels were traversed. Figure 2 The blue window shown in (a) is skipped as the block without the CIF sample without the changed pixel, while the red window is one of the blocks overlapping with the CIF sample with the changed pixel, and the target recognition segmentation is performed on it. The segmentation results are shown in Figure 2. Figure 2 (b) Figure 2 The segmentation results of (b) are differentially calculated to obtain the change detection results, such as Figure 2 (c) shown.

[0080] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to the above embodiments, ordinary technicians in the field should understand that the specific implementation methods of the present invention can still be modified or replaced by equivalents. Any modification or equivalent replacement that does not depart from the spirit and scope of the present invention should be covered by the scope of protection of the claims of the present invention.

Claims

1. A method for detecting rapid changes in multi-source satellite remote sensing images, characterized in that: The steps include: Step 1: Geometric registration: Using the first moment image Img1 as the reference image, a feature extraction algorithm is used to extract geometric features from the first moment image Img1 and the second moment image Img2, and feature matching is performed. The geometric mapping relationship between the second moment image Img2 and the first moment image Img1 is obtained by establishing an affine transformation model. Based on the geometric mapping relationship, a geometric transformation is performed on the second moment image Img2 to obtain the registered image Img'2 at the second moment compared to the first moment, so that the two images are in the same spatial coordinate system. Step 2, IR-MAD transformation: Perform an iterative weighted multivariate detection algorithm IR-MAD transformation on the first moment image Img1 and the registered image Img'2 obtained in step 1 to find the unchanged pixel sample PIF and the changed pixel sample CIF. The unchanged pixel sample PIF and the changed pixel sample CIF record the position information of the unchanged pixel in the image and the position information of the changed pixel in the image respectively; Step 3: Relative radiation correction: Based on the unchanged pixel sample set PIF obtained in step 2, a relative radiation correction model is constructed, and the relative radiation correction model is solved to obtain radiation correction model parameters. The relative radiation correction model parameters are used to perform radiation correction on the registration image Img'2 obtained in step 1 to obtain the corrected image Img at the second moment. 2_new ; Step 4: Target segmentation: For the first moment image Img1 and the second moment corrected image Img obtained in step 3 2_new , set the slice window size Size, slice into blocks, and obtain the slice sequence Tile1 of the first moment image Img1 and the corrected image Img at the second moment 2_new Slice sequence Tile2; traverse the slices, when the first moment image Img1 and the second moment corrected image Img 2_new When slices C1 and C2 at the same position contain the changed pixel sample CIF obtained in step 2, target segmentation is performed on slices C1 and C2 respectively, otherwise it is skipped; the trained target segmentation network model is used to segment slices C1 and C2, and the target type class1 and target shape vector shape1 in slice C1, and the target type class2 and target shape vector shape2 in slice C2 are extracted, and stored in the segmentation results Tar1 and Tar2 respectively; Step 5: Differential detection: perform differential calculation on the target shape vectors shape1 and shape2 in the segmentation results Tar1 and Tar2 obtained in step 4, and perform threshold segmentation to obtain the changed target shape vector shape change , the target type class1 corresponding to Tar1 and the class2 corresponding to Tar2 are respectively used as the target type class before and after the change before ,class after , store the change result C = [shape change ,class before ,class after ].

2. The method for detecting rapid changes in multi-source satellite remote sensing images according to claim 1, wherein: Step one specifically includes: taking the first image Img1 as the reference image, using the feature extraction algorithm SIFT operator to extract geometric features from the first image Img1 and the second moment image Img2 respectively, and generating a 128-dimensional SIFT feature descriptor; setting the neighborhood radius to calculate the distance between each pair of key points of the SIFT feature descriptor, and returning the point corresponding to the minimum distance as the best matching point; according to the coordinates (x1, y1) of the best matching point on the first moment image Img1 and the coordinates (x2, y2) of the second moment image Img2, substitute them into the affine transformation model to obtain the affine transformation parameters; according to the affine transformation parameters, perform a geometric transformation on the second moment image Img2 to obtain the registration image IMg'2 of the second moment compared with the first moment, so that the two images are in the same spatial coordinate system.

3. The method for detecting rapid changes in multi-source satellite remote sensing images according to claim 1, wherein: Step 2 specifically includes: performing linear transformations U and V based on the spectral vector values ​​of the two images Img1 and Img'2, and calculating the MAD variable based on the difference between the two transformation vectors; calculating the invariant probability P of each pixel based on the MAD no_change As the iteration weight, P no_change Repeat the iteration as the next iteration weight until the two typical correlation coefficients ρ N-i+1 The difference is less than the threshold λ or the number of iterations reaches the threshold T; no_change The set of pixel positions greater than the threshold p is used as PIF, and the remaining pixels are used as changed pixel samples CIF.

4. The method for detecting rapid changes in multi-source satellite remote sensing images according to claim 1, wherein: According to the unchanged pixel sample set PIF obtained in step 2, the pixel reflectance set I1 of the first moment image Img1 corresponding to PIF and the pixel reflectance set I2 of the registration image Img'2 obtained in step 1 are obtained. A relative radiation correction model is constructed based on I1 and I2, and the radiation correction model parameters are calculated using the least squares method; the linear regression parameters are solved using the least squares method; according to the solved linear regression parameters, the registration image Img'2 obtained in step 1 is subjected to radiation correction to obtain the corrected image Img at the second moment 2_new .

5. The method for detecting rapid changes in multi-source satellite remote sensing images according to claim 1, wherein: In step 4, the slice window size is 512 or 1024.

6. The method for detecting rapid changes in multi-source satellite remote sensing images according to claim 1, wherein: The target segmentation network model trained in step 4 is the YOLOv8m-seg network model.

7. The method for detecting rapid changes in multi-source satellite remote sensing images according to claim 1, wherein: The target shape vector in step 4 is used to record the XY set of vertices of each target polygon.

8. The method for detecting rapid changes in multi-source satellite remote sensing images according to claim 1, wherein: Target types in steps 4 and 5 include buildings, water bodies, vegetation, and bare land.

Images