Valley topography extraction method and water system extraction method and program

Abstract

To provide a valley landform extraction method, a water system extraction method, and a program capable of objectively creating a valley landform and a water system without depending on a scale of a map or subjectivity of a creator.SOLUTION: A valley landform extraction method and a water system extraction method comprise steps of: calculating an average elevation for n×n cells included in a digital elevation model (hereafter DEM(Lmin)); creating a low-resolution digital elevation model (L) that assigns the average altitude as an elevation of a square cell with a side length L (Lmin×n); performing a step of calculating plane curvature of a center cell for each 3×3 cell group contained in the DEM(L); calculating a normal plane curvature, which is the plane curvature multiplied by a cell side length L; determining a cell with a normal planar curvature equal to or less than a first threshold as a depression area; superimposing the depression area extracted for each of a plurality of DEMs(L); creating a polygon consisting of contiguous depression areas; and extracting a polygon whose length in a direction of a valley line is greater than or equal to a second threshold as a valley area.SELECTED DRAWING: Figure 1

Description

[Technical field]

[0001] The present invention relates to a method and program for extracting valley topography and a water system from a digital elevation model (DEM). [Background technology]

[0002] Conventionally, it has been common to obtain mountain topographical representations and to extract valley topography and water system maps from digital elevation models in order to assess the risk of slope disasters.

[0003] A river system map is created by visually reading the curvature of contour lines on a map and drawing a continuous valley line from the most upstream point that is judged to be a valley to the downstream. The standard for judging a valley is generally a contour line with a folding angle of 90 degrees or 53 degrees.

[0004] In this method, the results of valley extraction depend on the creator, so different valley topography and river system maps are obtained depending on the creator. Also, creating a wide-area river system map takes a considerable amount of time.

[0005] In addition to this method, there is also a method to calculate the upstream catchment area using a digital elevation model, and obtain a valley topography and a water system map for catchment areas above a certain threshold. This threshold is determined arbitrarily by the creator.

[0006] In this method, even if a valley is not deep, if it has a large water volume, it will be considered a river system, so strictly speaking, the river system map does not follow the criteria of valley topography. In addition, horizontal surfaces created by artificial structures are also included in the river system map.

[0007] Furthermore, both of the above mentioned two methods have the problem that the resulting river system charts differ depending on the scale of the map.

[0008] The technology described in Patent Document 1 extracts slopes and other features considered to be water-collecting topography from a digital elevation model, but it does not solve the problems described above. [Prior art documents] [Patent documents]

[0009] [Patent Document 1] JP 2013-221886 A [Non-patent literature]

[0010] [Non-Patent Document 1] Kawata, Y. and Niki, N. 2001. Extraction algorithm of 3D curvature features. MEDICAL IMAGING TECHNOLOGY 19: 142-153 Summary of the Invention [Problem to be solved by the invention]

[0011] In view of the above circumstances, an object of the present invention is to provide a valley topography extraction method, a river system extraction method, and a program that can create valley topography and river systems objectively, regardless of the scale of the map or the creator's subjectivity. [Means for solving the problem]

[0012] To achieve the above object, an aspect of the present invention is a method for extracting valley topography using a digital elevation model (Lmin) in which the earth's surface is divided into square cells with one side length Lmin and an elevation is assigned to each cell, the method comprising: an aggregation step of calculating an average elevation for n×n cells included in the digital elevation model (Lmin) and creating a low-resolution digital elevation model (L) by assigning said average elevation as the elevation of a square cell with one side length L (Lmin×n); a plan curvature calculation step of extracting 3×3 cells (hereinafter, cell groups) included in the digital elevation model (L) and calculating the plan curvature of the center cell for all extracted cell groups, for each of the multiple digital elevation models (L) created in the aggregation step; The method for extracting a valley topography comprises a normal plane curvature calculation step in which a process of calculating a normal plane curvature obtained by multiplying a surface curvature by the length L of a cell's side is performed for each of the plurality of digital elevation models (L) created in the aggregation step; a depression determination step in which a first threshold value for the normal plane curvature for determining that the central cell of the group of cells is a depression area, the first threshold value being common to the plurality of digital elevation models (L), and a cell having the normal plane curvature equal to or less than the first threshold value is determined to be a depression area; and a valley area extraction step in which the depression areas extracted for each of the plurality of digital elevation models (L) are overlaid to create a polygon consisting of consecutive depression areas, and a polygon whose valley line direction length is equal to or greater than a second threshold value is extracted as a valley area.

[0013] Another aspect of the present invention for achieving the above-mentioned object is a method for extracting a water system using a valley area obtained by the above-mentioned valley topography extraction method, comprising a water system extraction step of dividing the valley area into square analysis cells with one side measuring 1, calculating the cumulative flow rate by assigning a constant value to each analysis cell, and extracting a water system consisting of analysis cells whose cumulative flow rate exceeds a fourth threshold, wherein the fourth threshold is the cumulative flow rate at the head of the V-shaped valley topography.

[0014] Another aspect of the present invention for achieving the above object resides in a program for causing a computer to execute each step of the above valley topography extraction method.

[0015] Another aspect of the present invention for achieving the above object resides in a program for causing a computer to execute each step of the above aqueous extraction method. Effect of the Invention

[0016] According to the present invention, there is provided a valley topography extraction method and a river system extraction method and program that can objectively create valley topography and river systems regardless of the map scale or the creator's subjectivity. [Brief description of the drawings]

