Three-dimensional space-time kernel density estimation method based on adaptive bandwidth

Abstract

The invention discloses a three-dimensional space-time kernel density estimation method based on adaptive bandwidth, and the method comprises the steps: determining the boundary contour of a research region in a geographic space, determining the lower bound and upper bound of the occurrence time of a target event, constructing an effective space-time cube, and collecting a target event sample in the interior, thereby obtaining a target event data set. And inputting the data set to calculate and fit the real background distribution of the target event, namely a joint probability density model. And carrying out time-by-time slicing processing on the joint probability density model to obtain a space condition probability density model. And respectively carrying out qualitative analysis and quantitative analysis on the two models according to requirements, and analyzing a spatial distribution rule and a time evolution characteristic of a high-incidence area (hot area) of a target event. The adaptive bandwidth change is introduced into the traditional algorithm, so that the algorithm provided by the invention has better fitting result, more accurate prediction performance, more concentrated hot area identification capability and more robust extreme event impact absorption capability.

Description

Technical Field

[0001] This invention belongs to the field of three-dimensional spatiotemporal kernel density estimation technology, specifically relating to a three-dimensional spatiotemporal kernel density estimation method based on adaptive bandwidth. Background Technology

[0002] In existing technologies, the traditional Spatiotemporal Kernel Density Estimation (TSTKDE) method is a commonly used and mature approach for analyzing and evaluating geographic information events (or points of interest). It uses a limited number of event records to determine the optimal fixed bandwidth on both the spatial and temporal sides, and then fits the true probability density distribution behind the event using kernel methods. This helps experts, scholars, and institutions explore and analyze the historical development and changes of these events, classify risk levels, provide key decision-making basis for downstream tasks, and predict future development trends. However, this traditional algorithm has drawbacks. It always uses a fixed bandwidth for model fitting, which inevitably leads to oversmoothing estimates in densely distributed data areas, resulting in lost details and blurred hotspots, while in sparsely distributed data areas, it produces undersmoothing estimates with noise, ultimately leading to unreliable fitting results. This limitation stems from the fixed bandwidth selection of traditional algorithms, inevitably arising when processing non-uniform data, regardless of whether the application is one-dimensional, two-dimensional, or three-dimensional, or whether the chosen fixed bandwidth is the "optimal bandwidth." Real-world events are often non-uniformly distributed in space and time. Some international experts have recognized this and proposed that 3D TSTKDE should innovate towards adaptive bandwidth to better address the biases caused by traditional algorithms in handling non-uniform problems. However, research on spatiotemporal kernel density estimation methods with adaptive bandwidth is very limited both domestically and internationally. Furthermore, these methods are computationally intensive in 3D applications, making them difficult to implement practically, lacking engineering feasibility, and lacking concrete real-world data application examples for empirical verification.

[0003] Therefore, a new three-dimensional spatiotemporal kernel density estimation method is needed to solve the above-mentioned technical problems. Summary of the Invention

[0004] The purpose of this invention is to provide a three-dimensional spatiotemporal kernel density estimation method based on adaptive bandwidths (STKDE with adaptive bandwidths, hereinafter referred to as ASTKDE) to solve the problem that traditional spatiotemporal kernel density estimation methods in the prior art produce unreliable fitting results when processing non-uniformly distributed data.

[0005] The technical solution of this invention to solve its technical problem is as follows:

[0006] A three-dimensional spatiotemporal kernel density estimation method based on adaptive bandwidth includes the following steps:

[0007] S1: Determine the geographical boundary of the study area and the lower and upper time bounds of the target events. Based on the boundary, upper, and lower time bounds, obtain the effective spatiotemporal domain. Collect target event samples within this domain to obtain the target event dataset. Fit the true background distribution of the target events to the target event dataset, thus obtaining the joint probability density model. The specific joint probability density function is:

[0008] ;

[0009] ;

[0010] ;

[0011] ;

[0012] ;

[0013] ;

[0014] ;

[0015] in, The joint probability density function is used; the boundary contour is a closed polygon on the XY plane of geographic space. The upper bound of time is The upper bound of time is . (By a closed polygon) and upper bound of time, lower bound of time The enclosed area is the effective spatiotemporal domain. , composed of closed polygons Minimum bounding rectangle and upper and lower time bounds The enclosed cube is a spacetime cube. In the effective spatiotemporal domain Internal collection Event Samples Each event sample That is, each event sample contains geospatial coordinates. , and time coordinates , For the first The horizontal coordinate of each event sample in geographic space; For the first The vertical coordinate of each event sample in geographic space; For the first The coordinates of an event sample in time and space; It is an indicator function, for each event sample. The contributed kernel density falls into the effective spatiotemporal domain Take 1 if the condition is met, otherwise take 0; , For each event sample The kernel function based on adaptive bandwidth is specifically the Gaussian kernel function. In This can be either in the general formula or It could also be This uses the assumption of XY isotropy, that is, any event After it occurs, the impact is consistent in all directions across space. For the general formula ; For the first Spatial adaptive bandwidth for each event sample For the first Time-side adaptive bandwidth for each event sample; for The optimal fixed bandwidth on the spatial side for each event sample; for Optimal fixed bandwidth on the time side for each event sample; For global space-side bandwidth; For global time-side bandwidth, let's directly set... , ; for The leading spatiotemporal kernel density estimation model corresponding to each target event sample The geometric mean of the kernel density values. It uses a fixed optimal bandwidth and The leading spatiotemporal kernel density estimation model at sample points The kernel density value at that location; This is the edge correction factor;

[0016] S2: Time-by-time slicing is performed on the joint probability density model to obtain the spatial conditional probability density model. The specific spatial conditional probability density function is as follows:

[0017] ;

[0018] in, It is the spatial conditional probability density function;

[0019] S3: Perform qualitative and quantitative analysis on the joint probability density model and the spatial conditional probability density model according to the requirements.

[0020] Preferably, the voxel differentiation rules of the joint probability density model include: the ASTKDE-Large model and the ASTKDE-Base model.

[0021] Preferably, the differential voxel parameters corresponding to the ASTKDE-Large model are: X = 695, Y = 735, T = 365, and the corresponding voxel spatiotemporal scale is 50 m × 50 m × 10 days; the differential voxel parameters corresponding to the ASTKDE-Base model are: X = 348, Y = 368, T = 183, and the corresponding voxel spatiotemporal scale is 100 m × 100 m × 20 days.

[0022] Preferably, the joint probability density model employs concurrent computing and load balancing during the calculation process.

