System and program for stereoscopic visualization of crustal deformation

Abstract

To obtain a crustal deformation stereoscopic image processing system that can accurately color-code images based on elevation differences that take into account vertical and horizontal components, even when crustal deformation causes displacement in both vertical and horizontal directions. [Solution] A red stereoscopic map Gbi (image) is generated based on a new DEM of a predetermined area (S12), and a grayscale red stereoscopic map GBgi is generated (S16). On the other hand, the new DEM (post-earthquake) and the pre-earthquake (old) DEM (old DEM) are differentiated (simple difference) to highlight the areas of change (S24). Then, these areas of change are color-coded according to the degree of subsidence and uplift (S26).

Description

[Technical Field]

[0001] This invention relates to a system for processing stereoscopic images of crustal deformation. [Background technology]

[0002] After an earthquake, aerial photographs are taken to assess the extent of damage, including fires, liquefaction, building collapses, structural damage, landslides, and damage to rivers, ports, airports, and other infrastructure. In addition, changes are sometimes analyzed using observational data from multiple satellite SARs before and after an earthquake [methods using Synthetic Aperture Radar (SAR)].

[0003] This method allows for obtaining displacement amounts on the order of several centimeters over a wide range by analyzing the results using differential interferometric synthetic aperture radar (hereinafter simply referred to as "interferometric SAR"), which compares the SAR measurement results from two different time periods.

[0004] Specifically, this method utilizes the phase difference obtained from observation results at two different time periods to understand changes in topography. Generally, the results of interferometric SAR are represented as a striped SAR interferogram. The SAR interferogram is created by first dividing the phase difference into multiple ranges (hereinafter referred to as "phase difference ranges"), classifying the interferometric SAR results (i.e., phase differences) into these phase difference ranges, and assigning a specific color (or grayscale) to each phase difference range within the same phase difference range.

[0005] Because it can acquire displacement data over a wide area at once and periodically, interferometric SAR is used in a variety of applications, and various technologies utilizing interferometric SAR have been proposed to date.

[0006] For example, Patent Document 1 (Patent No. 7335733: Building Damage Estimation System: Kokusai Kogyo) proposes an invention that calculates the displacement of a target area with high accuracy and over a wide range by utilizing the measurement results of interferometric SAR and a satellite positioning system (GNSS: Global Navigation Satellite System).

[0007] Specifically, based on the phase and intensity of radio waves at two different time points obtained by synthetic aperture radar, an interference level representing the coherence of radio waves at those two time points is calculated. Then, the interference level ratio is calculated between the target interference level at the time when building damage is to be estimated and the reference interference level calculated based on the two time points during normal conditions.

[0008] Then, by comparing this interference level ratio with a threshold for damage estimation, the presence or absence of building damage is estimated, the estimated damage locations are detected, and the scale of damage at those estimated damage locations is estimated.

[0009] The aforementioned estimation of damage locations is performed by calculating n types of interference levels based on observations from synthetic aperture radar at n wavelengths (where n is a natural number greater than or equal to 3), calculating interference level ratios for each of the n types of interference levels to detect damage locations, and then estimating the scale of damage in n stages ("high," "medium," "low") for each damage location based on the overlap of these damage locations based on the n types of interference level ratios.

[0010] Patent Document 2 (Patent No. 4339289: NEC System Technology Corporation) By performing stereo processing on multiple images of a predetermined region at a first time point (before the earthquake) and a second time point (after the earthquake), 3D data from the first time point (before the earthquake) and the second time point (after the earthquake) is extracted.

[0011] Then, changes in the target are determined from the orthomosaic images and orthomosaic 3D data at the first time point (before the earthquake) and the second time point (after the earthquake).

[0012] The determination of this change is made by comparing the heights of the ortho 3D data at the first time point (before the earthquake) and the ortho 3D data at the second time point (after the earthquake), and by comparing the elevation of the ground surface with the height.

Prior Art Documents

Patent Documents

[0013]

Patent Document 1

Patent Document 2

Summary of the Invention

Problems to be Solved by the Invention

[0014] However, the one in Patent Document 1 detects the damage estimation location based on the phases and intensities of the radio waves at two times obtained by synthetic aperture radar, so there is a problem with the accuracy of the estimated location.

[0015] That is, it is not possible to accurately estimate the earthquake damage locations considering the vertical and horizontal components.

[0016] On the other hand, although it is described in Patent Document 2 that the heights of the ortho 3D data at the first time point (before the earthquake) and the ortho 3D data at the second time point (after the earthquake) are compared, and the elevation of the ground surface is compared with the height for determination, there is no description regarding the estimation of the earthquake damage locations displaced in the vertical and horizontal directions.

[0017] For example, in the Noto Peninsula earthquake, crustal movements (mainly rigid XYZ direction movements) and landslides (blocky or fragmented movements) occurred simultaneously on the ground surface.

[0018] The present invention has been made in view of the above problems, and even when the crust moves vertically and horizontally, it can accurately colorize with elevation differences considering vertical and horizontal components, so that crustal movements caused by natural disasters such as earthquakes can be shown with higher accuracy than conventional techniques. The purpose is to obtain a crustal movement stereoscopic visualization image processing system.

Means for Solving the Problems

[0019] The crustal movement stereoscopic visualization image processing system according to the present invention includes: (A) a storage unit that stores an old DEM at a predetermined time and a new DEM whose acquisition time is newer than the predetermined time; (B) means for matching the horizontal components of the point cloud of the old DEM and the point cloud of the new DEM; (C) means for obtaining a simple elevation difference value between the old DEM and the new DEM after the horizontal components are matched; (D) means for generating a simple elevation difference image with a first color indicating no difference if the simple elevation difference value is within a predetermined value, a second color gradient indicating uplift if it is greater than or equal to the predetermined value, and a third color gradient indicating subsidence if it is less than or equal to the predetermined value, and displaying it on the screen; (E) means for generating a red stereoscopic image based on the new DEM; (F) means for generating a red grayscale image obtained by grayscale processing of the red stereoscopic image; (G) It is characterized by having means for generating a crustal displacement stereoscopic image obtained by multiplying and synthesizing the simple elevation difference image and the red grayscale image and displaying it on the screen.

[0020] Also, the crustal movement stereoscopic visualization image processing program according to the present invention causes a computer to, (A) means for storing an old DEM at a predetermined time and a new DEM whose acquisition time is newer than the predetermined time in a storage unit; (B) means for matching the horizontal components of the point cloud of the old DEM and the point cloud of the new DEM; (C) means for obtaining a simple elevation difference value between the old DEM and the new DEM after the horizontal components are matched; (D) Means for generating and displaying a simple elevation difference image on a screen, the first color indicating no difference if the simple elevation difference value is within a predetermined value, the second color gradient indicating uplift if it is above the predetermined value, and the third color gradient indicating subsidence if it is below the predetermined value. (E) Means for generating a red stereoscopic image based on the new DEM, (F) Means for generating a red grayscale image by converting the red 3D image to grayscale, (G) Means for generating a crustal displacement stereoscopic image by multiplying and combining the simple elevation difference image and the red grayscale image, and displaying it on a screen. To perform its function as such. [Effects of the Invention]

[0021] As described above, according to the present invention, even if the Earth's crust is displaced in both vertical and horizontal directions, the elevation difference, taking into account the vertical and horizontal components, can be accurately colored. In other words, even in cases involving widespread crustal deformation caused by natural disasters such as earthquakes, the location and scale of landslides can be recognized with greater precision. Furthermore, it can eliminate the misidentification problem, such as "shading," that can occur with simple elevation difference methods in phenomena like the Noto Peninsula earthquake, where crustal deformation (rigid body movement) and localized landslides or collapses occur simultaneously. [Brief explanation of the drawing]