[0017] [Figure 1] FIG. 1 is a flow chart of a method for extracting a valley topography and a method for extracting a water system. [Diagram 2] FIG. 2 is a diagram showing a portion of DEM(1) and DEM(2). [Diagram 3] FIG. 2 is a diagram showing a portion of DEM(1) and DEM(2). [Figure 4] This shows the elevation of a 3x3 cell in a DEM. [Diagram 5] FIG. 1 is a diagram showing the relationship between the plan curvature and the normal plan curvature in a V-shaped valley topography and the valley topography. [Figure 6] FIG. 13 is a diagram for explaining aggregation of DEMs when n is an odd number. [Figure 7] FIG. 13 is a diagram for explaining aggregation of DEMs when n is an even number. [Figure 8] FIG. 13 is a diagram for explaining aggregation of DEMs when n is an even number. [Figure 9] FIG. 13 is a diagram for explaining aggregation of DEMs when n is an even number. [Figure 10] FIG. 13 illustrates an example of a depression region. [Figure 11] FIG. 1 illustrates valley regions and slopes. [Figure 12] FIG. 13 is a diagram for explaining a method for extracting a valley region based on length. [Figure 13] This is a diagram showing cumulative discharge in a V-shaped valley. [Figure 14] FIG. 1 is a diagram illustrating a method for removing short valleys from a water system. [Figure 15] FIG. 13 is a diagram for explaining a method for thinning a water system. [Figure 16] FIG. 13 is a diagram showing an example of the extraction method of valley topography and water systems from a DEM. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0018] A valley topography extraction method and a water system extraction method according to an embodiment of the present invention will be described. Fig. 1 is a flow diagram of the valley topography extraction method and the water system extraction method. First, in the valley topography extraction method, a digital elevation model with a lower resolution is created from a digital elevation model with the highest resolution that serves as a reference.

[0019] A digital elevation model is data in which the earth's surface is divided into square cells with one side of Lmin, and an elevation is assigned to each cell. In this embodiment, a case where Lmin is 1m will be described. Hereinafter, a digital elevation model with a cell length of L (meters) will be referred to as DEM(L).

[0020] Figure 2(a) shows an example of a portion of DEM(1). An elevation is assigned to each cell, and the elevation of each cell is hij (i = 1, 2, ; j = 1, 2, ). The average elevation is calculated for the n × n cells included in this DEM(1), and this average elevation is assigned as the elevation of a square cell with a side length of L (Lmin × n). A low-resolution digital elevation model (L) consisting of cells of length L obtained in this way is created (aggregation step; step S1 in Figure 1).

[0021] For example, Figure 2(b) shows DEM(2) where n is 2 and length L is 2. First, calculate the average elevation for the 2 x 2 cells contained in DEM(1). For example, let H11 be the average elevation of the elevations h11, h12, h21, and h22 of the 2 x 2 cells in the upper left.

[0022] Then, the above-mentioned average elevation H11 is assigned as the elevation of one cell with a side length of 2m, which is the same as that of the 2 x 2 cells. Similarly, the average elevation is calculated for the 2 x 2 cells included in DEM(1) and set as the elevation of the cell with a length of 2m, to create DEM(2).

[0023] In the example of FIG. 2, the case has been described where the boundaries (sides) of the cells of DEM(2) coincide with the boundaries of DEM(1), but the resolution is not limited to this type of reduction.

[0024] Figure 3(a) shows the state where the boundaries of the cells of DEM(2) (shown by dotted lines) are shifted by half a degree from the boundaries of the cells of DEM(1). Figure 3(b) shows an excerpt of DEM(2) from Figure 3(a).

[0025] DEM(2) is offset from each side of DEM(1) by 0.5 (half the length of the side of DEM(1)). When the boundaries of the low-resolution and high-resolution cells do not match, the average elevation is calculated from the elevations of multiple cells of DEM(1) that overlap one cell of DEM(2).

[0026] For example, the elevation H11 of the top left cell of DEM(2) can be calculated from the elevations h11, h12, h13, h21, h22, h23, h31, h32, and h33 of the nine cells in DEM(1) that overlap with that cell. The detailed calculation will be explained later with reference to Figure 9.

[0027] Although not shown, if the boundary shift is other than half the cell length of DEM(1), an appropriate weighted average can be calculated according to the position of the boundary shift between DEM(1) and DEM(2).

[0028] In this way, depending on the position of the boundary, multiple different DEMs (2) are obtained from DEM (1), as shown in Figures 2(b) and 3. In this embodiment, a DEM (2) whose boundary coincides with DEM (1), and a DEM (2) whose boundary is shifted by half the cell length of DEM (1) are created from DEM (1).

[0029] It is preferable to create multiple low-resolution DEMs, not just one. In the embodiment described below, DEM(2) to DEM(500) were created at intervals of 1 m based on DEM(1). In addition, two patterns were created for each DEM(L), one with the same boundary as described above, and one with the boundary shifted by half the cell length of DEM(1).

[0030] Next, the plan curvature is calculated for multiple DEMs (L) (plan curvature calculation step; step S2 in FIG. 1). If the plan curvature is Kplan, the cross-sectional curvature is Kcross, the maximum slope angle on the valley line is θ, and the elevations of the 3×3 cells (hereinafter also referred to as "cell group") with sides of length d shown in FIG. 4 are Z1 to Z9, then they can be found as follows. The calculation of Kplan and Kcross is a known method described in Non-Patent Document 1, and a detailed explanation will be omitted.

number

[0031] Both plan curvature and cross section curvature are one of the curvatures of a cross section line extending in a direction perpendicular to the slope direction. Plan curvature is the curvature perpendicular to the maximum slope angle term. Cross section curvature is the curvature in the direction of maximum slope angle when the cross section being examined is taken parallel to the elevation direction. Negative values ​​of plan curvature and cross section curvature indicate an upwardly concave topography, and positive values ​​indicate an upwardly convex topography.