[0023] This concurrent computing and load balancing strategy not only utilizes the inherently divide-and-conquer nature of the kernel density estimation (KDE) algorithm, but also stimulates computer hardware performance, greatly improves computing power, and significantly shortens computing time. It is very important in high-dimensional (three-dimensional) spatiotemporal kernel density estimation calculations, ensuring the engineering feasibility of the algorithm.

[0024] Preferably, the specific process of performing qualitative and quantitative analysis on the joint probability density model in step S3 is as follows: obtaining statistical significance. The local density peaks of the region are obtained by importing the joint probability density model into the volume rendering software. The volume rendering gradient shading is performed based on the kernel density value of each voxel. The obtained local density peaks are also imported into the scene. The probability density hotspots in the model are analyzed by displaying the volume rendering and these peaks together.

[0025] Preferably, the volume rendering gradient coloring rule is: red (salience) → Yellow (significance) → Blue (significance) → Completely transparent, with a gradient transition coloring based on the salience thresholds corresponding to red, yellow, and blue.

[0026] Preferably, the specific process of performing qualitative and quantitative analysis on the joint probability density model in step S3 is as follows: In RStudio software, the joint probability density model is plotted as a three-dimensional isosurface and cut by year. The model segments of each year are projected onto the XY plane to form the maximum projection range. The local density peak points in the year are also projected onto the XY plane, and the area of ​​the maximum projection range is calculated. The specific results of the annual slices and the corresponding statistical situation of the hot zone area with different significance are analyzed.

[0027] Preferred color: red (significance) ), yellow (significance) ), blue (significance) Coloring is done according to the significance thresholds corresponding to red, yellow, and blue.

[0028] Preferably, the specific process of performing qualitative and quantitative analysis of the spatial conditional probability density in step S3 is as follows: calculating the statistical significance of each time slice using the spatial conditional probability density model. Local density peak points within the range are connected across adjacent time slices according to a set minimum distance threshold to form the migration trajectory of the hot zone over time, and the migration trajectory of the hot zone over time is analyzed.

[0029] Preferably, the optimal fixed bandwidth on the space side Optimal fixed bandwidth on the time side It can be obtained through leave-one-out cross-validation, interpolation, log-likelihood, or empirical rules.

[0030] The beneficial effects of this invention are as follows: By determining the geographical boundary of the study area and the lower and upper time bounds of the target event occurrence, an effective spatiotemporal cube is constructed. Target event samples are collected within this cube to obtain a target event dataset. This dataset is then used to calculate and fit the true background distribution of the target events, i.e., the joint probability density model. Further time-slicing of the joint probability density model yields a spatial conditional probability density model. Qualitative and quantitative analyses are performed on these two models according to requirements to analyze the spatial distribution patterns and temporal evolution characteristics of high-incidence areas (hotspots) of target events. By introducing adaptive bandwidth variation into traditional algorithms, the invented algorithm achieves better fitting results, more accurate prediction performance, more concentrated hotspot identification capabilities, and a more robust ability to absorb the impact of extreme events. Attached Figure Description

[0031] Figure 1 This is a flowchart of the three-dimensional spatiotemporal kernel density estimation method based on adaptive bandwidth according to the present invention;

[0032] Figure 2 This is an architecture diagram of the ASTKDE algorithm of this invention;

[0033] Figure 3 This is a comparison chart of the ISE results of TSTKDE and ASTKDE in the high-density upper tail and even global regions of the synthesized data in Experiment 1 of this invention;

[0034] Figure 4This is a cube diagram of the fixed threshold τ stability of TSTKDE and ASTKDE before and after absorbing the impact of extreme rainstorm events in Experiment 2 of this invention;

[0035] Figure 5 This is the design diagram of Experiment 3 of the present invention regarding the 13 prediction schemes for data that are not available in the near future;

[0036] Figure 6 This is a graph showing the prediction results of TSTKDE and ASTKDE for the future short term without available data in Experiment 3 of this invention;

[0037] Figure 7 This is a comparison chart of the hot zone identification results of TSTKDE and ASTKDE in Experiment 4 of this invention;

[0038] Figure 8 This is a graph showing the results of the sensitivity analysis of ASTKDE to changes in input parameters in Experiment 5 of this invention;

[0039] Figure 9 This is a volume rendering of the ASTKDE spatial conditional probability density model of road collapse disasters in Zhengzhou City from 2014 to 2023, as described in Embodiment 1 of the present invention.

[0040] Figure 10 This is a plot of the hot zone of the ASTKDE joint probability density model for road collapse disasters in Zhengzhou City from 2014 to 2023, showing the area statistics of each statistically significant region. (This is from Embodiment 1 of the present invention.)

[0041] Figure 11 These are perspective views (left) and top views (right) of the hot zone (α<0.01) trajectory of the ASTKDE spatial conditional probability density model for road collapse disasters in Zhengzhou City from 2014 to 2023, as described in Embodiment 1 of the present invention.

[0042] Figure 12 This is an overlay analysis diagram of historical images of the built-up area of ​​Zhengzhou (1990-2020) and the ASTKDE spatial conditional probability density model of road collapse disaster hotspots in Embodiment 1 of the present invention;

[0043] Figure 13 This is an overlay analysis diagram of the ASTKDE spatial conditional probability density model regarding the distribution of special soil types and road collapse disasters in Zhengzhou City, as described in Embodiment 1 of the present invention. Detailed Implementation

[0044] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention.

[0045] like Figure 1 , Figure 2As shown, this invention discloses a three-dimensional spatiotemporal kernel density estimation method based on adaptive bandwidth, comprising the following steps:

[0046] S1: Determine the geographical boundary of the study area and the lower and upper time bounds of the target events. Based on the boundary, upper, and lower time bounds, obtain the effective spatiotemporal domain. Collect target event samples within this domain to obtain the target event dataset. Fit the true background distribution of the target events to the target event dataset, thus obtaining the joint probability density model. The specific joint probability density function is:

[0047] ;

[0048] ;

[0049] ;

[0050] ;

[0051] ;

[0052] ;

[0053] ;