[0022] [Figure 1] This is a flowchart illustrating the overview of the crustal deformation stereoscopic image processing system of this embodiment 1. [Figure 2] This is an explanatory diagram illustrating the overview of the crustal deformation stereoscopic image processing system of this embodiment 1, showing the progression of images. [Figure 3] This is an explanatory diagram of a red stereoscopic map after an earthquake. [Figure 4] This is an explanatory diagram of a grayscale red stereoscopic map. [Figure 5] This is an explanatory diagram of the image after color adjustment following the multiplication blending process. [Figure 6]This is a program block diagram of the crustal deformation stereoscopic image processing system of this embodiment 1. [Figure 7] This flowchart (1) supplements the processes of differential processing and coordinated display of terrain change images. [Figure 8] This flowchart (2) supplements the processes of differential processing and coordinated display of terrain change images. [Figure 9] This is an explanatory diagram illustrating an example of a table representation of pre-earthquake patterns (PTai) and post-earthquake patterns (PTbi). [Figure 10] This is an explanatory diagram illustrating the relationship between the elevation difference and the gradient scale between the pre-earthquake pattern PTai and the post-earthquake pattern PTbi. [Figure 11] This is a three-view drawing of a trapezoidal terrain. [Figure 12] This is an explanatory diagram of vectors resulting from the displacement of a trapezoidal terrain. [Figure 13] This is an explanatory diagram illustrating the difference in the amount of uplift before and after an earthquake. [Figure 14] This is an explanatory diagram for a simple difference diagram. [Figure 15] This is an explanatory diagram of the shift difference diagram. [Figure 16] This is a composite image of a grayscale red relief map and a shift difference diagram. [Figure 17] This is an explanatory diagram of the image after color correction. [Figure 18] This is an explanatory diagram of images showing topographic changes after an earthquake caused by the processing of this embodiment. [Figure 19] This is an explanatory diagram showing how local areas are color-coded according to their degree of subsidence and uplift. [Figure 20] This is a flowchart illustrating the crustal deformation stereoscopic image processing system of Embodiment 2. [Figure 21] This is an explanatory diagram of the image when using the difference obtained from DSM. [Figure 22] This is a magnified view of a portion of Figure 21. [Figure 23] This is an explanatory diagram of another example image resulting from this processing. [Figure 24]This is an explanatory diagram of the image when combined with orthophoto. [Figure 25] This is an explanatory diagram illustrating the process of generating a red 3D image. [Figure 26] This is a schematic diagram illustrating the estimation of crustal deformation components using the unchanged area Pmi (stable region SC) of Embodiment 2. [Figure 27] This diagram illustrates that, during a major earthquake, simple DEM (Digital Emission Model) differences cannot separate wide-area seismic deformation from localized seismic movement (landslides, collapses). [Figure 28] This diagram illustrates the minimization of the sum of squared elevation differences using the optimal shift amount conversion parameter (α★). [Figure 29] This is an explanatory diagram of the Kumamoto-type optimization. [Figure 30] This is an explanatory diagram for the optimization of the Noto type. [Figure 31] This is an explanatory diagram illustrating the choice between the Kumamoto-type and Noto-type. [Figure 32] This is an explanatory diagram illustrating the correction by the stable region in this second embodiment. [Figure 33] This is an explanatory diagram of a simple elevation difference image obtained using conventional methods. [Figure 34] This is an explanatory diagram of the difference image after fluctuation correction using the maintenance method. [Figure 35] This diagram illustrates how the sediment balance (erosion rate and deposition rate) can be determined using this method. [Figure 36] This is an explanatory diagram illustrating the process of canceling crustal deformation by moving the unchanging area Pmi (SC or reference area). [Figure 37] This is an explanatory diagram illustrating the search for the optimal shift amount. [Figure 38] This is an explanatory diagram illustrating the effect of alignment on eliminating crustal deformation. [Figure 39] This flowchart supplements the processing of coordinated display of terrain change images. [Figure 40] This is a diagram illustrating the specific process of specifying the reference region and searching for the optimal shift. [Modes for carrying out the invention]

[0023] This embodiment makes it possible to more precisely recognize the location and scale of landslides, even when they are accompanied by widespread crustal deformation caused by natural disasters such as earthquakes. Crustal deformation (earthquakes, etc.) can be categorized into two types: the Noto earthquake type, where vertical deformation is dominant, and the Kumamoto earthquake type, where horizontal deformation is dominant. Furthermore, these two types can sometimes coexist. Regardless of the type, this is a crustal deformation stereoscopic image processing system that can display crustal deformation with higher accuracy than conventional technologies.

[0024] For phenomena such as the Noto Peninsula earthquake (where vertical movement was dominant), where crustal deformation (rigid body movement) and localized landslides or collapses occur simultaneously, this method eliminates the misidentification problem, such as "shading," that occurs with simple elevation difference methods, and accurately extracts the location and displacement of landslides and collapses. An embodiment that eliminates this misidentification problem, such as shading, will be described in Embodiment 2. Embodiment 1 describes the case where the vertical component is dominant.

[0025] Regardless of the type of earthquake described above, information is used from areas that have not experienced crustal deformation (also called unaffected areas Pmi) near areas that are uplifting or sinking due to self-sliding collapses, etc. (referred to as local areas or altered areas) (even Pmi areas are affected). In this embodiment, the background topography is matched (approximated) by aligning (approximating) the mesh within the unchanged area Pmi (new DEM) outside the collapsed area (changed region) (background) with the mesh within the area of ​​the old DEM corresponding to this Pmi. Subsequently, coloring is performed based on elevation difference values, enabling precise estimation of the location of landslides and collapses.

[0026] (Elimination of "shading effect" through proprietary correction) Simply calculating the elevation difference between the old and new DEMs can lead to a "shading effect" where uplift and subsidence are mistakenly identified as being caused by widespread crustal deformation (resulting in a red and blue shaded image when uplift is represented in red, subsidence in blue, and no displacement in white).

[0027] To solve this, in Embodiment 2, which will be described later, a two-step method is employed: first, when it is determined that the image is "shaded," a shift is performed in the horizontal component, i.e., in the X and Y directions, to eliminate the shading effect; and if the image still remains red or blue and does not become white, a further shift is performed in the vertical component, Z direction. This allows for accurate identification of localized changes, such as landslides, after matching the background topography.

[0028] (Combination of crustal deformation model and 3D image) Furthermore, this embodiment employs a unique method of multiplying and combining a "simple elevation difference image," which is colored based on simple elevation difference values, with a "red relief map image" that emphasizes terrain relief, after reducing its saturation to gray. As a result, both elevation changes and terrain relief are combined into a single image, generating a "terrain displacement stereoscopic image" that is very visually easy to understand. This is a forward-thinking method that will greatly contribute to understanding the situation during disasters and to rapid interpretation by experts.

[0029] The embodiments shown below illustrate devices and methods for realizing the technical concept (structure, arrangement) of the invention, and the technical concept of the present invention is not limited to those described below. The technical concept of the present invention can be modified in various ways within the scope of the claims. Furthermore, it should be noted that the drawings are schematic, and the configuration of the devices and systems may differ from those in reality.

[0030] First, let me briefly describe the outline of Embodiment 2. During the Noto Peninsula earthquake, crustal deformation (mainly rigid movement in the XYZ directions) and landslides (movement in chunks or fragments) occurred simultaneously on the ground surface. While it should be possible to precisely extract the distribution of crustal deformation and landslides caused by the earthquake using high-precision (resolution approximately 50 cm) laser-measured DEMs (digital elevation models) acquired before and after the earthquake, this is not easy.