[0032] In this way, in the plan curvature calculation step, cells contained in a DEM are extracted, and the plan curvature of the center cell of each extracted cell group is calculated using the above formula (1). This calculation is then performed for all of the multiple DEMs created in the aggregation step.

[0033] Next, the normal plane curvature is calculated (normal plane curvature calculation step; step S3 in Fig. 1). The normal plane curvature is obtained by multiplying the plane curvature by the length of the side of the cell. Specifically, the normal plane curvature is calculated by multiplying the plane curvature of the cell included in the DEM(L) by the length L of the side of the cell, and this calculation is performed for all the DEMs created in the aggregation step.

[0034] Here, the relationship between plan curvature, normal plan curvature, and valley topography will be explained using Figure 5. Figures 5(a) to (c) show a group of cells included in a DEM with a certain resolution and their elevations. Figures 5(d) to (f) show the elevation difference and height difference of Figures 5(a) to (c). Figure 5(g) shows a cross section in the X direction of Figures 5(a) to (c). Figure 5(h) shows a cross section in the Y direction at the valley line of Figures 5(a) to (c). Figure 5(i) shows a cross section in the X direction when Figure 5(a) is doubled in the horizontal and vertical directions.

[0035] The simplest V-shaped valley topography is one in which two slopes with the same inclination angle and height join at a valley line. This type of V-shaped valley topography is expressed by a specific pattern of elevation in the group of cells shown in Figures 5(a) to 5(c). In Figures 5(a) to 5(c), the dotted lines indicate contour lines, and the numbers and symbols in each cell indicate elevation. The cell in row i and column j is also written as Cij. In addition, the three cells (C12, C22, and C32) located in the center of the group of cells and aligned in the Y direction contain valley lines.

[0036] In order for the shape of all the contour lines within a cell to be uniform, firstly, the difference in elevation between the central cell, which forms the valley line, and its adjacent cells must be constant. Secondly, in order for the intervals between the contour lines to be uniform, the difference in elevation between the cells in the downstream direction must be constant.

[0037] To illustrate the first condition, in the pattern of Fig. 5(a), as shown in Fig. 5(d), the elevation difference between the central cell C12 and its neighboring cells C11 and C13 is 2a. The elevation difference between the central cells C22 and C23 and their neighboring cells is also 2a. In other words, the elevation difference is constant. Similarly, in the pattern of Fig. 5(b), as shown in Fig. 5(e), the elevation difference is constant at a, and in the pattern of Fig. 5(c), as shown in Fig. 5(f), the elevation difference is constant at 1 / 2a.

[0038] To illustrate the second condition, in the pattern of Fig. 5(a), as shown in Fig. 5(d), the elevation difference between the upstream cell C21 and the downstream cell C11 is a. The elevation difference between the other upstream and downstream cells is also a. In other words, the elevation difference is constant. Similarly, in the pattern of Fig. 5(b), as shown in Fig. 5(e), the elevation difference is constant at a, and in the pattern of Fig. 5(c), as shown in Fig. 5(f), the elevation difference is constant at a.

[0039] The V-shaped valley shapes shown in Figures 5(a) to 5(c) are determined by the elevation difference a in the valley line direction (Y direction) and the magnification k of the elevation difference in the direction perpendicular to the valley (X direction), as shown in Figures 5(g) and 5(h). Here, k is positive when the elevation of the slope decreases toward the valley line. In other words, k is positive in the case of a V-shaped valley shape.

[0040] The k of the V-shaped valley shape shown in FIG. 5(a) to FIG. 5(c) is 2, 1, and 1 / 2, respectively. The maximum inclination angle θ on the valley line is tan -1 (a / L), and the inclination angle Φ of the side wall is tan -1 (ka / L). In each of the V-shaped valley topographies in Figures 5(a) to 5(c), the curvature in the maximum dip direction in the central cell is 0, and the cross-sectional curvature Kcross is calculated from the elevation difference between the central cell and the cells on the left and right, and the cell width L. Therefore, the plan curvature in the patterns shown in Figures 5(a) to 5(c) can be expressed as equation (2).

number

[0041] For example, when the opening width and depth of the contour lines are equal (when the contour line folding angle is approximately 53 degrees) as in Figure 5(a), the plan curvature Kplan is -4 / L when k = 2. When the contour line folding angle is 90 degrees as in Figure 5(b), the plan curvature Kplan is -2 / L when k = 1. These plan curvatures have the same value even when the DEM is rotated on a plane.

[0042] As shown in the above formula (2), the plan curvature Kplan depends on the cell length L. For example, as shown in Figure 5(i), if the shape of the valley is expanded twice in the horizontal and vertical directions, the plan curvature will be halved. On the other hand, the value obtained by multiplying the plan curvature by the cell width remains constant. This value is called the normal plan curvature, and is also written as Kplan,normal.

[0043] The absolute value of the normal plan curvature does not depend on the width of the valley or the cell size, but is large when the center of the cell is in the center of the V-shaped valley and becomes smaller as the cell moves away from the valley line. By utilizing this property, it is possible to extract V-shaped valleys of different scales using a single threshold value by calculating the normal plan curvature at different cell lengths using a group of cells.

[0044] For example, in the V-shaped valley topography of the patterns in Figures 5(a) and 5(b), Kplan,normal is calculated for each different cell length L (different scales) as follows. L Kplan, normal 1 -4 -2 2 -4 -2 3 -4 -2 ...