[0054] in, The joint probability density function is used; the boundary contour is a closed polygon on the XY plane of geographic space. The upper bound of time is The upper bound of time is . (By a closed polygon) and upper bound of time, lower bound of time The enclosed area is the effective spatiotemporal domain. , composed of closed polygons Minimum bounding rectangle and upper and lower time bounds The enclosed cube is a spacetime cube. In the effective spatiotemporal domain Internal collection Event Samples Each event sample That is, each event sample contains geospatial coordinates. , and time coordinates , For the first The horizontal coordinate of each event sample in geographic space; For the first The vertical coordinate of each event sample in geographic space; For the first The coordinates of an event sample in time and space; It is an indicator function, for each event sample. The contributed kernel density falls into the effective spatiotemporal domain Take 1 if the condition is met, otherwise take 0; , For each event sample The kernel function based on adaptive bandwidth is specifically the Gaussian kernel function. In This can be either in the general formula or It could also be This uses the assumption of XY isotropy, that is, any event After it occurs, the impact is consistent in all directions across space. For the general formula ; For the first Spatial adaptive bandwidth for each event sample For the first Time-side adaptive bandwidth for each event sample; for The optimal fixed bandwidth on the spatial side for each event sample; for Optimal fixed bandwidth on the temporal side and optimal fixed bandwidth on the spatial side for each event sample. Optimal fixed bandwidth on the time side It can be obtained through leave-one-out cross-validation (LOOCV), interpolation (Plug-in), log-likelihood (LIK), and empirical rules (Silverman's rule of thumb or other empirical rules). For global space-side bandwidth; For global time-side bandwidth, let's directly set... , ; for The leading spatiotemporal kernel density estimation model corresponding to each target event sample The geometric mean of the kernel density values. It uses a fixed optimal bandwidth and The leading spatiotemporal kernel density estimation model at sample points The kernel density value at that location; This is the edge correction factor;

[0055] The voxel differential rules of the joint probability density model include: the ASTKDE-Large model and the ASTKDE-Base model. The differential voxel parameters corresponding to the ASTKDE-Large model are: X = 695, Y = 735, T = 365, and the corresponding voxel spatiotemporal scale is 50 m × 50 m × 10 days. The differential voxel parameters corresponding to the ASTKDE-Base model are: X = 348, Y = 368, T = 183, and the corresponding voxel spatiotemporal scale is 100 m × 100 m × 20 days.

[0056] The joint probability density model employs concurrent computation and load balancing during the calculation process. Specifically, the dataset is evenly divided into multiple subsets, meaning the computational task is evenly divided into multiple sub-parts, ensuring that the computational load of each sub-part is comparable, thus achieving load balancing. Then, one thread is created for each sub-part, and multiple threads perform parallel computation on the CPU simultaneously. Finally, the multiple computational results are assembled into the final result. This concurrent computation and load balancing strategy leverages the inherent divide-and-conquer nature of the kernel density estimation (KDE) algorithm while also maximizing computer hardware performance, significantly improving computing power and reducing computation time. This is crucial in high-dimensional (three-dimensional) spatiotemporal kernel density estimation calculations and is one of the important means to ensure the algorithm's practical engineering feasibility.

[0057] S2: Time-by-time slicing is performed on the joint probability density model to obtain the spatial conditional probability density model. The specific spatial conditional probability density function is as follows:

[0058] ;

[0059] in, It is the spatial conditional probability density function;

[0060] S3: Perform qualitative and quantitative analysis on the joint probability density model and the spatial conditional probability density model according to the requirements.

[0061] The specific process of qualitative and quantitative analysis of the joint probability density model is as follows: obtaining statistical significance. The local density peaks in the region are identified. The joint probability density model is imported into the volume rendering software. Volume rendering gradient shading is applied based on the kernel density value of each voxel. The obtained local density peaks are also imported into the scene. Through volume rendering and joint display of these peaks, the probability density hotspots in the model are analyzed. The corresponding volume rendering gradient shading rule is: red (significance). → Yellow (significance) → Blue (significance) → Completely transparent, with a gradual transition coloring based on the significance thresholds corresponding to red, yellow, and blue. Alternatively, by plotting a 3D isosurface of the joint probability density model in RStudio software and cutting it by year, projecting the model fragments of each year onto the XY plane to form the maximum projection range, and also projecting the local density peaks within the year onto the XY plane, and calculating the area of ​​the maximum projection range, the specific results of the yearly slices and the corresponding statistical situation of the hot zone areas with different significance are analyzed. The corresponding isosurface coloring rule is: The isosurface coloring rule is: Red (significance) ), yellow (significance) ), blue (significance) Coloring is done according to the significance thresholds corresponding to red, yellow, and blue.

[0062] The specific process of qualitative and quantitative analysis of spatial conditional probability density is as follows: The statistical significance of each time slice is calculated using the spatial conditional probability density model. Local density peaks within the range are connected across adjacent time slices based on a set minimum distance threshold to form the migration trajectory of the hot zone over time. This trajectory is then analyzed. In fact, these hot zone trajectories are not isolated but interconnected, exhibiting clustering, splitting, and disappearance on a large scale. These collective behaviors reveal the temporal evolution characteristics of the event's hot zone. Further analysis of these hot zone trajectories is necessary.

[0063] This application innovatively introduces a point-by-point varying, data-driven adaptive bandwidth to fit a joint probability density model. Compared to traditional algorithms, this model outperforms them in several ways: 1) Whether in high-density tail regions or across the entire domain, the joint probability density model exhibits a smaller ISE (integral squared error), resulting in a better fit. This advantage is more pronounced with increasingly uneven data distribution and is applicable to various optimal bandwidth selection strategies. 2) Under extreme event impacts, the joint probability density model is more sensitive and focused. It can migrate more hotspots into the impact window, focusing more "attention," and retain more hotspot skeletons in historical windows, resulting in fewer "forgotten" hotspots and generating more accurate predictions within the impact window—something traditional algorithms cannot achieve. 3) The joint probability density model outperforms traditional algorithms in predicting future data that is not yet available in the short term. It does not overfit and provides more reliable future predictions. 4) By assigning smaller bandwidth to dense regions, more details of the target event can be characterized, facilitating downstream hotspot analysis tasks—a deficiency in traditional algorithms. In addition, through 1) model variants (voxel differential rules), two model variants with different numbers of voxels (parameters) are given, which can facilitate different users to select according to their own hardware computing power and application scenarios; 2) concurrent computing and load balancing significantly improve computing power and effectively solve the problem of large computing volume and large computing time required in high-dimensional (here, three-dimensional) scenarios; 3) XY isotropy can be abstracted into XY anisotropy. At the code implementation level, it can not only support XY isotropy (XY isotropy is often assumed in geospatial analysis), but also XY anisotropy. This makes the innovative algorithm of this invention applicable not only to geospatial and temporal events, but also to events in all other research fields that can be abstracted into three-dimensional points.

[0064] The following four experiments will illustrate the superior performance and scientific validity of the "Adaptive Bandwidth-Based 3D Spatiotemporal Kernel Density Estimation Algorithm (hereinafter referred to as ASTKDE)" compared to the traditional algorithm (based on fixed optimal full-band spatiotemporal kernel density estimation, hereinafter referred to as TSTKDE) from different perspectives.