[0031] With simple elevation difference calculations, if there is a large horizontal movement component, it can result in a shading effect (where the side in the direction of movement appears to have risen and the opposite side appears to have subsided due to the building's movement). When there is a large vertical component, the displacement from crustal movement overlaps overall, making it difficult to identify landslides and other geological events.

[0032] (Estimation of crustal deformation and identification of landslides, and identification of completely destroyed houses) To solve this problem, when comparing DEMs before and after an earthquake, we consider widespread crustal deformation as rigid body movement with XYZ components that do not involve rotation.

[0033] Select an area (Pmi) that is free from landslides and other forms of sediment movement. Simply calculate the difference and apply the colormap. For example, a gradient color scheme could be created where +5m is red (second color), 0m is white (first color), and -5m is blue (second color). The value of 5m should be changed as needed. Details about gradients will be explained later.

[0034] If the DSM difference becomes shaded during this process, the XY components are adjusted and shifted to eliminate the shading. The point where no shadow is cast at all is defined as the horizontal movement vector. In that case, if red or blue components remain throughout, a vertical shift is performed to bring it closer to white. Let the amount of that shift be defined as the vertical displacement.

[0035] Afterward, to make the landslide topography easier to see, the area is grayed out and a red relief map is overlaid, and then the area is interpreted and extracted. Note that in areas with a total house destruction distribution, the DSM difference may show many negative values. Furthermore, by using not only red-white-blue but also repeating retardation colors in the color map, a more detailed analysis using moiré patterns becomes possible.

[0036] <Embodiment 1> This embodiment 1 uses a Digital Elevation Model (DEM) of a 1m mesh (or 10cm, 50cm, 2m, 5m, or 10m) for each of two different time periods (for example, before and after an earthquake) to represent the amount of subsidence and uplift in the affected area after an earthquake (or typhoon, natural disaster) using hues proportional to the vertical (Z) and horizontal components of elevation (for example, blue for subsidence, red for uplift). This is then multiplied and combined with a grayscale red image to obtain a predetermined color tone (see Figure 5). Based on the vertical (Z) components of multiple unchanged areas Pmi (arbitrary) outside the affected area (local area), the previous vertical (Z) component is adjusted, and then multiplied and combined again with the grayscale red image to obtain a predetermined color tone (see Figure 5).

[0037] In other words, the image is made clearer by multiplying and compositing it with a grayscale red relief map. The operator specifies the areas to be read from the old and new DEMs.

[0038] Embodiment 1 (also referred to as Scenario) is described as an example of the Noto type where the vertical component is dominant, and Embodiment 2 is described as a scenario in which the Kumamoto type, where the horizontal component is dominant, and the Noto type occur.

[0039] (Embodiment 1) The overview of the crustal deformation stereoscopic image processing system of this embodiment 1 will be explained using Figure 1. Figure 1 is a flowchart illustrating the overview of the crustal deformation stereoscopic image processing system of this embodiment 1, and will be explained as the amount of change after an earthquake.

[0040] As shown in Figure 1, the computer uses a new DEM memory 10B that stores the post-earthquake (new) DEM (referred to as the new DEM) for a predetermined region (for example, Kumamoto or Noto) and an old DEM memory 10A that stores the pre-earthquake (old) DEM (old DEM) for the predetermined region to perform the following processing.

[0041] In other words, as shown in Figure 1, the computer stores the new DEM (post-earthquake) in the new DEM memory 10B and the old DEM in the old DEM memory 10A (S10, S20). The mesh of the new DEM is referred to as mbi, and its mesh number is referred to as mBi. The mesh of the old DEM is referred to as mai, and its mesh number is referred to as mBi. Mai and mbi are collectively referred to as mi. First, we will explain the process of generating a grayscale red image (also known as a DEM red grayscale image).

[0042] As shown in Figure 1, the computer reads a specified area (for example, 500m wide, 200m wide; this area may be larger) from the new DEM (post-earthquake) in memory 10B for the new DEM (S12). The specified area in the new DEM is denoted as Ebi, and the specified area in the old DEM is denoted as Eai (collectively referred to as Ei).

[0043] Then, based on this new DEM (post-earthquake), a red stereoscopic map Gbi (also called a DEM red stereoscopic image; corresponding to Ebi) is generated (Figure 2(a)) and stored in the red memory 50 (not shown in Figure 1) (S14). This generation method is described in Japanese Patent No. 3670274 and is a three-dimensional representation method of terrain that adjusts the amount of slope using the saturation of red and the degree of ridges and valleys using brightness, based on the digital elevation data (DEM). Details of this will be described later.

[0044] Note that the red stereoscopic map Gbi shown in Figures 2(a) and 3 is an example of a case where there is a lot of crustal deformation (vertical direction). Then, the new DEM (post-earthquake) red stereoscopic map Gbi (image) is converted to grayscale (also called grayscale: Figure 2(d)) to generate a grayscale red stereoscopic map GBgi (also called the DEM red grayscale image) (S16: see Figure 4). In the diagram, one example of a change (also called a fluctuation) that we are focusing on is labeled Fhi.

[0045] The computer also reads the pre-earthquake (old) DEM (old DEM) for the designated area Eai from the old DEM memory 10A (S22). Then, the computer performs a simple difference operation (simple difference image Gdi) between the new DEM (post-earthquake) of the designated area Ebi, which was read in step S12, and the pre-earthquake (old) DEM (old DEM) of the designated area Eai, which was read in step S22, and highlights the areas of change (including Fhi) (S24; see Figure 2(b)). The aforementioned differential data (also called differential data) consists of, for example, the old DEM's mesh number mAi (i=1, 2, ...), the old DEM's coordinates CAi (XAi, YAi, ZAi), the new DEM's mesh number mBi, the new DEM's coordinates CBi (XBi, YBi, ZBi), and the difference Δhi, as shown in Figure 9. In Figure 9, the difference Δhi is shown for the X, Y, and Z directions. ΔZhi is the simple elevation difference. Also, the simple elevation value is denoted as ZAi in the old version and ZBi in the new version.

[0046] Then, a coloring process is performed according to the simple elevation difference ΔZhi (S26). For example, areas that have subsided are given a blue gradient according to the degree of subsidence, and areas that have risen are given a red gradient according to the degree of rise. Processing in S32 adjusts the height so that the area around the changed part (especially the Fhi)Cgi of interest becomes white (see Figure 2(c)). In other words, (1) Generate a simple difference image Gdhai (Figure 2(b)) between the elevation value ZAi (old) of the pre-earthquake (old) DEM (Eai: old DEM) and the elevation value ZBi (new) of the post-earthquake (new) DEM (Ebi) (in mesh units (e.g., 1m mesh)). (2) The closer the simple elevation difference ΔZhi between the old and new meshes gets to 0 (no change), the whiter the meshes will be. Specifically, if the simple elevation difference ΔZhi is between 0 and +5m, the meshes will gradually change from white to red, and if it is +5m or more, the meshes will have a red gradient. Also, if the difference is between 0 and -5m or more, the meshes will gradually change from white to blue, and if it is -5m or more, the meshes will have a blue gradient. The color assignment uses a gradient bar Vai, as shown in Figure 2(b), for example.

[0047] The operator then displays an image on the screen that has been color-coded according to the degree of subsidence (simple elevation difference ΔZhi is 0 or greater) or uplift (elevation difference ΔZi is 0 or less), and determines the hue (S28). If it is determined that the overall hue (including the change point Fhi) is incorrect (including other change points), an unchanged area Pmi (also called a stable area SC) is specified, which is assumed to be unchanged, located a certain distance away (nearby, further away area: wide area) from the change point Fhi (preferably a specific change point of interest) (S30). A stable area SC that is roughly the same elevation as the old DEM of the change point Fhi is preferred. Note that the unchanged area Pmi (stable area SC) is an area where landslides, etc., have not occurred. The unchanged area Pmi specifies an area that is larger (in terms of length and width) than the changed area Fhi. If the changed area Fhi is 50m wide and 80m long, it is preferable that the unchanged area Pmi be approximately 70m wide and 100m long.

[0048] Then, a process to shift the Z value (described later) is performed (S32), and the process returns to step S24 to perform the simple difference process again (S32). In the simple difference image Gdi mentioned above, the unchanged area Pmi (or stable area SC) assumed to be unchanged is preferably specified as a different unchanged area Pmi, and the difference is performed multiple times for each specification (preferably until convergence occurs).

[0049] Then, in step S28, if the operator determines that the hue is correct (OK), the color-adjusted change area colored image Gdi' (see Figure 2(c)) is created, and this color-adjusted change area colored image Gdi' and the grayscale red stereoscopic map GBgi obtained in step S16 are multiplied and combined to adjust the color tone (S34: see Figure 5). In this embodiment, the image after adjusting the color tone is called the change area red-blue stereoscopic image Ghi.

[0050] Therefore, for example, the change point Fhi is visually highlighted, making it possible to calculate how much soil and sediment will be produced. Next, a specific program block for performing the process shown in Figure 1 above will be explained using Figure 6. However, in Figure 6, the horizontal component adjustment unit 90 that performs the horizontal component alignment in Embodiment 2 will be shown and explained.

[0051] As shown in Figure 6, the computer main unit 10 includes a new DEM memory 10B that stores the new DEM (after the earthquake), an old DEM memory 10A that stores the old DEM (before the earthquake), a red stereoscopic mapping unit 40, a grayscale conversion unit 60, a horizontal component adjustment unit 90, a difference unit 100, a change area coloring unit 112 (blue for subsidence, red for rise), a hue confirmation unit 116, a difference image reconstruction unit 120, a multiplication and blending unit 130, a color adjustment unit 150, and the like.

[0052] The red stereoscopic mapping unit 40 reads the DEM from the new DEM memory 10B, and generates a red stereoscopic map Gbi (see Figure 3) by adjusting the slope amount with the saturation of red and the ridge / valley degree with the brightness, and stores it in memory 50.

[0053] The grayscale conversion unit 60 converts the red stereoscopic map Gbi (specified area Ebi: after earthquake: see Figure 3) in memory 50 into grayscale for each mesh according to the degree of redness, to obtain the grayscale red stereoscopic map GBgi shown in Figure 4, and stores it in the grayscale red image memory 70.

[0054] The horizontal component adjustment unit 90 aligns the point clouds of the old DEM in the old DEM memory 10A and the new DEM (post-earthquake) in the new DEM memory 10B, using the mesh mi (mai, mbi). (Although horizontal component alignment is not shown in Figure 1, it is performed before discretization.) Note that mai is the mesh mai before the earthquake, and mbi is the mesh after the earthquake.

[0055] The aforementioned point cloud alignment involves identifying corresponding points in the point cloud of the mesh mi (mai, mbi) using the least squares method (CCICP: Classification and Cobined ICP).

[0056] The differential unit 100 calculates the simple elevation difference ΔZhi between the elevation value Zbi of the post-earthquake (new) mesh mbi and the elevation value Zai of the pre-earthquake (old) mesh mai, and stores it in memory 110.

[0057] The color-coded change area 112 applies a gradient (blue for subsidence, red for rise) each time the simple elevation difference ΔZhi is calculated, and stores it in memory 114.

[0058] The hue confirmation unit 116 displays the simple difference image on the screen 205 via the display processing unit 160 each time the hue is stored in the memory 114, and if the hue is deemed OK, it activates the multiplication blending unit 130.

[0059] The multiplication and synthesis unit 130 generates an image by multiplying the simple difference image in memory 114 with the grayscale red stereoscopic map GBgi in memory 70 for grayscale red images, stores it in memory 140, and displays it on the screen by the display processing unit 160.

[0060] The color adjustment unit 150 corrects the hue, saturation, contrast, etc., of the multiply composite image in memory 140 to a color tone suitable for its intended use and displays it on the screen.

[0061] The difference image reconstruction unit 120 determines a shift amount each time a non-changed area Pmi, which is assumed to be unchanged and is close to the changed area of ​​the simple difference image displayed on the screen, is specified, while sequentially specifying the non-changed area Pmi to the simple difference image Gdi or Gdi', and then causes the difference unit 100 to perform simple difference again based on this shift amount. Furthermore, steps S24 (differentiation), S26 (coloring of changed areas), S28 (hue determination), S30 (Pmi), and S32 (Z-shift) in Figure 1 are collectively referred to as the terrain change image coordinated display process.

[0062] In Figure 6, the horizontal component adjustment unit 90, the difference unit 100, the change area coloring unit 112, the difference image construction unit 120 (corresponding to the calculation of the shift amount), the hue confirmation unit 116, and the color adjustment unit 150 correspond to the terrain change image coordinated display processing.

[0063] Next, we will supplement the explanation of the process of coordinated display of terrain change images using Figures 7 and 8. The horizontal component adjustment unit 90 shown in Figure 7 adjusts the point clouds of the old DEM (pre-earthquake) in the old DEM memory 10A and the new DEM (post-earthquake) in the new DEM memory 10B in three-dimensional space by moving the respective mesh mi in the horizontal direction (X, Y) little by little (for example, 10cm, 20cm, ... 50cm for a 1m mesh, unit: least squares method).

[0064] The horizontal component adjustment unit 90 generates a set of point clouds in this area as the pre-earthquake pattern PTai (function) and the post-earthquake pattern PTbi (function) of Pmi (S50). Since it is defined in three-dimensional space, the overall model is a crustal displacement model.

[0065] These pre-earthquake pattern PTai and post-earthquake pattern PTbi may be specifically presented in a table, for example, as shown in Figure 9.

[0066] As shown in Figure 9, the pre-earthquake pattern PTai consists of the mesh numbers mAi (i=1, 2, 3, ...) of the old DEM (pre-earthquake) and the coordinates CAi (XAi, YAi, Zai), etc., of the old point cloud dai (including 1) of mBi within the old DEM (pre-earthquake) mesh (time information may also be associated with it).

[0067] The post-earthquake pattern PTbi (i=1, 2, 3, ...: Pmi) consists of the mesh mBi number (i=1, 2, 3, ...) of the new DEM (post-earthquake) and the coordinates CBi (XBi, yBi, ZBi) within mBi (time information may also be associated with it).

[0068] Then, after adjusting the horizontal component, the terrain change image is coordinated display processing is performed (S51). This topographic change image-coordinated display process is explained using the trapezoidal embankment shown in Figures 10 and 11 as Fhi.

[0069] For example, if a trapezoidal embankment moves diagonally from west to east and a displacement occurs in the Z direction, see Figure 11. The pre-earthquake pattern PTai (i=1, 2, 3, ...) and the post-earthquake pattern PTbi shown in Figure 10(a) are generated in three-dimensional space.

[0070] However, Figure 10(a) shows the crustal displacement model (Fhi) in cross-section. ΔZi in Figure 10(a) represents the vertical ground change (corresponding to the simple elevation difference ΔZhi) before and after the earthquake in the ZX plane when viewed from direction A in Figure 11(a).

[0071] Furthermore, Figure 11 shows a trapezoidal topography, with Figure 11(a) showing the top surface (however, Figure 11(b) shows the view from direction A in Figure 11(a), and Figure 11(c) shows the view from direction B in Figure 11(a). The dotted line indicates the position before the earthquake (old), and the solid line indicates the position after the earthquake (new).

[0072] However, Figure 11 and Figure 10(a) are not identical, and Figure 11 is used to facilitate understanding of Figure 10. Figure 10(a) shows a case where the trapezoidal topography is displaced diagonally from west to east (horizontally diagonally) and the displacement occurs in the -Z direction. The section Lap between the pre-earthquake pattern PTai on the west side shown in Figure 10(a) and the post-earthquake pattern PTbi is wide when viewed from direction A in Figure 11(a) (ZX plane).

[0073] On the other hand, the section Lak between the pre-earthquake pattern PTai and the post-earthquake pattern PTbi, which are in the eastward direction, is narrow when viewed from direction A in Figure 11(a) (ZX plane). In the interval Lap, ΔZi is written as ΔZmi, and in the interval Lak, ΔZSi is written as ΔZSi.

[0074] Furthermore, in Figure 10(a), the diagonal displacement direction (three-dimensional) from the pre-earthquake pattern PTai to the post-earthquake pattern PTbi is indicated by a diagonal arrow (hereinafter, displacement direction arrow aLi).

[0075] In other words, if we illustrate Figure 10(a) three-dimensionally as shown in Figure 12, the displacement direction arrow ALi, which points from the upper western end paui (three-dimensional coordinates) of the pre-earthquake pattern PTai to the upper western end pbui (three-dimensional coordinates) of the post-earthquake pattern PTbi, can be represented by a vector defined by its length and direction. This vector is calculated by the difference image reconstruction unit 120 (shift (S32) in Figure 1) in Figure 6. This calculation will be described later.

[0076] In other words, this vector represents the displacement direction (diagonal) considering the vertical and horizontal components. In this embodiment, it is referred to as the vertical-horizontal displacement component vector VLi. By determining each pattern for each mesh, calculating the vertical-horizontal displacement component vector VLi for each of these point clouds, and taking its average, the displacement direction of this trapezoidal terrain can be determined. It is preferable to illustrate this vector on the image.

[0077] This result is stored in memory 110 (S56). Then, each time the coordinated display processing of terrain change images (including the differential processing) is completed, the change area coloring section 112 shown in Figure 6 performs gradient processing (S58).

[0078] Then, it is determined whether the Pmi is designated as an area with no change. If it is, the process returns to step S51 and the terrain change image coordinated display process (including the differential processing) is performed again (S60).

[0079] Furthermore, if it is determined in step S60 that no change area Pmi has been specified, the image is stored in memory as Gdi' (the first one being a simple difference image) and displayed on the screen (S63), as shown in Figure 8. At this time, a sub-screen is displayed on the screen to indicate whether the hue is NO or OK (not shown).

[0080] The operator then determines whether the hue of this simple difference image, Gdi', is correct. If it is incorrect, they input a NO command using the computer's mouse or other means. The change area coloring unit 112 shown in Figure 6 monitors whether a NO command or an OK command has been issued. In other words, as shown in Figure 8, it determines whether the hue is OK (S64).

[0081] Furthermore, we will provide additional explanation regarding the coordinated display processing of terrain change images. This embodiment will be explained using the front and rear (old and new) patterns shown in the ZX plane of Figure 10(a).

[0082] The topographic change image coordinated display processing defines the unchanged area Pmi (function) as having an elevation difference of ±0 (dotted line Lhi in Figure 10(a)), and approximates the pre-earthquake pattern PTai to the unchanged area Pmi (so that the difference is ±0). In this embodiment, the shift amount is calculated to perform this approximation. The shift amount will be described later.

[0083] Then, the shift amount is used to correct the elevation ZBi of the entire mesh mbi (mb1, mb2, ...) of the new DEM area Ebi (if it is positive than Pmi, ZBi is subtracted by the shift amount; if it is negative, the shift amount is added). Then, the simple elevation difference ΔZhi between the corrected new DEM and the old DEM is calculated using a discretization process, and a gradient process is performed in step S58. This difference relative to ±0 (the elevation of the unchanging area Pmi) is also called the elevation difference relative to ±0 ΔLhi (i=1, 2, ...).

[0084] The aforementioned elevation difference ΔLhi relative to ±0 indicates the degree of subsidence or uplift relative to an elevation difference of ±0 (elevation of the unchanging area Pmi: indicated as Lhi in Figure 10). Subtract the difference from the elevation difference ±0 (elevation of the unchanging area Pmi: Pmi is defined in the new DEM (post-earthquake)). If this subtracted value is negative, this negative value is calculated as the elevation difference Δhi (-Δhi) relative to the elevation difference ±0 (elevation of the unchanging area Pmi). In other words, it is settling relative to Pmi.

[0085] Furthermore, if the subtraction value is positive, this positive value is calculated as the elevation difference Δhi (+Δhi) relative to the elevation difference ±0 (elevation of the unchanging area Pmi). In other words, it is uplifted relative to Pmi. The results of this elevation difference calculation process relative to ±0 are shown in Figure 10(b). In other words, the line Lhi representing ±0 (elevation of the unchanging area Pmi) is defined in the three-dimensional coordinate system, and the calculated elevation difference Δhi is defined in this three-dimensional coordinate system (see Figure 10(b)).

[0086] Figure 10(b) is shown in the ZX plane (coordinate system). However, Figure 10(b) shows the line (dotted line Lhi) representing ±0 (elevation of the area Pmi with no change). In Figure 10(b), the shaded area below the ±0 line (dotted line Lhi), which represents the elevation of the area Pmi with no change, is called the trapezoidal settlement displacement range Ui, and the shaded area above it is called the uplift displacement range Ri.

[0087] Figures 10(a) and 10(b) show the intervals Lak, Lap, Lbp, La1, La2, Lb1, Lb2, Lb3, Lb4, and Lb6 for the old and new patterns, respectively.

[0088] To further explain using Figure 10(b), in the area where only subsidence (vertical direction (negative)) occurred, the elevation difference ΔLhi relative to ±0 (elevation of the unchanged area Pmi) is the same as the vertical ground deformation width -ΔZi.

[0089] In this embodiment, Figure 10(b) is referred to as a cross-sectional view of the vertical and horizontal component displacements with respect to elevation in the area Pmi where no change occurs. Furthermore, the interval Lak in Figure 10(a) is the sum of (-ΔZi) + (-ΔZmi) (Ui).

[0090] On the other hand, the section Lbk in Figure 10(b) is a raised section (Ri) with respect to the line (dotted line Lhi) of ±0 (elevation of the area Pmi with no change). Figure 10(b) shows ΔZmi(-) from the pre-earthquake elevation line Lai to demonstrate that the horizontal displacement range and direction can be estimated (determined) due to the diagonal displacement of the trapezoidal topography.

[0091] In other words, the -Δhi due to sedimentation is (-ΔZi) + (-ΔZmi), and the Z component of Ui is (-ΔZi) + (-ΔZmi) + (ΔZi) = -ΔZmi. In this embodiment, in Figure 10(b), the Z component of Ui is denoted as -Δhqi, and the Z component of the pre-earthquake elevation line Lai from Lhi is denoted as -Δhpi.

[0092] On the other hand, in Figure 10(b), the +Δhi due to the uplift is ΔZsi.

[0093] Section La1 in Figure 10(b) represents the displacement range in the X direction (east) on the western side of the trapezoidal topography, where -Δhqi is the vertical component (Zi) due to the displacement (oblique) on the western side of the trapezoidal topography caused by subsidence of -Δhpi (-ΔZi).

[0094] On the other hand, section La3 in Figure 10(b) represents the displacement range in the X direction (east) on the eastern side of the trapezoidal topography, where +Δhpi is the vertical component (Zi) due to the displacement (diagonal) on the eastern side of the trapezoidal topography caused by subsidence of -Δhpi (-ΔZi).

[0095] Figure 10(c) shows the colored ranges (M1, M2, ... M11). M3 represents the gradient range on the west side of the old and new sections, and M9 represents the gradient range on the east side.

[0096] The intensity of the color is determined by the depth of -Δpqi((-Δhpi)+(-Δhqi).

[0097] If the vertical and horizontal components are displaced, the original geological features (Ah1, Ah2) will be displaced as shown in Figure 13, Bh1, Bh2, before and after the earthquake.

[0098] For example, if the ground rises relative to the previous state, the resulting image is shown in Figure 14 (this is also called a simple difference image). However, Figure 14 is the opposite of Figure 10, showing the case where the ground subsided before the earthquake and rose after the earthquake.

[0099] As shown in Figure 14, if an entire area is uplifted, the entire area will be light red depending on the amount of uplift (+Z component), while areas that are similar to the old and new areas will appear whitish.

[0100] The operator determines whether the hue of this simple difference image Ghai (also called a simple difference diagram) is correct. If it is determined to be incorrect, an area Pmi is specified that is away from the local region and assumes there is no displacement (Z) (see Figure 14).

[0101] The area Pmi shown in Figure 14 is free of trees and not particularly rugged, so a larger area than Fhi is preferable.

[0102] In Figure 14, the leftmost part is a valley-like terrain, so it is preferable to make the area Pmi smaller. This coloring process (gradient) of the changed area makes the area around the local region Fhi appear whitish, as shown in Figure 15. This adjustment is performed by changing the area Pmi until the hue is deemed correct (Pm2 and Pm3 in Figure 15). If the hue is determined to be correct (hue OK), the grayscale red relief map and the shift difference map are combined (see Figure 16).

[0103] Next, color adjustment is performed. This makes the landslide area (localized region) clearly visible, as shown in Figure 17 (image is Ghg).

[0104] Figure 18 shows an example of images of an area before and after an earthquake, processed according to the above embodiment, clearly illustrating the uplift and subsidence before and after the earthquake.

[0105] Furthermore, as shown in Figure 19, local areas may be color-coded according to their degree of subsidence and uplift.

[0106] (Embodiment 2) Figure 20 is a flowchart illustrating the crustal deformation stereoscopic image processing system of Embodiment 2. This embodiment estimates crustal deformation, extracts landslides, and identifies completely destroyed houses, etc. A Digital Surface Model (DSM) is used to identify completely destroyed houses.

[0107] The computer stores the DEM (Digital Image Module) from before the earthquake (old), the DEM (Digital Image Module) from after the earthquake (new), the DSM (Digital Image Structure) from before the earthquake (old), and the DSM (Digital Image Structure) from after the earthquake (new) in its memory.

[0108] Then, as shown in Figure 20, the computer reads the pre-earthquake (old) DEM (Eai) and the post-earthquake (new) DEM (Ebi) (S100, S104). The pre-earthquake (old) DEM (Eai) and the post-earthquake (new) DEM (Ebi) are stored in memory. A simple differential processing is performed on these DEMs (S110).

[0109] The image obtained through this process is called a simple difference image (GDI).

[0110] On the other hand, in order to detect the lateral displacement vector, a simple difference operation is performed between the pre-earthquake (old) DSM (Eai) and the post-earthquake (new) DSM (Ebi) (S112). The image obtained by this operation is also called the new-old DSM simple difference image.

[0111] Furthermore, the post-earthquake (new) DEM and post-earthquake (new) DSM are combined to create a red 3D map (S104). Then, this red 3D image is converted to grayscale (S108).

[0112] On the other hand, a simple elevation difference image (a simple difference image of the old and new DEMs) and a simple difference image of the old and new DSMs are combined, and a gradient process is applied to this combined image (S106). In other words, a red, white, and blue image proportional to the elevation difference is obtained. This is displayed on the screen. At this time, a sub-screen (not shown) is displayed on the screen to input whether or not it should be shaded.

[0113] The operator then determines whether the red, white, and blue image is shaded or not. If it is determined not to be shaded, the operator inputs "not shaded" on the sub-screen. This causes the computer to shift and correct the X and Y axes of the red, white, and blue image (S114). The amount of this X,Y shift is determined as shown in equations 1 and 2 described later (calculated by setting ΔZ=0).

[0114] Furthermore, if the input indicates that the degree of blue or red in the red-white-blue mesh is incorrect, a Z-shift process is performed (S116). This Z-shift process is calculated using equations 1 and 2 described later (with ΔX and ΔY = 0). Note that ΔX is also written as dx, ΔY as dy, and ΔZ as dz.

[0115] Furthermore, if the white color of the mesh in the red-white-blue image is determined to be appropriate, a multiplicative composite image of the gray 3D map and the red-white-blue image (difference) is obtained (S124).

[0116] This allows for the interpretation and extraction of landslide blocks from the DEM (S128). It also allows for the interpretation and extraction of completely destroyed houses from the DSM (S126). The aforementioned X and Y shift correction is performed until it is determined that the image is not shaded, and the horizontal component vector at that point can be extracted (S120).

[0117] Furthermore, the aforementioned Z-shift processing is performed until red, blue, or white is deemed correct, and the Z-correction amount (Z-shift amount) at that time is extracted as the vertical variation amount (S122).

[0118] The image obtained through this process is shown in Figure 21, and a magnified view of a part of it is shown in Figure 22. Then, the vertical displacement amount from step S122 and the horizontal displacement amount (horizontal displacement vector) from step S120 are combined (S201). Then, the pre-earthquake DEM and pre-earthquake DSM are corrected using this composite vector (S204). Figure 23 shows an example of an image obtained by this process in another embodiment (moire pattern). Figure 24 shows an image when combined with an orthophoto image.

[0119] Next, the process of generating the red 3D image will be explained using Figure 25. As shown in Figure 25, a composite image of grayscale representation is obtained by multiplying and combining the ground-level opening image Dp (with ridges emphasized in white) and the underground-level opening image Dq (with bottoms emphasized in black) (Dh = Dp + Dq).

[0120] Then, an image is generated by assigning red to the gradient-enhanced image Dr according to its gradient, and a red 3D image KGi is obtained by combining this image with the composite image Dh, in which the ridge is enhanced in red.

[0121] Specifically, as shown in Figure 25, a composite image Dh with a grayscale representation is obtained by multiplying and combining the ground-level opening image Dp (with white emphasis on ridges) and the underground-level opening image Dq (with black emphasis on bottoms), and a slope-enhanced image Dr is obtained in which red is emphasized as the slope increases compared to the slope image Dra.

[0122] Then, by combining this gradient-enhanced image Dr with the composite image Dh, a red-colored stereoscopic image KGi is obtained in which the ridges are emphasized in red.

[0123] The red relief map may be colored in various ways, such as blue, green, light green, red, purple, vermilion, orange, or yellow, depending on the area and season. However, for seas, lakes, rivers, etc., it is preferable to use shades of blue or brown.

[0124] In the above embodiment, we described an example of high-speed super-resolution processing using a DEM of the ground surface, but high-speed super-resolution processing may also be performed using a DEM of the seabed.

[0125] Furthermore, the computer provides a means for obtaining a ground-level opening image (Dp) in which brighter colors are assigned to the greater the ground-level opening value, a ground-level opening image (Dq) in which darker colors are assigned to the greater the ground-level opening value, and a slope-enhanced image (Dr) in which red is emphasized as the slope value increases. A means for obtaining a first composite image (Ki) by superimposing a ground-level opening image (Dp), a ground-level opening image (Dq), and a slope-enhanced image (Dr), A means for reading image data from a ground aperture image (Dp) and obtaining a data assigned to the a* channel for each reading, A means for reading image data of the underground opening image (Dq) and obtaining b data assigned to the b* channel for each reading, A means for reading image data of a gradient-weighted image (Dr), and assigning each reading to an L* channel to obtain L data, A means of obtaining Lab image data (Li) of ground opening (Dp), underground opening (Dq), and slope-weighted image (Dr) by defining these data in the L*a*b* space each time a data, b data, and L data are obtained. Means for generating a second composite image (Lab color red image KLi) by combining a Lab image (Li) and the first composite image (Ki), It may be made to perform its function as such.

[0126] The above topographic change image enhancement process will be explained in more detail. Before performing the topographic change image enhancement process, first prepare DEMs (Digital Element Models) before and after the earthquake (e.g., topography of a peninsula) (S20, S10 in Figure 1). First, the raw difference (simple difference) is taken (S24 in Figure 1), and the operator views the simple elevation difference image Gdi, which includes "crustal deformation + landslides / collapse".

[0127] (Terrain change image enhancement processing) Then, the operator specifies several reference areas (unchanged areas Pmi (stable areas)) where landslides or collapses have not occurred. Initially, these are located near the change point Fhi (one at roughly the same elevation). In that context, "the visible difference = only crustal deformation" can be considered.

[0128] Then, in that reference region, the difference between the preceding and succeeding DEMs becomes the minimum squared error. Translation in the X and Y directions (horizontal component) A fixed offset in the Z direction (vertical component) Find the optimal dx, dy, and dz. Then, we shift all the "later DEMs" by the calculated (dx, dy, dz) and take the difference again to obtain the simple elevation difference image G'di.

[0129] This cancels out crustal movements, leaving only localized sediment movement (collapse and deposition) visible. This allows the following processes to be performed. Limiting to the area of ​​soil movement (e.g., the polygon mask of a collapsed slope), the volume (m³) 3 Calculate ). When the amount of variation differs depending on the region, calculations are performed separately for each region.

[0130] As a premise: □Pre-earthquake DEM:pre_dem.tif □Post-earthquake DEM:post_dem.tif □(Optional) Sediment Movement Area Mask (1: Target, 0: Other): landslide_mask.tif *If it's not there, just comment out the final volume calculation section. □The reference area is specified as a "rectangle of pixel coordinates". Example: Provide a list of (row_min, row_max, col_min, col_max)

[0131] Summary of key points for use Data preparation Place opre_dem.tif and post_dem.tif in / content / in Colab. (Mounting Google Drive is also OK) If necessary, create a landslide_mask.tif file from the collapsed slope polygon and adjust it to the same resolution and range. How to determine the reference area Looking at the rough_diff diagram in the O code, □ "No slope collapse or soil movement has occurred." □"But crustal movements are happening." Specify several areas using row and column indices. A rough rectangle is fine to start with. You can change it as many times as you like if you don't like it. Search range for shift amount I've set omax_shift_xy to ±5, but If the expected horizontal displacement is large, expand the range to ±10, ±20, etc. The wider the search range, the more computationally intensive it becomes.

[0132] When both vertical and horizontal are large In this implementation, XY is calculated by brute-force searching for integer pixels, and Z is calculated by analytically finding the optimal value for each shift. We optimize it by bundling it together in this way. You can also do it step by step, like "aligning Z first, then X and Y," Ultimately, what we want is to minimize the RMSE (dx, dy, dz), Essentially, it's the same idea.

[0133] Cancellation of crustal deformation components By applying the optimized (dx, dy, dz) in the reference region to the entire terrain in one go, We are removing the component that can be considered almost as rigid body displacement. As a result, the differential DEM (ideally) only retains localized collapse and deposition.

[0134] In other words, the terrain change image coordinated display processing, which is step S32 (Z shift) shown in Figure 1 and the difference image reconstruction unit 120 shown in Figure 6, performs the following processing more specifically. This process is the core of the present invention, and it begins by ignoring areas of sediment movement. Focusing only on the unchanged area Pmi (stable area SC) where landslides and other sediment movements are judged not to have occurred, the pure crustal deformation component is estimated by 3D transforming the entire DEM (old DEM) so that the elevation difference before and after the earthquake is minimized in this unchanged area Pmi (stable area SC) (see Figure 26).

[0135] During a major earthquake, simple DEM differences (simple differences) cannot separate wide-area seismic deformation (uplift, subsidence, and tilting) from localized sediment movement (landslides and collapses) (see Figure 27), making it extremely difficult to accurately grasp the location and amount of sediment movement. This embodiment models the deformation as a 3D Affin transformation combining crustal deformation, tilt, and height offset. It then finds (calculates) five parameters (α) that best represent this deformation.

[0136]

number

[0137] In other words, the parameter (α) that minimizes the error within the unchanging area Pmi (stable region) is the parameter that minimizes the error within the stable region. ★ Search for ). Optimal shift amount conversion parameter (α ★ ) is the transformed pre-earthquake DEM (Z) within the stable region (SC). - pre') and post-earthquake DEM(Z - The objective function (E(α)) is defined as minimizing the sum of squares of the elevation differences (pre) (see Figure 28).

[0138]

number

[0139] Here, the optimization order for vertically dominant types (also known as the Noto type, for example) and horizontally dominant types (also known as the Kumamoto type) will be explained using Figures 29 (vertically dominant type) and 30 (horizontally dominant type). In the case of a vertically dominant type, step 1 involves optimizing the height and slope components (Figure 29(a)). In this case, we fix ΔX=ΔY=0 and optimize the error E for (ax, ay, a0) using the least squares method. Then, in step 2, re-optimize all parameters (Fig. 29(b)). That is, with the result of step 1 as the initial value, minimize the error E for all parameters α.

[0140] On the other hand, in the case of the horizontal-dominant type (Fig. 30), in step 1, estimate the initial value of the horizontal movement. Estimate the initial values of (ΔX, ΔY) using the cross-correlation method or feature point matching within the non-changing area Pmi(SC). Then, in step 2, perform the overall optimization of all parameters. With the result of step 1 as the initial value, minimize the error for all parameters α.

[0141] That is, since the parameter stability is different between the vertical-dominant type and the horizontal-dominant type, as shown in Fig. 31, judge the type of crustal movement (vertical-dominant type or horizontal-dominant type) in advance (d100), and switch the order of converging the parameters (d200 or d300), so as to derive a more robust and highly accurate solution.

[0142] That is, as shown in Fig. 32, prepare the data (S1). Read the DEMs (Z - pre, Z - post) before and after the earthquake, and unify the coordinate system and mesh size. Then, set the stable region SC(Pmi) (S2). This setting excludes the known land movement amounts from the rough differences and aerial photo interpretation. Next, judge the type of variation (S3). For example, judge whether it is the Noto type or the Kumamoto type from GNSS data or the like.

[0143] Then, estimate the optimal parameters (α ★ ). That is, execute the optimization sequence according to the type of variation. Next, create the corrected DEM (Z - pre ★ ). That is, use the optimal parameters (α ★ ) to convert Z - pre over the entire area. Then, the simple difference DEM(ΔZ) is calculated (S6).

[0144] This calculation is ΔZ=Z - post-(Z - pre ★ ) It is calculated as follows. Next, sediment transport is extracted and quantified (S7). A threshold is applied to ΔZ to extract the transport area. This allows for the calculation of erosion and depositional values.

[0145] This allows us to extract true topographic changes from the noise (background) (see Figures 33 and 34). Figure 33 shows an image obtained using a conventional simple difference method, while Figure 34 shows an image of the difference ΔZ after deformation correction using the method of this embodiment. In other words, the method of this embodiment effectively removes wide-area crustal deformation components, and only the sediment transport signal is clearly extracted. Table 1 shows the sediment balance (erosion and deposition) for areas 1-3, as shown in Figure 35.

[0146] [Table 1]

[0147] Furthermore, this will be explained using Figure 36. Figure 36 is an explanatory diagram illustrating the process of canceling crustal deformation by moving the unchanging area Pmi (SC or reference area). In Figure 36(a), the dotted mesh represents the unchanging area Pmi (SC or reference region) (Poset-DEM), and the thick mesh represents the old DEM (pre-DEM). The old DEM (pre-DEM) is the designated area Eai and is a large area, but in Figure 36(a), it is the same size as the unchanging area Pmi (Poset-DEM).

[0148] Then, as shown in Figure 36(b), the optimal (dx, dy, dz) is searched for while sequentially defining this unchanging area Pmi (SC or reference region) in the old DEM (pre-DEM) and minimizing the RMSE (Root Mean Squared Error) of the old DEM (pre-DEM) with respect to the unchanging area Pmi (SC or reference region).

[0149] For example, define the central mesh of the unchanged area Pmi as the mesh of the old DEM (for example, starting from the leftmost mesh), and move the area corresponding to the old DEM Pmi one mesh at a time (this area is moved in the XY direction first, then in the Z direction). Figure 36(c) shows the relationship between RMSE and shift amount. The circle indicates the convergence point.

[0150] Then, the found (dx, dy, dz: vectors) are applied to the entire mesh of the new DEM (post-earthquake) to generate a "DEM corrected for crustal deformation" (Figure 36(d)). The old DEM can also be used. The search for the optimal shift amounts (dx, dy, dz) is specifically described in Figure 37.

[0151] By performing this process and taking the simple difference between the old and new DEMs, we obtain the aligned image shown in Figure 38. In other words, the area around Fhi becomes white, and Fhi stands out. Note that the aligned image in Figure 38 also displays contour lines. However, Fhi is a different example from that in Figure 2.

[0152] Next, the process of displaying coordinated terrain change images will be explained using the flowchart in Figure 39. However, step S30 in Figure 1 (defining the unchanged area Pmi(SC) in the simple difference image Gdi) will be shown. This Pmi(Pm1) is placed near Fhi by the operator.

[0153] Then, the difference image reconstruction unit 120 calculates the shift amount (α ★A screen (not shown) is displayed to prompt the user to input an initial value a0 for calculating the value (S324). This initial value a0 (also called height offset) is the initial value in the Z direction and is determined and input according to the color of the gradient of elevation (red) or subsidence (blue) in the simple difference image Gdi in Figure 2(d).

[0154] For example, if the area is generally uplifted, enter +10m or +9m, ... or 1m, .... If the area is generally subsiding, enter -10m or -9m, ... -1m, ... (start with the largest value). Then, the unchanged area Pmi(SC) is defined in the old DEM (S328).

[0155] Then, the parameter (α) minimizes the error between the model (function) of the area PAi in the old DEM corresponding to Pmi as defined in this old DEM and the model (function) of Pmi (see Figure 36(c)). ★ We search for (α) (using equations 1 and 2) (S330). ★ ) corresponds to the shift amount.

[0156] Then, it is determined whether convergence has occurred (S332: see Figure 36(d)). If it has not converged, (α ★ ) is divided into two parts using the bisection method, and this is updated (α ★ ) and return the process to step S330 (S334). Also, if it is determined that convergence has occurred in step S332, update (α ★ ) Corrects all meshes of the new DEM in area Ebi (S336).

[0157] Then, in S24 of Figure 1, a simple elevation difference process is performed (S338). The aforementioned convergence is preferably determined when the length becomes, for example, 0.5m or less (0.4m is also acceptable). Figure 40 also shows the specification of the reference region and the search for the optimal shift.

[0158] Although the explanation described approximating the old DEM to Pmi, it is also acceptable to approximate Pmi to the old DEM. [Explanation of Symbols]

[0159] 10. Computer main unit 10B New DEM Memory 10A Memory for older DEMs 40 Red stereoscopic mapping unit 40, 60 Grayscale conversion section 90 Horizontal adjustment section 100 Differentiation section 130 Multiplication and Composition Section

Claims

1. (A) A storage unit that stores the old DEM from a predetermined time and the new DEM acquired at a time more recent than the predetermined time, (B) Means for matching the horizontal components of the point cloud of the old DEM and the point cloud of the new DEM, (C) A means for determining the simple elevation difference between the old DEM and the new DEM after matching the horizontal components, (D) Means for generating and displaying a simple elevation difference image on a screen, the first color indicating no difference if the simple elevation difference value is within a predetermined value, the second color gradient indicating uplift if it is above the predetermined value, and the third color gradient indicating subsidence if it is below the predetermined value. (E) Means for generating a red stereoscopic image based on the new DEM, (F) Means for generating a red grayscale image by converting the red 3D image to grayscale, (G) Means for generating a crustal displacement stereoscopic image by multiplying and combining the simple elevation difference image and the red grayscale image and displaying it on a screen, A crustal deformation stereoscopic image processing system characterized by having the following features.

2. moreover, (H) A means for calculating a model of the unchanging area (Pmi) based on the mesh within the unchanging area (Pmi) set in the simple elevation difference image, (I) A means for calculating a model of an area (PAi) based on the mesh within the area (PAi) of the old DEM that corresponds to the area (Pmi) without change, (J) A means for determining the minimum error vector when the model of the old DEM area (PAi) is approximated within the predetermined value and converges to the model of the change area (Pmi), (K) A means for correcting the coordinates of all meshes of the new DEM with the minimum error vector, recalculating the simple elevation difference value between the old DEM and the new DEM, and generating a simple elevation difference image (G'di) in which the first color component is more numerous, thereby emphasizing the changed areas. A crustal deformation stereoscopic image processing system according to claim 1, characterized by having the following features.

3. The crustal deformation stereoscopic image processing system according to claim 1 or 2, characterized in that the first color is "white", the second color is "red", and the third color is "blue".

4. The crustal deformation stereoscopic image processing system according to claim 1, characterized in that the area of ​​no change (Pmi) is larger than the area of ​​change (Fhi) in the simple elevation difference image, and the means (H) to (K) are executed again each time this area is set to a different area from the previous area of ​​no change (Pmi).

5. The storage unit includes an old DEM from a predetermined time, a new DEM acquired at a later time than the predetermined time, an old DSM from the same region as the old DEM from the predetermined time, and a new DSM. (L) Means for matching the horizontal components of the old DSM and the new DSM, (M) A means for determining the simple elevation difference value of the DSM between the old DSM and the new DSM after matching the horizontal components, (N) Means for generating a DSM elevation difference image colored based on the DSM elevation difference values, (O) Means for generating a DSM red stereoscopic image based on the new DSM, (P) A means for generating a DSM red grayscale image by converting a DSM red 3D image to grayscale, (Q) A means for generating a ground surface displacement stereoscopic image by multiplying and combining the DEM elevation difference image, the DEM red grayscale image, and the DSM red grayscale image, A crustal deformation stereoscopic image processing system according to claim 1, characterized by having the following features.

6. The crustal deformation stereoscopic image processing system according to claim 1, characterized in that, if the horizontal component is shaded, the X and Y directions are shifted by a predetermined amount.

7. Computers, (A) Means for storing in the memory unit old DEM from a predetermined time and new DEM acquired at a time more recent than the predetermined time, (B) Means for matching the horizontal components of the point cloud of the old DEM and the point cloud of the new DEM, (C) Means for determining the simple elevation difference between the old DEM and the new DEM after matching the horizontal components, (D) Means for generating and displaying a simple elevation difference image on a screen, the first color indicating no difference if the simple elevation difference value is within a predetermined value, the second color gradient indicating uplift if it is above the predetermined value, and the third color gradient indicating subsidence if it is below the predetermined value. (E) Means for generating a red stereoscopic image based on the new DEM, (F) Means for generating a red grayscale image by converting the red 3D image to grayscale, (G) Means for generating a crustal displacement stereoscopic image by multiplying and combining the simple elevation difference image and the red grayscale image, and displaying it on a screen. A crustal deformation stereoscopic image processing program that performs the function of [this].

8. Furthermore, the computer, (H). Means for calculating a model of the unchanging area (Pmi) based on the mesh within the unchanging area (Pmi) set in the simple elevation difference image, (I) Means for calculating a model of an area (PAi) based on the mesh within the area (PAi) of the old DEM that corresponds to the area (Pmi) without change, (J) A means for determining the minimum error vector when the model of the old DEM area (PAi) is approximated within the predetermined value and converges to the model of the change area (Pmi), (K) A means for correcting the coordinates of all meshes of the new DEM with the minimum error vector, recalculating the simple elevation difference value between the old DEM and the new DEM, and generating a simple elevation difference image (G'di) in which the first color component is more numerous, thereby emphasizing the areas of change. A crustal deformation stereoscopic image processing program according to claim 7, which performs the function of being a crustal deformation stereoscopic image processing program.

9. The crustal deformation stereoscopic image processing program according to claim 7 or 8, characterized in that the first color is "white", the second color is "red", and the third color is "blue".

10. Computers, The crustal deformation stereoscopic image processing program according to claim 7, characterized in that the area of ​​no change (Pmi) is larger than the area of ​​change (Fhi) in the simple elevation difference image, and the means (H) to (K) are executed again each time this area is set to a different area from the previous area of ​​no change (Pmi).

11. The storage unit includes an old DEM from a predetermined time, a new DEM acquired at a later time than the predetermined time, an old DSM from the same region as the old DEM from the predetermined time, and a new DSM. Computers, (L). Means for matching the horizontal components of the old DSM and the new DSM, (M) Means for determining the simple elevation difference value of the old DSM and the new DSM after matching the horizontal components. (N). Means for generating a DSM elevation difference image colored based on the DSM elevation difference values. (O) Means for generating a DSM red stereoscopic image based on the new DSM, (P) Means for generating a DSM red grayscale image by converting a DSM red 3D image to grayscale, (Q) Means for generating a ground surface displacement stereoscopic image by multiplying and combining the DEM elevation difference image, the DEM red grayscale image, and the DSM red grayscale image, A crustal deformation stereoscopic image processing program according to claim 7, which performs the function of being a crustal deformation stereoscopic image processing program.

12. Computers, The crustal deformation stereoscopic image processing program according to claim 7, which, if the horizontal component is shaded, performs the function of a means to shift the X and Y directions by a predetermined amount.

Images