[0045] If Kplan,normal for a cell is -4 or less, then the cell can be extracted as a V-shaped valley topography with a pattern like that of Figure 5(a) or a V-shaped valley topography with a smaller folding angle (strictly speaking, it is a depression because it is a single cell) regardless of the cell length or scale. Similarly, if Kplan,normal for a cell is -2 or less, then the cell can be extracted as a V-shaped valley topography with a pattern like that of Figure 5(b) or a V-shaped valley topography with a smaller folding angle regardless of the cell length or scale.

[0046] A depression is determined based on the normal plan curvature with such properties (depression determination step; step S4 in Figure 1). Specifically, a first threshold value for the normal surface curvature is appropriately set for determining that the central cell in the group of cells is a valley, as in "2" or "4" above. The first threshold value is common to multiple DEMs (L), as described above. Cells contained in a DEM whose normal plan curvature is equal to or less than the first threshold value are determined to be valleys (depressions). Then, using the same first threshold value, extraction of such valleys is performed for all DEMs (L).

[0047] However, when a low-resolution DEM is obtained from a reference DEM (DEM(1) in this embodiment), the normal plan curvature changes depending on the side length of the cell used in averaging (the length determined by n in the n×n cell used in averaging). This is because the ratio of the relative height of the cell in the center of the V-shaped valley topography to the surrounding wall surfaces and the length of the cell becomes smaller.

[0048] When the average elevation of n × n cells in DEM(1) with Lmin = 1 is calculated to create a new low-resolution DEM, the distribution of normal plan curvature of the cells that show a V-shaped valley topography differs depending on whether n is an even number or an odd number.

[0049] Using FIG. 6, we will explain the case where n is an odd number. FIG. 6(a) shows DEM(1) of a V-shaped valley topography with the valley line extending along the fifth column (j=5) centered on cell C55. FIG. 6(b) shows DEM(L) where L=n×Lmin is obtained by aggregating DEM(1) into n×n cells. For example, the elevation of C11 in DEM(L) is (k-1)an, and if n is 3, it becomes 3(k-1)a. Note that C11 in DEM(L) is the average elevation of 3×3 cells C11, C12, C13, C21, C22, C23, C31, C32, and C33 surrounded by thick lines in DEM(1). The same applies to the other cells in DEM(L).

[0050] The case where n is an even number will be explained using Figures 7 to 9. The 9x9 cells shown with solid lines in Figure 7 are the same as those in Figure 6(a), but only a portion of the elevation is shown. The 3x3 cells shown with dotted lines represent DEM(2) with L=2 (2xLmin) obtained by aggregating DEM(1) into 2x2 cells, and this DEM(2) represents a V-shaped valley topography with the valley line running along the second column (j=2). DEM(2) is positioned with a shift of half the cell length L=1 from DEM(1).

[0051] Figure 8(a) is an extract of DEM(2) shown in Figure 7, showing the elevations calculated from the cells of DEM(1). Figure 8(b) shows a generalization of Figure 8(a), that is, a lower-resolution DEM(L) for n being an even number.

[0052] The calculation of the elevation in Figure 8(a) is performed as follows. The 9x9 cells shown by the solid lines in Figure 9(a) are the same as those in Figure 6(a), but only a portion of the elevation is shown. The dotted lines in Figure 9(a) represent DEM(1'), which consists of 6x6 cells, placed in a position overlapping DEM(2) shown in Figure (7). DEM(1') shown by the dotted lines in Figure 9(a) is placed offset by half the length of a cell from the solid line DEM(1).