[0065] Experiment 1: Comparison and evaluation of fitting effect. The ISE (integral squared error) performance of ASTKDE compared with TSTKDE in high-density tail (local, or top n% voxel) regions and even global (100% voxels) was evaluated and verified by “simulating the background distribution of events in an urban environment with artificially synthesized data”.

[0066] A generalized probability density function (pdf) was designed to simulate the distribution of real-world events in an urban environment. By assigning different specific parameters to this generalized pdf, three instantiated background distribution probability density functions (pdf1, pdf2, and pdf3) were generated. These three pdfs exhibit increasing "non-uniformity," meaning the original distribution of the data becomes increasingly uneven. The formula for this generalized pdf is shown in Table 1, and it is expressed in two parts: "temporal distribution" and "spatial distribution."

[0067] Table 1 Generalized PDF Formula

[0068]

[0069] As shown in the table above, among which,

[0070] 1) On the time side, a three-peak superposition method is used. Weighted To simulate the periodic changes in the intensity of events over time; a linearly increasing function. Weighted This is used to simulate the escalating and deteriorating trend in the urban environment over time, where 'a' is the slope and 'b' is the offset, and it needs to satisfy... The system uses three Gaussian distributions superimposed to simulate the temporal impact of the alternation between flood season and non-flood season (the distribution of different rainfall amounts or other influencing factors over time); and a linearly increasing function is used to simulate the impact of irreversible environmental degradation in the urban environment over time (e.g., gradual management disorder, infrastructure damage and aging with increasing service life, or other factors that are constantly aggravating the deterioration) on the event. The weights for each of the three peaks must also satisfy the condition that the sum of the weights is 1. , This indicates that each peak follows a Gaussian distribution, using... express, Indicates the first The average value of each peak (the location of the peak). Indicates the first The bandwidth (standard deviation) of each peak, since the time side follows a one-dimensional Gaussian distribution, i.e., the mean... It is a scalar, bandwidth It is also a scalar.

[0071] 2) On the spatial side, a method of superimposing principal components and global noise is adopted. The principal components include a static distribution and a dynamic distribution: a central static Gaussian distribution. and a dynamic distribution that uses a superposition of three Gaussian mixtures and varies with time. The noise follows a global standard distribution. . The weights of the principal components, The weights are for the noise, and must satisfy the following conditions: More specifically: the central static Gaussian distribution is located at the exact center of the geographical area of ​​the study region. Indicates its spatial location, The bandwidth (standard deviation) is used to simulate a static, large-scale trend where the center is dense and the edges are sparse, and it does not change over time. Three time-shifting Gaussian mixture distributions are used to simulate the positional changes of hotspots caused by the continuous movement and migration of events over time or due to the dynamic influence of certain factors, reflecting a dynamic process of change. Control the first The spatial location changes of the distribution Control the first The dynamic changes in the bandwidth (standard deviation) of each distribution, with these three peaks determined by weights. Control, and have . The weights are those of a centrally statically distributed system. The weights are dynamically distributed and satisfy the following conditions: Since this is a two-dimensional Gaussian function in space, the average value... All are vectors, bandwidth (standard deviation) All are matrices.

[0072] By superimposing the temporal and spatial distributions above, the resulting PDF can effectively simulate the complexity and variability of event distribution in an urban environment. By assigning different parameters to this generalized PDF, we obtain PDF1, PDF2, and PDF3, ensuring that these three PDFs become increasingly uneven. Specifically, the three PDFs use progressively smaller bandwidths to emphasize this uneven distribution characteristic; and a "single-variable method" is used for constraint, that is, only the bandwidth of the dynamic three peaks of the three PDFs is controlled to decrease, while other parameters remain consistent, avoiding interference caused by other unnecessary changes. The specific values ​​of each variable are shown in Tables 2 and 3.

[0073] Table 2. Parameter configurations for pdf1, pdf2, and pdf3 on the time side.

[0074]

[0075] Note: Based on the "single variable" principle mentioned in the text, the time distributions of pdf1, pdf2, and pdf3 are assumed to be consistent, i.e., they share a common set of parameters. The values ​​of μ and σ in the table represent the upper and lower bounds of the time interval (…). , The proportion of ).

[0076] Table 3. Parameter configurations of pdf1, pdf2, and pdf3 on the spatial side.

[0077]

[0078] Note: The increasingly uneven distribution of pdf1, pdf2, and pdf3 is shown in the table. The decision is made based on the "single variable" principle described in the text, while all other parameters remain unchanged. (Average vector) This table uses polar coordinates and has the following structure: ,in, This is a proportional value corresponding to half the total length in the X or Y direction. Let be the angle of rotation. Then, in rectangular coordinate form: , where is the planar location of the mean of a two-dimensional Gaussian distribution. Covariance matrix. use The way of writing, The value represents the proportion of 1 / 2 of the total length in the corresponding X or Y direction. It is a 2×2 identity matrix. (See table) , , They are located at 33% and 54% of the time axis, respectively.

[0079] Monte Carlo simulations were performed sequentially on pdf1, pdf2, and pdf3, generating 500 points each time. The same 500 points were then fitted using both TSTKDE and ASTKDE models from this application, yielding their respective joint probability density models. Then, starting from the top 1% of the regions with the highest density in both the TSTKDE and ASTKDE models, the ISE (integral mean square error) with the real background pdf was calculated, and then iterated through 2%, 3%, and so on up to 100% of the region. This not only validated the upper tail region (high-density hotspot) but also the global region (100%). This process was repeated 200 times for each of pdf1, pdf2, and pdf3, and the results were statistically analyzed to obtain the following... Figure 3 The experimental results shown are illustrated below. Blue represents the traditional Spatiotemporal Kernel Density Estimation Method (TSTKDE), and red represents the Spatiotemporal Kernel Density Estimation Method (ASTKDE) based on adaptive bandwidth proposed in this application. The plot is shown as the mean curve ± 95% confidence interval (mean ± 95% CI). As can be seen:

[0080] 1) ASTKDE has a smaller ISE across the entire PDF, both globally and at the head (high-value upper-lower region). This indicates that the algorithm proposed in this application outperforms traditional algorithms in all cases and can produce more reliable fitting results.

[0081] 2) From pdf1 to pdf3 (horizontal comparison), the data distribution becomes increasingly uneven, and the gap between the two models becomes larger and larger. Regardless of which optimal bandwidth selection method is chosen (the experiment used interpolation (PI), log-likelihood (LIK), and smooth leave-one-out cross-validation (SCV) as the leading optimal bandwidth), ASTKDE can produce more reliable fitting results than TASTKDE. The more uneven the data, the more significant this advantage becomes.

[0082] 3) With the PDF fixed, under different leading bandwidth conditions (vertical comparison), it can be seen that no matter which method is chosen to select the optimal bandwidth as the leading bandwidth, ASTKDE outperforms TSTKDE. This shows that the proposed ASTKDE can be based on multiple optimal bandwidth assumptions and can produce more reliable fitting than TSTKDE, demonstrating its universality and superiority, rather than being applicable only to a certain specific bandwidth.

[0083] Experiment 2: Comparative Evaluation of the Ability to Absorb the Impact of Extreme Events. Based on real data of road collapse disasters in Zhengzhou City from 2014 to 2023, this experiment analyzes the performance of ASTKDE and TSTKDE in the face of uncertainty under the impact of extreme events. It should be noted that Zhengzhou, Henan Province, experienced the "7.20" catastrophic rainstorm disaster in 2021, providing an ideal basis for this comparison of uncertainty absorption performance in "absorbing the impact of extreme events".

[0084] Uncertainty Design Experiment: 1. Data uncertainty: A yearly n-out-of-n (sampling with replacement) bootstrap method is used to simulate data-level uncertainty. This method simulates omissions and duplications during the data collection process and preserves the annual and monthly distribution characteristics of the data. 2. State and conceptual uncertainties are designed: a steady state (2014-2020) before the rainstorm (the 2021 Zhengzhou, Henan "7•20" catastrophic rainstorm disaster) and an unsteady state (2014-2021) during the rainstorm impact are defined. Models are built using TSTKDE and ASTKDE in the steady state (2014-2020) and the unsteady state (2014-2021) before the rainstorm, respectively. The 1) migration / transportation location changes and 2) quantity / volume size changes of the disaster event hotspots generated by the two models under these two states are evaluated to quantitatively assess the uncertainty exhibited by ASTKDE and TSTKDE under extreme event impact. The quantitative evaluation method is as follows: Hot zones are defined using two different methods: 1. Hot zones with fixed probability density cumulative mass (accumulating from the regions with the highest probability density until the total mass reaches 5% probability density, abbreviated as HDR-5%), 2. Hot zones with a fixed threshold (calculating the TSTKDE and ASTKDE models respectively on the 2014-2020 data without sampling of the full dataset, and taking the 95th percentile as the fixed threshold). and Any area exceeding its respective threshold is considered a hot zone. HDR-5% focuses on the highest density hot zones at the head by limiting quality. These hot zones are few in number but carry a high probability mass of monomers and are typical examples. It observes their positional migration / transport changes under extreme event impacts. Fixed threshold τ hot zones are numerous and have greater fluctuations. It focuses more on group behavior and observes the increase / decrease in the number / volume of hot zones above a certain fixed value under extreme event impacts.

[0085] The experiment was repeated 200 times. For each calculation result of ASTKDE and TSTKDE, the stability cubes of the HDR-5% hot zone and the fixed threshold τ were recorded. Specifically, for each result, hot zones were set to 1 and non-hot zones to 0, forming a (0,1) mask result. Then, the 200 0-1 mask results were summed together and divided by the number of repetitions (200) to approximate the probability of each voxel becoming a hot zone in these 200 iterations. In this way, the HDR-5% hot zone stability cubes and the fixed threshold τ hot zone stability cubes of TSTKDE and ASTKDE were formed.

[0086] Using the TSTKDE and ASTKDE algorithms, a model was built using steady-state sampling points from 2014 to 2020. The probability of events occurring after sampling within the impact window (2021) was then predicted, and the prediction score was calculated using the log-likelihood method. Detailed statistical results are shown in Table 4.

[0087] Table 4. Statistical Table of Uncertainty Results for TSTKDE and ASTKDE under Extreme Event Impacts

[0088] Note: In the table, "Hot Zone Expansion Ratio (Before / After Heavy Rain)" refers to the number / volume of hot zones in the unsteady state (2014-2021) divided by the number / volume of hot zones in the steady state before the heavy rain (2014-2020). "Hot Zone Ratio in the 2021 Impact Window" refers to the ratio of the number / volume of hot zones falling within the 2021 impact window to the total number / volume of all hot zones in the current model. Note that the "Hot Zone Retention" for HDR-5% is 0, and the "Hot Zone Ratio in the 2021 Impact Window" is 100%, so the corresponding p-value is meaningless and will not be given.

[0089] Table 4 clearly and intuitively shows that:

[0090] 1) For the HDR-5% hot zone, after the rainstorm impact, whether it is ASTKDE or TSTKDE, the hot zone with the first 5% mass all moved to the 2021 impact window. This well demonstrates the strong impact of the "7.20" extreme rainstorm event in 2021 on the entire system.

[0091] 2) It can be seen that at a fixed threshold Under these conditions, ASTKDE's heat region retention ratio and crossover ratio are higher than TSTKDE's, and the proportion of heat regions in the 2021 impact window is also higher than TSTKDE's. This indicates that ASTKDE exhibits "two-sided optimization": it can generate more aggregated and concentrated heat regions within the 2021 impact window, or allocate more "attention" to the impact window, while also retaining more "memory" of heat region skeletons from past historical windows (2014-2020), which TSTKDE cannot achieve. In short, it can retain more certainty amidst considerable uncertainty. This proves that ASTKDE's ability to absorb "extreme event shocks" is stronger than TSTKDE's. Under extreme event shocks, ASTKDE's model is sensitive, focused, and volatile. In contrast, TSTKDE is more sluggish and mediocre, unable to lock onto smaller, more concentrated heat regions.

[0092] 3) ASTKDE's prediction for 2021 (the impact year window) is about 8% better than TSTKDE. This shows that ASTKDE's "two-end optimization" is effective, improving the prediction accuracy for the impact window, and does not indicate overfitting.

[0093] like Figure 4 As shown in the visualization, TSTKDE and ASTKDE are displayed at a fixed threshold. Below are four stability cube plots under steady-state (2014-2020) and non-steady-state (2014-2021) conditions. It can be seen that ASTKDE can achieve a fixed threshold in both steady-state and non-steady-state conditions following extreme rainstorm impact. The hot zones are concentrated in a smaller area, while the hot zones generated by TSTKDE appear scattered and blurry.

[0094] Experiment 3: Comparison and evaluation of predictive ability and overfitting. A rolling window was used to establish 13 prediction experiments to further evaluate and demonstrate the predictive performance of ASTKDE over TSTKDE for events that will occur in the near future, and to verify whether overfitting occurred.