[0053] The elevation of each cell CIJ (indicated by row I and column J) of DEM(1') shown by dotted lines in Figure 9(b) is the average elevation of the four cells of DEM(1) shown by solid lines in Figure 9(a) that overlap with cell CIJ. For example, cell C11 (I=1, J=1) of DEM(1') shown in Figure 9(b) overlaps with four cells C22 (i=2, j=2), C23 (i=2, j=3), C32 (i=3, j=2), and C33 (i=3, j=3) of DEM(1) shown in Figure 9(a). The elevation of cell C11 of DEM(1') is (k-1)5a / 2, which is the average elevation of the four cells C22, C23, C32, and C33 of DEM(1).

[0054] Then, by calculating the average elevation of the 2 x 2 cells included in DEM(1') shown by the dotted line, the elevation of DEM(2) shown in Figure 8(a) can be obtained. For example, if the average of the elevations of cells C11, C12, C21, and C22 in DEM(1') in Figure 9(b) is calculated, it becomes 2(k-1)a, which is the elevation of cell C11 in DEM(2) in Figure 8(a).

[0055] The method of calculating DEM(2) arranged with a shift of half the cell length from DEM(1) as shown in FIG. 9 can also be applied to the case of calculating the elevation Hij in FIG.

[0056] As shown in Figure 6(b) and Figure 8(b) above, the elevation differs depending on whether n in the n x n cells used for reducing the resolution is an even number or an odd number. If the normal plan curvature is corrected taking into account this difference between even and odd numbers, the normal plan curvature in the cell at the center of the V-shaped valley can be expressed as follows:

number

[0057] m is an integer. When n is odd, the absolute value of the normal plan curvature depends on the number of aggregates, and quickly converges to a minimum value of 1.5k. For example, a V-shaped valley topography with L=1 and k=1 shows a Kplan value of approximately 1.55 when aggregated to L=3. On the other hand, when n is even, the valley line passes between two cells of the minimum cell size Lmin, the two elevations on either side of the valley line are the same, and the normal plan curvature is always 1.5k.

[0058] For example, when the folding angle is 90 degrees, that is, when k = 1, the absolute value of the normal plan curvature obtained from formula (3) is 1.5 when n is an even number, and when n is an odd number, it also converges to the minimum value of 1.5 × 1. Therefore, if -1.5, which is the normal plan curvature of a V-shaped valley topography with a folding angle of 90 degrees, is set as the first threshold, cells with normal plan curvatures equal to or less than the first threshold are extracted as depression areas with folding angles of 90 degrees or less.

[0059] By changing the first threshold value for the normal plane curvature for determining that the central cell in a group of cells is a depression area according to n of the n x n cells used when reducing the resolution, it is possible to extract depression areas taking into account the effects of averaging even in a low-resolution DEM obtained by averaging high-resolution DEMs.

[0060] The first threshold determined in this way differs depending on the n of the n x n cells used when reducing the resolution, but it can be used commonly for multiple DEM(L)s that have been reduced in resolution using the same n x n cells.

[0061] The aggregation step, plan curvature calculation step, normal plan curvature calculation step, and depression determination step described above can be summarized as follows. First, a DEM(L) is created by lowering the resolution of the high-resolution reference DEM(1) in multiple stages. Next, the normal plan curvature is calculated for each cell contained in each DEM(L).

[0062] Next, the normal plan curvature when the contour line folding angle is 90 degrees and 53 degrees, which is judged to be an ideal V-shaped valley topography, is calculated from the above formula (3) and set as the first threshold. For example, if the folding angle is 90 degrees, k = 1 (see Figure 5 (b)), and the value obtained from formula (3) is set as the first threshold. Similarly, if the folding angle is 53 degrees, k = 2 (see Figure 5 (a)), and the value obtained from formula (3) is set as the first threshold.

[0063] Next, cells whose normal plane curvature calculated for each DEM(L) is equal to or less than the first threshold value are extracted as depression regions. Figure 10 shows an example of the extracted depression regions.

[0064] In Figures 10(a)-(f), a first threshold of -1.5 is used, and cells with normal plan curvatures below the first threshold, i.e., cells determined to be depression areas, are shown as black squares. The background of Figure 10(a) shows contour lines spaced at 10m intervals, and the background of Figures 10(b)-(f) shows a 2m mesh map of the amount of slope.

[0065] In Figure 10(a), the length L of a cell's side is 10 m, i.e., the cells determined to be depression areas in DEM(10) are shown as black rectangles. In Figure 10(b), DEM(2) to DEM(10), which are spaced 1 m apart, are grouped together, and the cells determined to be depression areas in each DEM are shown overlapping with black rectangles. Figures 10(c) to (f) show the cells determined to be depression areas in the same way as in Figure 10(b) as overlapping with black rectangles for each group: DEM(11) to DEM(20), DEM(41) to (50), DEM(91) to DEM(100), and DEM(191) to DEM(200).

[0066] In FIG. 10(a), the depression areas are scattered compared to FIG. 10(b)-(c), but in FIG. 10(b)-(c), it can be seen that the depression areas are extracted in a continuous manner. In FIG. 10(d)-(f), the cells become large and extend beyond the ridge. Therefore, the upper limit of the cell length is appropriately set so that the cell does not extend beyond the ridge (step S5 in FIG. 1). In the figure, for example, the upper limit of the cell length is set to L=20 m. Note that whether or not the depression area extends beyond the ridge may be determined by a person based on the positional relationship between the depression area and the contour line, or by a computer using image processing technology.

[0067] Next, a valley region is determined from the depression region obtained as described above. Considering that a valley topography is a topography in which a depression region continues in the direction of a falling water line, the valley region can be extracted from the depression region based on the slope angle and length.

[0068] First, the depression areas extracted for multiple DEMs belonging to one group are overlaid (see Figure 10). Next, the slope angle is calculated for each 1m cell in the multiple DEMs belonging to one group. In other words, multiple slope angles calculated for each DEM (L) are obtained for one 1m cell. Next, from the depression areas obtained by overlaying, cells whose minimum slope angle is equal to or greater than a third threshold are extracted (step S4 in Figure 1), and a polygon consisting of the extracted contiguous cells is created.

[0069] Taking the group of cells in Fig. 5(a) as an example, cell C22 is determined to be a depressed area, and the inclination angle of cell C22 is calculated based on all of its surrounding cells C11 to C33. This method of calculating the inclination angle is a known method, so a detailed description will be omitted.

[0070] FIG. 11 shows depression areas and their slope angles. FIG. 11(a) shows the depression areas extracted by DEM(2) to DEM(10) superimposed on each other. For each 1m cell included in the depression area, the smallest slope angle among the slope angles calculated using DEM(2), DEM(3), ..., and DEM(10) is adopted as the slope angle to be compared with the third threshold. In the example shown in the figure, the cells are colored according to the range to which the smallest slope angle belongs, such as 0 degrees to 2 degrees, 3 degrees to 4 degrees, ..., etc. Similarly, FIG. 11(b) shows the cells for DEM(11) to DEM(20) colored according to the smallest slope angle.

[0071] The third threshold value, which is compared with the slope angle of the cell, is the slope angle that is deemed to be a valley. This third threshold value is set by visually judging the valley and its slope angle displayed on a map, as shown in Figure 11. For example, when a cell is judged to be a valley based on the contour lines overlaid on the DEM, the third threshold value is set so that the slope angle of that cell is used. In addition, when a cell is judged to be a ridge or forest road based on the contour lines overlaid on the DEM, the third threshold value is set so that the slope angle of that cell is not used. In the example in Figure 11, 7 degrees is set as the third threshold value for the slope angle.

[0072] Furthermore, the slope angle of a valley is considered to be steeper than the slope angles of ridges and forest roads that are not valleys. When a certain slope angle is set as the third threshold, it is considered that the area (number of cells) extracted as a valley changes significantly as a result of ridges and forest roads being excluded. Therefore, the third threshold can be determined based on the change in the area of ​​the extracted valley when the third threshold is changed.

[0073] Since the inclination angle at which a terrain is determined to be a valley varies depending on the region, it is preferable to set the third threshold value for the inclination angle for each region.

[0074] Next, a method for extracting valley regions based on length will be described with reference to Fig. 12 (valley region extraction step; step S6 in Fig. 1). First, multiple DEMs with reduced resolution are grouped according to the cell lengths. For example, DEMs with cell lengths of 2 to 10 m, 11 to 20 m, 21 to 30 m, ... are grouped. Then, in any of the DEMs in the group, cells with a normal plan curvature equal to or less than a first threshold and a tilt angle equal to or greater than a third threshold are extracted and made into polygons.

[0075] FIG. 12 shows an example of polygon P consisting of six cells extracted from a group (DEM(2) to DEM(10)) with cell lengths of 2 to 10 m. If the cells are the same size, a polygon with a large area will be a thin strip, but a small polygon will approximate an ellipse. Furthermore, the length of the longest cell contained in polygon P is 4 m. The area of ​​polygon P can be approximated to (4 / 2) × (17 / 2) × π = 17π, which is the area Sp of the ellipse.

[0076] On the other hand, in a 3x3 cell group where the length of one cell is the maximum cell length of polygon P, the shortest valley length is considered to be three cells. The area Sth of an ellipse inscribed in three cells is 3πL 2 / 4.

[0077] The area Sth divided by the square of the cell length L is approximately 2.4. The approximate area Sp of polygon P is L 2 If the value obtained by dividing Sp by L is 2.4 or more, then polygon P can be determined to be a valley longer than the shortest valley. 2 The value obtained by dividing the distance between the polygon P and the valley region is approximately 3.3, which is greater than or equal to the second threshold value of 2.4. Therefore, the polygon P can be determined to be a valley region having a sufficient length.

[0078] The steps for extracting valley regions from depression regions can be summarized as follows. First, for the depression regions of the grouped multiple DEMs (L), extract those with a minimum slope angle of the third threshold or more for each cell. Then, overlap the depression regions (cells) with slope angles of the third threshold or more extracted for each of the multiple DEMs (L), and extract continuous depression regions as polygon P. If the length of polygon P in the valley line direction is equal to or greater than the second threshold, the polygon P is extracted as a valley region. A specific method for determining whether the length of polygon P is equal to or greater than the second threshold is as follows. The area Sp of an ellipse that approximates polygon P is calculated, and the value obtained by dividing the area Sp by the square of the length of the largest cell included in polygon P is calculated. If the value is equal to or greater than the second threshold, 2.4, polygon P is a valley region, and if it is less than 2.4, it is excluded as a local depression rather than a valley region.

[0079] Next, a water system is created (water system extraction step; step S7 in FIG. 1). Specifically, the valley region created as described above is divided into cells of a size set for analysis (hereinafter referred to as analysis cells), a constant value is assigned to each analysis cell, and the accumulated flow rate is calculated from upstream to downstream. The analysis cell whose accumulated flow rate exceeds the fourth threshold is regarded as a water system. This water system is based on the analysis cell that first exceeds the fourth threshold. A known method can be used to calculate the accumulated flow rate, so a detailed explanation will be omitted.

[0080] Figure 13 shows the accumulated discharge resulting from a valley region overlapping near the starting point of a V-shaped valley. The analysis cell is a square cell with one side measuring l. L represents the width of the valley. The ratio of these, L / l, is taken as w. The dotted valley line represents the valley line, and the grey cell represents the valley zone. The figure shows the accumulated discharge for ratios of 3, 5, and 7. If a discharge of 1 is assigned to the analysis cell included in the valley region and the accumulated discharge is calculated, it will be as shown in the figure. The accumulated discharge at the valley line has been extracted to the lower side.

[0081] When a flow rate of 1 is assigned to the valley region overlapping near the starting point of the V-shaped valley topography as shown in FIG. 13, the cumulative flow rate v at the starting point of the V-shaped valley topography is expressed by Equation (4).

number

[0082] As shown in Figure 13, the cumulative flow rate at the start point of the V-shaped valley topography is 4 when w = 3, 14 when w = 5, and 30 when w = 7. Therefore, by using the cumulative flow rate expressed by formula (4) as the fourth threshold, cells whose cumulative flow rate is equal to or greater than the fourth threshold (the consecutive cells extracted at the bottom in the figure) become part of a river system whose base point (valley head) is the cell that first exceeded the fourth threshold (the cells with cumulative flow rates of "4," "14," and "30" in the example in the figure).

[0083] For example, when calculating the cumulative flow rate for the valley area obtained from the group DEM(2) to DEM(10), l = 5, L = 10, and w = 2, use 1 as the fourth threshold value from formula (4) (the decimal points of the calculation result of the formula are rounded down). Similarly, for the valley area obtained from the group DEM(11) to DEM(20), use l = 5, L = 20, w = 2, and 8 as the fourth threshold value (the decimal points are also rounded down), and for the valley area obtained from the group DEM(21) to DEM(30), use l = 5, L = 30, w = 6, and 21 as the fourth threshold value (the decimal points are also rounded down). In this way, a river system starting from the most upstream position of the valley area can be obtained for each DEM group.

[0084] As shown in Figure 14, the river systems obtained as described above may have accumulated flow rates above the fourth threshold near bends and confluences. The cells in the figure represent valley regions, and the shaded cells represent river systems. Cells marked with an "X" indicate cells that have accumulated flow rates above the fourth threshold near bends and confluences. Thus, even within a single valley region, many short valleys (cells marked with an "X") branch off from the main valley.

[0085] Therefore, the length l of the valley width cell is set as the fifth threshold, and valleys separated by valley order whose lengths are longer than the fifth threshold are extracted (water system extraction step; step S8 in Figure 1). The valley order is classified using the Strahler method, and when valleys of the same order merge, the order of the merged valley becomes the valley order +1. In the example of Figure 14, the valley L1 remains as a water system because it is longer than the fifth threshold l. L1, L3, L4, and Lb are primary valleys, and La is a secondary valley, and their lengths are judged individually. L1, La, and Lb remain as water systems because they are longer than the fifth threshold l. The valleys marked with "x" in L3 and L4 are below the fifth threshold, so are excluded from the water system. La and Lb ultimately remain as L2, which is a water system that forms one primary valley. The fifth threshold is, for example, 5 m for DEM(2) to DEM(10), 15 m for DEM(11) to DEM(20), and 25 m for DEM(21) to DEM(30).

[0086] Next, the water system is thinned (water system extraction step; step S8 in FIG. 1). Thinning will be explained with reference to FIG. 15. The cells in this figure represent valley regions, and among them, the hatched cells represent water systems. As shown in FIG. 15(a), the water system includes two flow paths that extend in one direction (the vertical direction in the figure) and have a width of one cell. One of these flow paths is left, and the other is removed. As a result, a thinned water system consisting of a single flow path can be obtained, as shown in FIG. 15(b).

[0087] The process described in Figures 14 and 15 removes the starting points of valleys that occur near the main stream. In other words, it is possible to extract only the starting points of valleys that are identified on the most upstream side of the valley. Using the starting points of the valleys extracted in this way, the cumulative flow rate is calculated again, and a river system that is in harmony with the valley width distribution map can be obtained (water system extraction step; step S9 in Figure 1).

[0088] Figure 16 shows the results obtained by extracting valley topography and water systems from DEM. In Fig. 16(a), a cell consisting of cells painted black is extracted as a depression area. Fig. 16(b) shows a valley area (polygon) obtained by extracting the depression area based on the second threshold value of the slope angle and length. The cells painted red have a slope angle of 0 degrees to 7 degrees, and the cells painted blue have a slope angle of 7 degrees to 40 degrees. It can be seen that the valley area with a slope angle of 7 degrees to 40 degrees matches the actual valley topography.

[0089] Figure 16(c) shows a river system consisting of green cells that show accumulated flow rates above the fourth threshold, calculated for valley regions with slopes between 7 and 40 degrees. Figure 16(d) shows the river system with short valleys removed in red. Figure 16(e) shows the thinned river system in yellow. Figure 16(f) shows the light blue river system obtained by calculating accumulated flow rates from the valley heads shown in blue cells.

[0090] The valley topography extraction method and program described above calculates the normal plan curvature for multiple DEMs (L) with reduced resolution. The normal plan curvature is a value that does not depend on the cell length L of the DEM. Therefore, the first threshold value set for the normal plan curvature can be commonly used for multiple DEMs (L) with different resolutions, and V-shaped valley topography can be objectively extracted as a depression area. In addition, the valley topography extraction method extracts valley regions based on the length of the depression region and the second threshold value, which makes it possible to exclude isolated depression regions and objectively extract valley regions having an appropriate length as a V-shaped valley topography. Furthermore, the valley topography extraction method extracts a polygon obtained by overlapping depression areas obtained for each of a plurality of DEM(L) as a valley area. For example, a V-shaped valley topography with a gentle cross-sectional slope cannot be extracted as a depression area with a DEM(L) with a short L, but may be extracted as a depression area with a DEM(L) with a long L. In the valley topography extraction method of the present invention, depression areas obtained from DEM(L) with a plurality of resolutions are overlapped to form a polygon (valley area), so that V-shaped valley topography of various shapes is unlikely to be overlooked, and valley areas can be extracted comprehensively. Therefore, the valley area extracted by the valley topography extraction method is adapted to the actual valley topography.

[0091] In addition, the valley topography extraction method creates a valley area (polygon) for a depression area whose slope angle is equal to or greater than a third threshold value, thereby making it possible to obtain a valley area that is more closely matched to the actual valley topography.

[0092] In addition, the valley topography extraction method extracts the valley region based on the area of ​​the polygon and the second threshold value, thereby obtaining a valley region that is more closely matched to the actual valley topography.

[0093] In addition, the valley topography extraction method uses formula (3) described in this specification as the absolute value of the first threshold value, which makes it possible to extract depression areas more objectively.

[0094] In addition, the valley topography extraction method determines the length of one side of the cells used for aggregation so that the depression area does not exceed the ridge. This makes it possible to obtain a depression area that is more suitable for the actual valley topography.

[0095] The method and program for extracting a river system extracts a river system from a valley area obtained by the valley topography extraction method. The extraction of the river system is performed based on the accumulated flow rate and the fourth threshold value. This makes it possible to obtain a river system that starts from the head of a V-shaped valley. In addition, since it is only necessary to calculate the accumulated flow rate within the range of the valley area obtained by the valley topography extraction method, the river system can be obtained more efficiently than by calculating the accumulated flow rate from the entire DEM. Additionally, drainage systems are extracted by setting the ratio w of the length l of the analysis cell to the width L of the valley area, and calculating the cumulative flow rate for each ratio w as shown in Figure 13. Even if the valley head cannot be detected at a certain ratio w, it may be possible to detect it at another ratio w. Therefore, even for V-shaped valley topography with various shapes, the valley head is unlikely to be overlooked, and it is possible to comprehensively extract the valley head and extract the drainage system starting from that valley head. In addition, since a water system is created for a valley area that matches the valley topography, it is possible to avoid extracting areas that are not valley topography but have a large water collection area as a water system as in the conventional technology, and it is possible to extract a water system that does not include horizontal surfaces created by artificial structures. In addition, since the drainage systems are created based on valley areas created from DEM(L) of different resolutions, different drainage systems are not created depending on the map scale.

[0096] In addition, the method for extracting river systems excludes short valleys based on the length of the valley and the fifth threshold. This makes it possible to exclude river systems that originate from short valleys. Furthermore, the river systems are thinned, and the cumulative flow rate is calculated again from the most upstream of the river system, and river systems consisting of analysis cells that exceed the fourth threshold are extracted. This makes it possible to obtain river systems based on the main valley topography in the region.

[0097] Each step in the valley topography extraction method and the water system extraction method may be implemented as a program, and the program may be executed by a computer, thereby realizing each step on the computer.

[0098] The program may be recorded on a computer-readable recording medium. Examples of computer-readable recording media include magnetic recording devices, optical disks, magneto-optical recording media, and semiconductor memories. The program may be downloadable via a communication means such as the Internet. [Industrial Applicability]

[0099] The maps can be used to represent mountainous terrain, and can also be used in industrial sectors to assess the risk of slope disasters.

Claims

1. A method for extracting valley topography using a digital elevation model (Lmin) in which the earth's surface is divided into square cells each having a side length of Lmin and an elevation is assigned to each cell, comprising the steps of: an aggregation step of calculating an average elevation for n×n cells included in the digital elevation model (Lmin) and creating a low-resolution digital elevation model (L) by assigning the average elevation as the elevation of a square cell whose side length is L (Lmin×n); a plan curvature calculation step in which a step of extracting 3×3 cells (hereinafter referred to as cell groups) included in the digital elevation model (L) and calculating the plan curvature of the center cell for all the extracted cell groups is executed for each of the multiple digital elevation models (L) created in the aggregation step; a normal planar curvature calculation step of calculating a normal planar curvature by multiplying the planar curvature by the length L of a side of a cell for each of the plurality of digital elevation models (L) created in the aggregation step; a depression determination step of setting a first threshold value for normal planar curvature for determining that the central cell of the group of cells is a depression area, the first threshold value being common to a plurality of digital elevation models (L), and determining that a cell having the normal planar curvature equal to or less than the first threshold value is a depression area; a valley area extraction step of overlaying the depression areas extracted for each of the plurality of digital elevation models (L) to create polygons consisting of continuous depression areas, and extracting the polygons whose lengths in the valley line direction are equal to or greater than a second threshold as valley areas. A method for extracting a valley topography, comprising:

2. 2. The method for extracting a valley topography according to claim 1, comprising the steps of: The valley region extraction step includes overlapping the depression regions extracted for each of the plurality of digital elevation models, extracting cells from the overlapping depression regions in which the minimum slope angle in the digital elevation models is equal to or greater than a third threshold value, and creating the polygon consisting of the extracted consecutive cells. A method for extracting a valley topography, comprising:

3. The method for extracting a valley topography according to claim 1 or 2, comprising: The valley region extraction step extracts the polygon as a valley region if a value obtained by dividing the area of ​​the polygon by the square of the length of the maximum cell included in the polygon is equal to or greater than the second threshold value of 2.

4. A method for extracting a valley topography, comprising:

4. A method for extracting a valley topography according to any one of claims 1 to 3, comprising: In the depression determination step, the absolute value of the first threshold value is determined by using the formula (1). A method for extracting a valley topography, comprising: [0010] Here, Kplan,normal is the first threshold, m is an integer, and k is a magnification of the elevation difference a in the valley line direction and the elevation difference in the direction perpendicular to the valley.

5. A method for extracting a valley topography according to any one of claims 1 to 3, comprising: In the aggregation step, the length L of one side of the cell used for aggregation is determined so that the depression area obtained in the depression determination step does not exceed the ridge. A method for extracting a valley topography, comprising:

6. A method for extracting a water system using a valley area obtained by the valley topography extraction method according to any one of claims 1 to 5, comprising the steps of: a water system extraction step of dividing the valley region into square analysis cells each having a side length of 1, calculating an accumulated flow rate by applying a constant value to each analysis cell, and extracting a water system consisting of analysis cells whose accumulated flow rate exceeds a fourth threshold value; The fourth threshold value is the cumulative flow rate at the head of the V-shaped valley.

1. A method for extracting an aqueous system comprising the steps of:

7. 7. A method for extracting an aqueous system according to claim 6, comprising the steps of: The water system extraction step extracts a water system consisting of analysis cells in which the length of a valley divided by a valley order is longer than a fifth threshold value.

1. A method for extracting an aqueous system comprising the steps of:

8. 8. A method for extracting an aqueous system according to claim 7, comprising the steps of: The water system extraction step includes thinning the water system, calculating the cumulative flow rate again from the most upstream of the water system, and extracting a water system consisting of analysis cells whose cumulative flow rate exceeds the fourth threshold value.

1. A method for extracting an aqueous system comprising the steps of:

9. A method for causing a computer to execute each step of the valley topography extraction method according to any one of claims 1 to 5. A program characterized by:

10. A method for causing a computer to execute each step of the aqueous extraction method according to any one of claims 6 to 8. A program characterized by:

Images