[0095] After absorbing the impact of the extreme rainstorm event on July 20, 2021, the model gradually returned to a new steady state (2022, 2023). We used the past 9 years as the training set to train the ASTKDE model and the TSTKDE model respectively, and used the unseen data for the next 4 months as the test set to verify the predictive performance of the two models for these 4 months (unseen future data). For the previous 9 years, we still used the same annual n-out-of-n (sampling with replacement) bootstrap method. Then, we rolled the prediction window, from January to April 2023 to January to April 2024, making a total of 13 predictions, with 100 samples for each prediction. The reason for choosing 4 months is that the field engineering construction cycle is approximately 4 months.

[0096] like Figure 5 As shown, the training and test sets for the above experimental methods were selected. Figure 6 As shown, ASTKDE produces more reliable and accurate estimates than TSTKDE in every prediction task. Its superior predictive performance on data not yet available in the near future also demonstrates that ASTKDE does not exhibit overfitting. By assigning a smaller bandwidth to the high-value upper tail region (hotspot), the hotspot also reveals more detail, facilitating downstream analysis tasks. This is an advantage that TSTKDE cannot achieve.

[0097] Experiment 4: Comparative evaluation of hot zone identification capabilities. Using the real data dataset of road collapse disasters in Zhengzhou from 2014 to 2023, the hot zone identification capabilities of TSTKDE and ASTKDE were compared under the same conditions.

[0098] We give the null hypothesis. Road collapse disasters in Zhengzhou from 2014 to 2023 within the effective temporal and spatial domain It follows a standard distribution (i.e., a distribution with no distributional characteristics). Then, in the effective spatiotemporal domain... Monte Carlo simulations were performed, with each simulation randomly generating the same number (1257) of sample points as the real dataset. Joint probability density models of TSTKDE and ASTKDE were then established. The statistic for each voxel becoming a hot zone is defined by the following formula:

[0099] ;

[0100] in, The joint probability density model of TSTKDE or ASTKDE derived from real datasets. Two probability density models, TSTKDE or ASTKDE, are generated using random points. The random number of times (the number of times the random process is repeated) is 999 times in this case. We select the significance level. To refuse This study argues that these voxels are the true hotspots identified by the TSTKDE or ASTKDE algorithms, breaking the null hypothesis of "random distribution" and revealing them as statistically significant hotspots with high frequency. This null hypothesis-based Monte Carlo verification method effectively eliminates noise and filters out more statistically significant hotspots. After 999 random iterations, the smaller the identified hotspot range / volume, the better the model's hotspot identification ability.

[0101] The results are as follows Figure 7 As shown, during 999 random trials, the hot zones identified by TSTKDE and ASTKDE gradually converged. ASTKDE had a smaller limit (5.05%), filtering out a smaller range of hot zones compared to TSTKDE (7.71%). This demonstrates the superior performance of the proposed ASTKDE algorithm in hot zone identification, locking onto more concentrated hot zones, while TSTKDE appeared more diffuse and ambiguous. Furthermore, it should be noted that ASTKDE assigns a smaller bandwidth to the hot zones, meaning the hot zones contain more detail, which is highly beneficial for downstream hot zone analysis tasks—capabilities lacking in TSTKDE.

[0102] Experiment 5: Sensitivity Analysis of ASTKDE. Experiments 1-4 compared and analyzed the advantages of the ASTKDE model in many aspects. However, the sensitivity behavior of the ASTKDE model itself needs to be further understood and analyzed to see how the model's behavior changes and whether it is robust when the key input parameters are changed or disturbed.

[0103] The experimental method is as follows: introduce bandwidth perturbation parameters. Then multiply by the base bandwidth respectively. , got 5 A 5-bandwidth perturbation grid, i.e., 25 bandwidth combinations, plus an edge correction factor. The on / off state generates 50 parameter combinations. Using pdf2 from Experiment 1 as the real background distribution, Monte Carlo simulations are performed, generating 500 points each time. The joint probability density model is then calculated for each of the 50 combinations. This process is repeated 200 times, with a baseline bandwidth of [missing information]. This represents the average bandwidth selected using plug-in interpolation over 200 iterations. Then, the ISE (Intersection over Union) of each combination within the edge region of the effective spatiotemporal cube is calculated, serving as the performance metric for the edge region. The area under the curve (AUC) of the intersection-over-union (IoU) curve between the identified and actual hot zones within the effective spatiotemporal cube for each combination is calculated, serving as the global performance metric. The edge region is defined as follows: the study area is simplified to a regular rectangle, and the three axes of the effective spatiotemporal cube... In the middle, respectively separated The most recent 20 voxels, and distance The edge region formed by the 10 most recent voxels, calculated to account for 29.7% of the total voxels and 5.47% of the total probability density, indicates a dense center and sparse edges. Since the true background distribution pdf2 is known, the true hot zones are also known. Here, hot zones are defined as the top 1%, 2%,...10% of their respective distributions, representing 10 distinct hot zones. Taking 1% as an example, the crossover ratio (CRO) is calculated between the true 1% hot zone and the respective 1% hot zones of the 50 ASTKDE-fitted models. Then, the CROs for 2%, 3%,...10% are calculated, resulting in 50 CRO curves. For each CRO curve, the area under the curve (AUC) is calculated as the indicator for that group, avoiding the bias caused by specifying only a single value (e.g., 5%) of the hot zone. The AUC (Area Under the Curve) method is a commonly used statistical method in machine learning.

[0104] like Figure 8 The results of the sensitivity analysis are shown in three subplots: left (a), middle (b), and right (c). The mean surface ± 95% confidence band surface of the response index is presented for all subplots. The mean surface is shown in... - The plane was projected and subjected to bilinear interpolation and gradient color display, visually illustrating the gradient changes. Several contour lines corresponding to the mean surface were also projected onto... - Plane. The magnitude of the gradient is also plotted for each grid point, with larger magnitudes indicated by red arrows.

[0105] for Figure 8The left figure (a) demonstrates that, in most cases, adding edge correction results in better performance at edges than not adding it, and the benefit of edge correction is bandwidth-dependent. The average plane is cut by the solid line at the bottom (where the pairing difference between q on and q off equals 0). - The projection on the plane serves as a watershed. Most grid points with a pairwise difference less than 0 are on the right side, indicating that in most cases, q on is better than q off, and in a few cases, q off is better than q on. This is actually the result of finite sampling under the influence of the classic "bias-variance trade-off" problem in parameterless estimation, which is perfectly acceptable.

[0106] for Figure 8 Figure (b) in the middle demonstrates that incorporating edge correction into the global performance results in a significantly better overall result than not. Since it's an AUC statistic, a higher value is better; qon minus qoff should be greater than 0 for optimal performance. (See Figure 5 in the middle.) In a 5-bandwidth grid, q on is better.

[0107] for Figure 8 Figure (c) on the right only shows the global ISE response surface of ASTKDE under 25 bandwidth perturbations of q on. It can be clearly seen that the global ISE average response surface is very smooth with bandwidth jitter / fluctuation, which proves that our proposed ASTKDE model is robust. That is, the model output is stable under bandwidth perturbation / fluctuation and there is no unstable situation of "the output jumps sharply when the bandwidth is slightly jittered".

[0108] Furthermore, we analyzed and demonstrated the space-side bandwidth. and time-side bandwidth To further reveal the behavior of the proposed innovative model, we examine the more sensitive details of which approach is more effective. By calculating the partial derivatives of the second-order (or first-order) difference for the 5×5 grid and the second-order difference for the inner 3×3 grid, as shown in Table 5, we can see that the proposed ASTKDE is more sensitive to spatial bandwidth, i.e., it exhibits a larger gradient in spatial bandwidth. However, it is only applicable to the baseline bandwidth. Within the surrounding small neighborhood, through observation Figure 8 As can be seen in the right figure (c), in Near the large bandwidth combination, the gradient of the temporal bandwidth actually exceeds that of the spatial bandwidth, resulting in a counterexample.

[0109] Table 5. Statistical table of spatial and temporal bandwidth gradients based on ASTKDE using finite difference.

[0110]

[0111] Note: In the gradient calculation process for the 5×5 bandwidth grid, second-order difference is used for the middle 3×3 grid, but first-order (forward or backward) difference is used for the edge grid to obtain approximate values. For greater accuracy, we additionally performed second-order difference gradient statistics for the middle 3×3 bandwidth grid. and These are the average values ​​of the absolute values ​​of the spatial and temporal bandwidth gradients, respectively. Here, the reference bandwidth... (rice), (Corresponding to 181.1 days), as mentioned above, is the average bandwidth obtained by using interpolation (plug-in) for all random processes.

[0112] Example 1:

[0113] By applying the proposed ASTKDE algorithm to urban road collapse disaster cases in Zhengzhou City from 2013 to 2023 (10 years), a joint probability model and a spatial conditional probability density model were generated. Qualitative and quantitative analyses of these two models respectively provided new insights into the spatial distribution and temporal evolution characteristics of high-incidence areas (hotspots) of road collapse disasters in Zhengzhou City, and predicted high-incidence areas in the short term. Furthermore, by overlaying the established models with the distribution of historical built-up areas and special soil types in Zhengzhou City from 1990 to 2020, new cutting-edge analytical results were obtained, demonstrating that the proposed algorithm model can also provide crucial basis for downstream analysis tasks, playing an important role and showing broad application prospects.

[0114] Taking 1257 damage records from 2014 to 2023 recorded by the Zhengzhou Municipal Urban Management Bureau as an example, the data was organized into three columns: X, Y, and T. The code was run to calculate a joint probability density model. The generated joint probability model was then subjected to the following steps: 1) generating local probability density peak points using the code; 2) volume rendering according to the shading rules, resulting in... Figure 9 3) Run the code, slice and project the data according to the year, calculate the area of ​​the hot zone with different statistical significance under the three levels, and obtain the results. Figure 10 As you can see: Figure 9 Figure 10Together, these findings depict the spatial distribution patterns and temporal evolution of road subsidence hotspots in Zhengzhou over the past 10 years. Spatially, the old city area serves as a significant boundary, while temporally, the extreme rainstorm disaster of July 20, 2021, acts as a crucial dividing line. The disaster ecosystem exhibits a three-stage change pattern: small at both ends and large in the middle. The entire disaster ecosystem has undergone a development trend of "development—extreme rainstorm impact—more severe." Before the extreme rainstorm, the hotspots were mainly concentrated within the old city area. Under the impact of the extreme rainstorm, the hotspots expanded dramatically, reaching a 10-year peak and extending beyond the constraints of the old city area. After the rainstorm, the area of ​​the hotspots did not return to pre-disaster levels but instead became even larger. This indicates that the road subsidence disaster situation in Zhengzhou has become more severe. The black dots represent local probability density peaks, which can be seen to be related to the annual flood season, meaning that rainfall significantly affects the temporal distribution of the hotspots.

[0115] like Figure 11 As shown, the spatial distribution model of road collapse disasters in Zhengzhou from 2014 to 2023 reveals how the spatial development and changes of the disaster have occurred: some areas were historically hot zones, then de-hot zones, and then became hot zones again in 2022 and 2023, reflecting the recurring nature of the disaster (activation and deactivation); torrential rains have profoundly altered the disaster ecosystem. The largest historical disaster spatiotemporal flow gradually disappeared after the torrential rains, replaced by new spatiotemporal flows that are generated and growing stronger. Therefore, in practical applications, relevant departments can use the joint probability distribution model obtained by the three-dimensional spatiotemporal kernel density estimation method to focus on deploying disaster prevention and control work in the new disaster spatiotemporal flow areas.

[0116] like Figure 12 As shown, by accumulating slices using a spatial conditional probability density model, the "temporal edge density distribution" in the figure was obtained, and contour lines were drawn. Overlaying this with images of Zhengzhou's historical built-up areas led to the conclusion that "areas built early on are often also areas prone to road collapse disasters." This reflects that the aging and leakage problems caused by the excessive service life of urban underground infrastructure (especially underground pipe networks) are a significant factor contributing to road collapse disasters in Zhengzhou. Furthermore, from several key areas in the figure (Area 1, Area 2, Area 3), it can be inferred that the maximum service life of underground pipelines in Zhengzhou is 30 years (this is a rough estimate, ignoring factors such as the type of underground pipelines and construction standards, but it is of great significance for urban macro-governance).

[0117] like Figure 13As shown, by accumulating slices using a spatial conditional probability density model, the "temporal marginal density distribution" in the figure was obtained, and contour lines were drawn. Analysis was then performed by overlaying this onto unfavorable soil types in Zhengzhou (including slightly and moderately collapsible loess, liquefiable soil, and soft soil). The results showed that the widely distributed collapsible loess in Zhengzhou is not the main cause of road collapses. This refutes the view held by some experts and scholars that the collapsible loess in Zhengzhou is the primary cause of frequent and severe road collapses.

[0118] pass Figure 12 , Figure 13 The spatial conditional probability density model yields a profound and cutting-edge new understanding of road collapse disasters in Zhengzhou, demonstrating the important role of the innovative algorithm model disclosed in this application for downstream analysis tasks or interdisciplinary applications in other fields.

[0119] The adaptive bandwidth-based three-dimensional spatiotemporal kernel density estimation method proposed in this application can not only analyze geospatial events (XY being the geospatial coordinates of the event, and T being time), but also, with the support for "XY anisotropy" calculation in the code, can be extended to other disciplines with three-dimensional events (three-dimensional points). For example, in physics, it can be used for three-dimensional statistics such as velocity, temperature, and mass of each object (particle); in chemistry, it can be used for three-dimensional statistics such as gold content, copper content, and iron content of each sample; in school statistics, it can be used for three-dimensional quantities such as height, weight, and age of each student sample; in computer graphics, it can be used to statistically analyze the joint distribution of R, G, and B values ​​of each pixel in an image, and so on.

[0120] Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

Claims

1. A three-dimensional spatiotemporal kernel density estimation method based on adaptive bandwidth, characterized in that, Includes the following steps: S1: Determine the geographical boundary of the study area and the lower and upper time bounds of the target events. Based on the boundary, upper, and lower time bounds, obtain the effective spatiotemporal domain. Collect target event samples within this domain to obtain the target event dataset. Fit the true background distribution of the target events to the target event dataset, thus obtaining the joint probability density model. The specific joint probability density function is: ； ； ； ； ； ； ； in, The joint probability density function is used; the boundary contour is a closed polygon on the XY plane of geographic space. The upper bound of time is The upper bound of time is , composed of closed polygons and upper bound of time, lower bound of time The enclosed area is the effective spatiotemporal domain. , composed of closed polygons Minimum bounding rectangle and upper and lower time bounds The enclosed cube is a spacetime cube. In the effective spatiotemporal domain Internal collection Event Samples Each event sample That is, each event sample contains geospatial coordinates. , and time coordinates , For the first The horizontal coordinate of each event sample in geographic space; For the first The vertical coordinate of each event sample in geographic space; For the first The coordinates of an event sample in time and space; It is an indicator function, for each event sample. The contributed kernel density falls into the effective spatiotemporal domain Take 1 if the condition is met, otherwise take 0; , For each event sample The kernel function based on adaptive bandwidth is specifically the Gaussian kernel function. In This can be either in the general formula or It could also be This uses the assumption of XY isotropy, that is, any event After it occurs, the impact is consistent in all directions across space. For the general formula ; For the first Spatial adaptive bandwidth for each event sample For the first Time-side adaptive bandwidth for each event sample; for The optimal fixed bandwidth on the spatial side for each event sample; for Optimal fixed bandwidth on the time side for each event sample; For global space-side bandwidth; For global time-side bandwidth, let's directly set... , ; for The leading spatiotemporal kernel density estimation model corresponding to each target event sample The geometric mean of the kernel density values. It uses a fixed optimal bandwidth and The leading spatiotemporal kernel density estimation model at sample points The kernel density value at that location; This is the edge correction factor; S2: Time-by-time slicing is performed on the joint probability density model to obtain the spatial conditional probability density model. The specific spatial conditional probability density function is as follows: ； in, It is the spatial conditional probability density function; S3: Perform qualitative and quantitative analysis on the joint probability density model and the spatial conditional probability density model according to the requirements.

2. The three-dimensional spatiotemporal kernel density estimation method based on adaptive bandwidth according to claim 1, characterized in that: The voxel differentiation rules of the joint probability density model include: ASTKDE-Large model and ASTKDE-Base model.

3. The three-dimensional spatiotemporal kernel density estimation method based on adaptive bandwidth according to claim 2, characterized in that: The differential voxel parameters corresponding to the ASTKDE-Large model are: X = 695, Y = 735, T = 365, and the corresponding voxel spatiotemporal scale is 50 m × 50 m × 10 days; the differential voxel parameters corresponding to the ASTKDE-Base model are: X = 348, Y = 368, T = 183, and the corresponding voxel spatiotemporal scale is 100 m × 100 m × 20 days.

4. The three-dimensional spatiotemporal kernel density estimation method based on adaptive bandwidth according to claim 1, characterized in that: The joint probability density model employs concurrent computing and load balancing during the calculation process.

5. The three-dimensional spatiotemporal kernel density estimation method based on adaptive bandwidth according to claim 1, characterized in that, The specific process of performing qualitative and quantitative analysis on the joint probability density model in step S3 is as follows: obtaining statistical significance. The local density peaks of the region are obtained by importing the joint probability density model into the volume rendering software. The volume rendering gradient shading is performed based on the kernel density value of each voxel. The obtained local density peaks are also imported into the scene. The probability density hotspots in the model are analyzed by displaying the volume rendering and these peaks together.

6. The three-dimensional spatiotemporal kernel density estimation method based on adaptive bandwidth according to claim 5, characterized in that, The volume rendering gradient coloring rule is: Red (salience) → Yellow (significance) → Blue (significance) → Completely transparent, with a gradient transition coloring based on the salience thresholds corresponding to red, yellow, and blue.

7. The three-dimensional spatiotemporal kernel density estimation method based on adaptive bandwidth according to claim 1, characterized in that, The specific process of qualitative and quantitative analysis of the joint probability density model in step S3 is as follows: In RStudio software, the joint probability density model is plotted as a three-dimensional isosurface and cut by year. The model fragments of each year are projected onto the XY plane to form the maximum projection range. The local density peak points in the year are also projected onto the XY plane, and the area of ​​the maximum projection range is calculated. The specific results of the annual slices and the corresponding statistical situation of the hot zone area with different significance are analyzed.

8. The three-dimensional spatiotemporal kernel density estimation method based on adaptive bandwidth according to claim 7, characterized in that, The isosurface coloring rule is: red (salience). ), yellow (significant) ), blue (significance) Coloring is done according to the significance thresholds corresponding to red, yellow, and blue.

9. The three-dimensional spatiotemporal kernel density estimation method based on adaptive bandwidth according to claim 1, characterized in that, The specific process of performing qualitative and quantitative analysis of the spatial conditional probability density in step S3 is as follows: The statistical significance of each time slice is calculated using the spatial conditional probability density model. Local density peak points within the range are connected across adjacent time slices according to a set minimum distance threshold to form the migration trajectory of the hot zone over time, and the migration trajectory of the hot zone over time is analyzed.

10. The three-dimensional spatiotemporal kernel density estimation method based on adaptive bandwidth according to claim 1, characterized in that: The optimal fixed bandwidth on the space side Optimal fixed bandwidth on the time side It can be obtained through leave-one-out cross-validation, interpolation, log-likelihood, or empirical rules.

Images