A rain anomaly identification and repair method based on bidirectional propagation

Abstract

The application discloses a rain anomaly identification and repair method based on bidirectional propagation, which comprises the following steps: establishing a rainfall propagation sequence, selecting bidirectional witness stations, forming bidirectional propagation constraints, and performing consistent gate identification and repair; according to the arrival time and spatial position of a rainfall event, a propagation arrival time plane is fitted, the propagation direction of the rainfall is determined, each rainfall station is projected to the direction to establish a propagation arrival sequence, and then upstream witness stations and downstream witness stations are selected according to the distance along the propagation direction and the lateral offset distance; the method makes the selection of the witness stations match the moving path of the target rainfall process, and solves the problem that the fixed geographical neighborhood in the prior art cannot reflect the actual propagation stage of a rain cluster, thereby causing systematic deviation in spatial consistency judgment.

Description

Technical Field

[0001] This invention relates to the field of hydrological monitoring data processing technology, specifically to a method for identifying and repairing rainfall anomalies based on two-way propagation. Background Technology

[0002] Rain gauges are crucial monitoring facilities for obtaining information on the spatiotemporal distribution of regional rainfall. Their observational data are widely used in flood forecasting, water resource allocation, flash flood warnings, and hydrological analysis. Due to factors such as sensor malfunctions, communication interruptions, equipment aging, power supply anomalies, and external environmental interference, rainfall data may exhibit issues such as sudden increases, abnormally low values, or consecutive missing measurements. If abnormal data is directly entered into subsequent operational systems without identification, it may affect the calculation of areal rainfall and the results of hydrological forecasts.

[0003] Existing rainfall anomaly identification and repair technologies mainly fall into two categories: one is based on the time series characteristics of the target rain gauge station itself, using fixed thresholds, sliding statistics, or differences between adjacent time points to determine anomalies, and then using historical averages or trend extrapolation for repair; the other is to introduce spatial neighbor station information, using geographically adjacent stations around the target rain gauge station to determine spatial consistency, and using inverse distance interpolation, spatial statistical interpolation, or neighbor station averages to generate repair estimates.

[0004] In practical applications, the aforementioned existing technologies still have the following problems: First, rainfall processes are spatiotemporally non-stationary, and fixed geographical neighborhoods cannot reflect the actual propagation sequence of rain clusters. The selection of geographically neighboring stations is based solely on fixed spatial distances, without considering the actual propagation direction of the rainfall process and the time difference of arrival between stations. This results in the selected neighboring stations potentially being at different stages of rainfall propagation, causing systematic biases in spatial consistency judgments and repair estimates. Second, repair estimates typically only utilize observational evidence from a single direction or fixed neighborhood of the target rain gauge station. When neighboring stations in different directions provide contradictory estimates for the target station, automatic repair may still be triggered, making it impossible to distinguish between real local rainfall enhancement and monitoring anomalies, resulting in the loss of effective extreme value information. Summary of the Invention

[0005] The present invention proposes a method for identifying and repairing rainfall anomalies based on bidirectional propagation, which can at least solve one of the technical problems in the background art.

[0006] To achieve the above objectives, the present invention adopts the following technical solution: A method for identifying and repairing rainfall anomalies based on bidirectional propagation, comprising: Based on the rainfall data of each rain gauge station within a rainfall event, the arrival time of the event at each rain gauge station is determined. Plane fitting is performed based on the spatial location of the rain gauge station and the arrival time of the event to determine the direction of rainfall propagation. Based on the projected position of the rain gauge station in the direction of rainfall propagation, the propagation arrival sequence of the rain gauge stations is established. Based on the propagation arrival order of rain gauges, the candidate rain gauges are constrained by the propagation direction and the lateral offset distance to obtain the upstream witness station set and the downstream witness station set. The propagation arrival time difference is determined based on the arrival times of events at the witness station and the target rain gauge station. The phase correction is performed on the rainfall sequence of the witness station, and the forward rainfall estimate and the reverse rainfall estimate are generated by fusion respectively. The bidirectional propagation constraint interval of the target rain gauge station is constructed by combining the internal dispersion. Based on the bidirectional propagation constraint interval and forward and reverse rainfall estimates, the anomaly types of the target rain gauge are identified, distinguishing between missing anomalies, low-value anomalies, and high-value anomaly candidates. For missing anomalies, direct fusion and repair are performed. For low-value anomalies, repair or retention is performed after confirmation based on upstream and downstream rainfall increments. For high-value anomaly candidates, the a priori and posterior windows of the witness station are determined based on the propagation arrival time difference between the target rain gauge and the witness station. Within each window, the rainfall enhancement is calculated and intensity and time delay matching judgment is performed. The propagation posterior confirmation result is obtained based on the judgment result. Based on the propagation posterior confirmation result and the fusion value of forward and reverse rainfall estimates, the repaired rainfall value is generated.

[0007] Preferably, the process of determining the arrival time of the event at the rain gauge station includes: calculating the cumulative rainfall and the total rainfall of the rainfall event at each rain gauge station up to the current sampling time; determining the moment when the cumulative rainfall first reaches a preset proportion of the total rainfall of the rainfall event as the arrival time of the event; and when the total rainfall of the rainfall event is lower than a preset minimum rainfall threshold, the corresponding rain gauge station does not participate in the subsequent propagation direction calculation.

[0008] Preferably, the process of determining the direction of rainfall propagation in a rainfall event includes: fitting a propagation arrival time plane using the least squares method and calculating the fitting residual for the arrival time of each rain gauge station; identifying and excluding rain gauge stations that deviate from the propagation arrival time plane based on the absolute deviation of the median of the residuals, and then refitting; determining the direction of rainfall propagation based on the spatial gradient of the refitted plane, and arranging each rain gauge station in ascending order of its projection position in the propagation direction to obtain the order of propagation arrival.

[0009] Preferably, the process of obtaining the upstream witness station set and the downstream witness station set includes: calculating the distance along the propagation direction and the lateral offset distance between the candidate rain gauge station and the target rain gauge station; screening candidate rain gauge stations whose distance along the propagation direction and lateral offset distance both meet the preset distance thresholds, and using them as the upstream candidate witness station set and the downstream candidate witness station set, respectively; calculating the witness weight for the candidate witness stations according to the distance attenuation method; sorting the candidate stations according to the witness weights and selecting the candidate stations with a preset upper limit, and obtaining the upstream witness station set and the downstream witness station set, respectively; and not performing automatic repair if the number of stations in either set is less than the preset minimum number.

[0010] Preferably, the process of phase correction of the rainfall sequence of the witness station includes: calculating the propagation arrival time difference based on the event arrival times of the target rain gauge station and the witness station; shifting the rainfall sequence of the witness station to the time phase of the target rain gauge station according to the propagation arrival time difference; and correcting the shifted rainfall according to the proportion of the total rainfall of the two stations or the proportion of historical responses. The fusion process of generating forward rainfall estimates and reverse rainfall estimates does not use the observed rainfall at the target rain gauge station at the current time.

[0011] Preferably, the process of fusing and generating forward and reverse rainfall estimates includes: when the witness station is located upstream of the propagation, the rainfall sequence of the upstream witness station is forward-shifted by the propagation arrival time difference to obtain the upstream corrected rainfall estimate; when the witness station is located downstream of the propagation, the rainfall sequence of the downstream witness station is backward-shifted by the propagation arrival time difference to obtain the downstream corrected rainfall estimate; the upstream corrected rainfall estimate and the downstream corrected rainfall estimate are respectively fused by weighted median according to the witness weight to obtain the forward and reverse rainfall estimates, and the internal dispersion of each direction is determined according to the deviation between the internal corrected rainfall and the estimated value of the corresponding direction.

[0012] Preferably, the process of constructing the bidirectional propagation constraint interval of the target rain gauge station includes: calculating the bidirectional fusion center value and the bidirectional estimation difference between the forward and reverse rain gauge estimates; calculating the bidirectional merged discrete quantity by combining the upstream discrete quantity, the downstream discrete quantity, and the bidirectional estimation difference; calculating the half-width of the constraint interval based on the bidirectional merged discrete quantity, the absolute uncertainty parameter, and the relative uncertainty parameter; and obtaining the lower and upper boundaries of the bidirectional propagation constraint interval by combining the bidirectional fusion center value.

[0013] Preferably, the process of identifying the anomaly type of the target rain gauge includes: calculating the bidirectional consistency based on forward and reverse rainfall estimation, and determining the automatic repair threshold value in combination with the number of valid witness stations; when the automatic repair threshold value meets the preset conditions, determining whether the original observed rainfall of the target rain gauge exceeds the bidirectional propagation constraint interval; when the original observed rainfall is missing, it is determined to be a missing measurement anomaly; when it exceeds the upper boundary of the bidirectional propagation constraint interval, it is determined to be a high-value anomaly candidate; when it is below the lower boundary of the bidirectional propagation constraint interval, it is determined to be a low-value anomaly.

[0014] Preferably, the process of calculating the rainfall enhancement and performing intensity and time delay matching judgment includes: calculating the magnitude by which the observed value of the target station exceeds the upper boundary of the two-way propagation constraint interval to obtain the high-value excess; determining the theoretical prior center time based on the propagation arrival time difference between the target rainfall station and the upstream witness station, and using the maximum period rainfall within the upstream prior window as the upstream rainfall enhancement; determining the theoretical posterior center time based on the propagation arrival time difference between the target rainfall station and the downstream witness station, and using the maximum period rainfall within the downstream posterior window as the downstream rainfall enhancement; filtering witness stations whose lateral offset distance exceeds a preset offset threshold, calculating the proportion of the upstream rainfall enhancement to the high-value excess and the peak time delay deviation, calculating the proportion of the downstream rainfall enhancement to the high-value excess and the peak time delay deviation, and when either proportion or time delay deviation meets the corresponding matching condition, determining that the propagation posterior confirmation is passed and retaining the original observation; otherwise, determining that the propagation posterior confirmation is not passed and repairing with the fused value.

[0015] Preferably, the process of generating and repairing rainfall values ​​includes: for low-value anomalies, determining whether there is an increase in rainfall at the corresponding propagation time upstream or downstream; if there is an increase in rainfall upstream or downstream, the low value is confirmed as an anomaly and repaired using a two-way fusion central value; if there is no increase in rainfall upstream or downstream, the low value is confirmed as a real reduction in rainfall, and the original observed rainfall is retained; for high-value anomaly candidates, when the post-propagation confirmation result is negative, the two-way fusion central value is used for repair; when the post-propagation confirmation result is positive or the evidence is insufficient, the original observed rainfall is retained and a verification mark is output; for missing measurement anomalies, the two-way fusion central value is directly used for repair.

[0016] As can be seen from the above technical solution, the rainfall anomaly identification and repair method based on bidirectional propagation of the present invention aims to address the problem that existing rainfall anomaly processing methods rely on fixed geographical neighbor stations or one-way estimation, which easily misjudge real local enhanced rainfall as anomalies and implement replacements.

[0017] 1. This invention fits the propagation arrival time plane based on the cumulative rainfall arrival time and spatial location of each rain gauge station in a rainfall event. The direction of rainfall propagation is determined by the spatial gradient of the plane, and each rain gauge station is projected onto this direction to establish the propagation arrival sequence. This makes the selection of upstream and downstream witness stations match the movement path of the target rainfall process, solving the problem that fixed geographical neighborhoods cannot reflect the actual propagation stage of the rain cluster, leading to systematic deviations in spatial consistency judgment.

[0018] 2. This invention selects witness stations from both upstream and downstream of the propagation path. Based on the propagation time difference between stations, the rainfall sequence of the witness stations is phase-corrected to generate forward and reverse rainfall estimates. Neither estimate uses the observed rainfall at the target rain gauge station at the current moment as input. The difference between the estimates in the two directions and the degree of dispersion within each direction are used to construct a bidirectional propagation constraint interval. This ensures that the reference for anomaly determination comes from independent observational evidence on both sides of the target station, overcoming the problems of lack of reverse verification and inability to distinguish between real local rainfall enhancement and monitoring anomalies in unidirectional estimation.

[0019] 3. For high-value data exceeding the constraint range, this invention further introduces a post-hoc confirmation mechanism for the propagation process. Based on the propagation time difference between stations, the a priori window of the upstream witness station and the posterior window of the downstream witness station are determined. By analyzing the intensity ratio and time delay deviation between the rainfall enhancement amount within the window and the excess amount at the target station, it is determined whether the high value has temporal evolution evidence of real rainfall propagation. Only isolated anomalies that do not have propagation interpretability are repaired, while the original observations of real local enhanced rainfall that have propagation interpretability are retained. This completely solves the technical problem of existing technologies that easily misjudge real local heavy rainfall as high-value anomalies and incorrectly replace them. Attached Figure Description

[0020] Figure 1 This is a flowchart of the method of the present invention; Figure 2 A flowchart for constructing a bidirectional propagation constraint interval is provided for this invention. Figure 3 This is a diagram illustrating the propagation sequence of the upstream and downstream witness stations of this invention; Figure 4 This invention uses a combination of forward estimation and backward regression to approximate the target rainfall map. Figure 5 This is a schematic diagram of the antecedent and consequent windows for high-value anomalies in this invention.

[0021] Figure 6 This diagram illustrates the mechanism by which localized enhanced rainfall is preserved in accordance with the present invention. Figure 7 This is a graph illustrating the anomaly identification and retention indicators of the present invention; Figure 8 This is a diagram showing the number of automatic repair actions in this invention. Detailed Implementation

[0022] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are some embodiments of the present invention, but not all embodiments.

[0023] Example 1, as Figure 1As shown in the figure, the rainfall anomaly identification and repair method based on bidirectional propagation described in this embodiment includes steps S1 to S4.

[0024] This method uses a single rainfall event as the processing unit. First, it determines the direction of rainfall propagation and establishes the station propagation arrival order based on the arrival time and spatial location of the cumulative rainfall at each rain gauge station within the event. Second, for the target rain gauge station to be identified, upstream and downstream witness stations are selected. Subsequently, the rainfall sequence at the witness stations is phase-corrected based on the inter-station propagation time difference. The upstream witness station generates a forward rainfall estimate, and the downstream witness station generates a reverse rainfall estimate. A bidirectional propagation constraint interval is constructed based on the estimates in both directions and their dispersion. The method is further defined as follows: if the bidirectional estimates meet the consistency condition and the observed rainfall at the target rain gauge station exceeds the constraint interval or is missing, the propagation is considered complete. When an anomaly is identified, the target rain gauge data is designated as anomaly candidates. Among them, missing anomalies and low-value anomalies are directly processed through the repair branch. Low-value anomalies and high-value anomalies undergo posterior verification of the propagation process, and the decision on whether to perform automatic repair is based on the posterior verification result. For high-value anomalies, the intensity and time matching relationship between the rainfall enhancement in the upstream and downstream posterior windows and the high-value excess at the target station is used to determine whether the high value has real rainfall propagation interpretability. For low-value anomalies, the existence of corresponding rainfall enhancement in the upstream and downstream areas is used to determine whether the low value at the target station is an isolated anomaly. Finally, the decision on whether to perform automatic repair is based on the posterior verification result of the propagation process.

[0025] Before implementation, this method first divides rainfall events. It then uses the continuous dry interval method to separate independent rainfall events from a continuous rainfall monitoring sequence. This method is applied to each rain gauge station to obtain the start and end times of their respective rainfall events. The start and end times of a regional rainfall event are taken as the union of the rainfall event times of all valid rain gauge stations within that event, resulting in a unified regional rainfall event time boundary. A regional rainfall event is only considered valid if the total rainfall at at least one rain gauge station within that event is greater than or equal to V. min Only when the rainfall event in the area proceeds to the subsequent step S1 processing; when the rainfall at each sampling time interval within a continuous ΔT hour is 0, the previous rainfall event is determined to have ended; when a sampling time interval with a rainfall greater than 0 occurs, a new rainfall event is determined to have begun; the value of ΔT ranges from 1 to 3 hours, and the specific value is determined based on the regional rainfall characteristics and monitoring time resolution; at the same time, a minimum total rainfall threshold V for the event is set. min When the total rainfall at all rain gauges during a rainfall event is lower than V min At that time, it was determined that the process did not constitute an independent rainfall event and was not included in subsequent processing; V min The recommended value range is 0.5 mm to 2.0 mm.

[0026] Step S1: Establish the rainfall propagation sequence: Based on the arrival time and spatial location of multiple rain gauges within a rainfall event, determine the direction of rainfall propagation for this rainfall event, and establish the rain gauge propagation arrival sequence according to the projection position of each rain gauge in the direction of rainfall propagation.

[0027] Step S1.1: Determine the arrival time of the event.

[0028] The cumulative rainfall ratio method is used to process the rainfall sequence of each rain gauge station within a single rainfall event to determine the event arrival time for each rain gauge station.

[0029] Suppose a rainfall event includes I rain gauges and K sampling times. The rainfall at rain gauge i during the k-th sampling time period is denoted as... The cumulative rainfall at rain gauge i up to the k-th sampling time. Calculate according to the following formula:

[0030] Where i represents the rain gauge station number; k represents the current sampling sequence number; and m represents the cumulative summation sequence number; This represents the m-th sampling time. Indicates that rain gauge i is at Rainfall at any given moment; This represents the cumulative rainfall at rain gauge station i up to the k-th sampling time. The formula accumulates the non-negative rainfall at each sampling time to calculate the cumulative rainfall sequence of rain gauge station i at each sampling time, which is used to subsequently calculate the total rainfall of the event and the arrival time of the event.

[0031] Total rainfall of the event at rain gauge i Calculate according to the following formula:

[0032] Where K represents the total number of samples included in this rainfall event; This represents the total rainfall of rain gauge i in this rainfall event. The formula takes the cumulative rainfall after all sampling of the rainfall event is completed as the total rainfall of the event. The output result is used to determine the arrival time of the event and can also be used to determine the rainfall response correction coefficient in the subsequent step S3.1.

[0033] The moment when the cumulative rainfall first reaches a preset percentage of the total rainfall for the event is defined as the event arrival time. The sampling sequence number at which the cumulative rainfall first reaches the preset percentage is also defined. Determined according to the following formula:

[0034] in, Indicates the proportion of events that arrive, satisfying ; This indicates that the cumulative rainfall has reached [a certain level] for the first time. The sampling sequence number is used to calculate the earliest sampling sequence number that satisfies the cumulative rainfall threshold. The output sampling sequence number is used to map the event arrival time.

[0035] Event arrival time of rain gauge i Determined according to the following formula:

[0036] in, Indicates the first Each sampling sequence number corresponds to a sampling time; this formula maps the sampling sequence number that meets the conditions to the real time, and outputs the arrival time of events at each rain gauge station. The output is used to fit the propagation arrival time plane in step S1.2, and is also input into step S3.1 to calculate the inter-station propagation arrival time difference.

[0037] Using the cumulative rainfall ratio to determine the event arrival time can reduce the impact of sporadic pre-event rainfall and individual noise values ​​on the propagation sequence; in this embodiment, Pick , The value of is related to the rainfall pattern characteristics and sampling frequency of the rainfall event; if If the value is too small, such as less than 0.3, it is easy to misjudge sporadic, invalid early rainfall or sensor noise as the start of an event, leading to distortion in the propagation direction fitting; if If the value is too large, such as greater than 0.7, effective information from the initial stage of rainfall will be ignored, causing the event start point to lag behind the actual arrival time of the rain cluster. This is especially problematic for slowly varying rainfall patterns, severely reducing the accuracy of propagation sequence establishment. Those skilled in the art can select the value within the range of 0.3 to 0.7 that best matches the event start and end times with the radar echoes of the rain cluster movement, based on the cumulative rainfall curves of historical rainfall events. .

[0038] When the total rainfall V at rain gauge i i When the value is 0, the station does not calculate the arrival time of the event. It does not participate in the propagation direction calculation; when V i If the rainfall is greater than 0 but lower than the preset minimum rainfall threshold for a single station, the event arrival time of that station will not be used in the propagation direction calculation; the preset minimum rainfall threshold for a single station and the minimum total rainfall threshold V for the event. min Take the same value.

[0039] Step S1.2: Determine the direction of rainfall propagation.

[0040] Rain gauge The two planar spatial coordinates are denoted as follows: and A propagation arrival time plane is established based on the spatial location of multiple rain gauges and the arrival time of the event.

[0041] in, and These represent the first spatial coordinate and the second spatial coordinate of rain gauge i, respectively. The constant term representing the time of arrival of the propagation plane; and These represent the coefficients of change of the event arrival time along the two spatial coordinate directions, respectively; Indicates rain gauge station The arrival time is fitted with the residual; this formula obtains the plane model parameters through least squares fitting, and outputs the plane parameters. , , and the fitting residuals at each station , The planar spatial gradient is used for subsequent robust screening of abnormal sites and is used to calculate the direction of rainfall propagation.

[0042] During the fitting process, the fitting residuals of each rain gauge station were used as a basis. i is used for robust screening to reduce the impact of a few abnormal arrival times on the propagation direction; specifically, the absolute deviation of the residual median can be used to identify rain gauges that deviate from the propagation plane, and n is used as the basis for selection. ε The absolute deviation of the residual median is used as the deviation judgment threshold, n ε It is recommended to use 3.0, and refit after excluding off-site stations; the number of valid rain gauge stations participating in the final fitting should be no less than 4. If there are fewer than 4, the calculation of the propagation direction of this rainfall event should be abandoned and the subsequent automatic repair steps should not be entered.

[0043] The direction of rainfall propagation is determined based on the spatial gradient of the plane at the time of arrival, and the component of the propagation direction per unit rainfall is calculated:

[0044] in, and These represent the components of the unit rainfall propagation direction in the two spatial coordinate directions, respectively. The formula normalizes the spatial gradient to obtain the unit propagation direction of the rain cluster movement. This direction points to the region where the arrival time of the event gradually increases, that is, the direction in which the rain cluster moves from the region that arrives earlier to the region that arrives later.

[0045] when and When both are 0, the gradient magnitude is 0, directly determining that this rainfall event cannot form a stable propagation direction, and it will not proceed to the subsequent automatic repair steps; when... Less than the preset stability threshold η stabAt that time, it was determined that this rainfall event could not form a stable propagation direction, and therefore would not proceed to the subsequent automatic repair steps; η stab A speed of 0.1 h / km is recommended.

[0046] Step S1.3: Establish the arrival sequence of rain gauge propagation.

[0047] Project each rain gauge station onto the direction of rainfall propagation and its perpendicular direction:

[0048]

[0049] in, Indicates rain gauge station The projected position along the direction of rainfall propagation; Indicates rain gauge station The lateral projection position perpendicular to the direction of rainfall propagation.

[0050] according to Arrange the rain gauges in ascending order to obtain the propagation arrival order corresponding to this rainfall event; for any two rain gauges, when one rain gauge's... When it is small, it is considered as a relative upstream station of propagation; when it is small... When the value is large, it is considered as a downstream station of relative propagation.

[0051] The direction of rainfall propagation was obtained after processing. Event arrival times at each rain gauge station Projected position along the propagation direction Lateral projection position And the order of arrival of rain gauges, for use in step S2 to select bidirectional witness stations; such as Figure 3 As shown, the projection positions of each rain gauge station are obtained by projecting along the direction of rainfall propagation.

[0052] Step S2: Select bidirectional witness stations: Based on the rain gauge arrival order obtained in Step S1, select upstream and downstream witness stations for the target rain gauge to be identified. The selection of witness stations is constrained by both the inter-station distance along the rain propagation direction and the lateral offset distance perpendicular to the rain propagation direction.

[0053] Step S2.1: Divide upstream and downstream candidate witness stations.

[0054] The target rain gauge station to be identified is denoted as... Other rain gauges are numbered as follows: Based on the projected position along the propagation direction obtained in step S1 and lateral projection position Calculate candidate witness stations With the target rain gauge Distance along the propagation direction and lateral offset distance between them:

[0055]

[0056] in, and These represent the target rain gauge stations. The projection position along the propagation direction and the lateral projection position; Indicating candidate witness station With the target rain gauge The distance between them along the direction of propagation; This indicates the lateral offset distance between the two.

[0057] when At that time, candidate witness station Located at the target rain gauge station The upstream of the spread; when At that time, candidate witness station Located at the target rain gauge station The downstream spread.

[0058] The upstream candidate witness station set has been determined as follows:

[0059] The downstream candidate witness station set has been determined as follows:

[0060] in, Indicates the target rain gauge station The collection of upstream candidate witness stations; Indicates the set of downstream candidate witness stations; Indicates the maximum permissible distance along the propagation direction; This represents the maximum allowed lateral offset distance. The above formula applies distance constraints to filter the set of candidate witness stations. By setting two distance limits, it is possible to avoid selecting rain gauge stations that are geographically close but significantly deviate from the same rainfall propagation path.

[0061] when When the set is empty, the upstream candidate set is empty, and the final upstream witness station set is empty. Empty; when When the value is empty, the downstream candidate set is empty, and the final downstream witness station set is empty. An empty set in any direction will result in an insufficient number of witness stations in that direction.

[0062] Step S2.2: Calculate the candidate witness weights.

[0063] For candidate witness station i that meets the distance condition, the witness weight is calculated based on the distance along the propagation direction and the lateral offset distance:

[0064] in, Indicating candidate witness station Relative to the target rain gauge Witness weight; This indicates the distance attenuation scale along the propagation direction; This indicates the lateral distance attenuation scale.

[0065] When the distance along the propagation direction or the lateral offset distance between the candidate witness station and the target rain gauge station increases, its witness weight decreases accordingly. Therefore, the witness weight reflects not only the relationship before and after the rain cluster propagation, but also whether the candidate witness station and the target rain gauge station are located on similar propagation paths.

[0066] and The value of is related to the spatial correlation scale of typical rainfall systems in the region; preferably, it can be determined by statistically analyzing the decay curve of the cross-correlation coefficient of hourly rainfall sequences of adjacent rain gauges in historical heavy rainfall events within the study area as a function of distance, taking the cross-correlation coefficient decay to e. -1 The corresponding distances along the propagation direction and the lateral distance are respectively used as... and The value is determined based on the following: when detailed statistical data is unavailable, it can be set empirically based on 2 to 4 times the average spacing between regional rain gauge stations.

[0067] In this embodiment, based on the watershed range, Pick , Pick , Desirable , Desirable The above parameters can be adjusted according to the density of rain gauges and the regional rainfall scale.

[0068] Step S2.3: Determine the set of bidirectional witness stations.

[0069] Each upstream candidate witness station set and downstream candidate witness stations Candidate stations are ranked according to witness weight. Sort by size from largest to smallest, and select no more than [number] from each direction. The candidate stations form the final set of upstream witness stations. and downstream witness stations .

[0070] in, This indicates the maximum number of witness stations that can be selected in each direction of propagation; Indicates the target rain gauge station The final set of upstream witness stations adopted; Indicates the target rain gauge station The final set of downstream witness stations adopted.

[0071] To ensure sufficient observational support for the bidirectional propagation constraints, the number of witness stations in both directions should meet the following requirements:

[0072] in, and These represent the number of upstream witness stations and the number of downstream witness stations, respectively. This indicates the minimum number of witness stations required for each direction of propagation.

[0073] In this embodiment, up to 4 witness stations can be selected in each direction, and at least 2 valid witness stations are required in each direction. When the number of witness stations in any direction is insufficient, automatic repair is not performed using rainfall estimation in a single direction. Instead, the data corresponding to the target rainfall station is marked as insufficient bidirectional witness evidence.

[0074] The processed output is a set of upstream witness stations. Downstream Witness Station Collection and the witness weight corresponding to each witness station This provides the data for step S3 to generate forward and reverse rainfall estimates.

[0075] Step S3: Form bidirectional propagation constraints: Based on the upstream and downstream witness stations determined in step S2, generate forward rainfall estimates and reverse rainfall estimates respectively, and construct the bidirectional propagation constraint interval for the target rainfall station.

[0076] like Figure 3 As shown, the process of constructing the bidirectional propagation constraint interval includes the following parts.

[0077] Step S3.1: Correct the propagation phase of rainfall at the witness station.

[0078] For the target rain gauge and its witness station The time difference of arrival is calculated based on the arrival times of the two events:

[0079] in, Indicates the target rain gauge station The time of arrival of the event; Witnessing Station The time of arrival of the event; Indicates rainfall was witnessed by the station. Dissipation to target rain gauge The arrival time difference.

[0080] Witness Station Rainfall sequence corrected to target rain gauge Time phase:

[0081] in, Indicates the moment to be determined; Witnessing Station Rainfall at the moment of phase correction; Witnessing Station For the target rain gauge At any moment Corrected rainfall estimates; Witnessing Station With the target rain gauge The rainfall response correction factor between them Determined based on the proportion of total rainfall from the two rain gauges or the proportion of historical responses, when reliable historical data is lacking. A value of 1 can be temporarily set, but it is recommended to adjust it based on the elevation difference between the two stations and limit it to the preset range [0.7, 1.3] to prevent the correction coefficient from being excessively amplified or reduced; when When the sampled values ​​do not correspond to the original sampling time, linear interpolation can be performed using adjacent sampled values.

[0082] when At that time, the witness station Located upstream in the propagation process, there are usually and Therefore, the rainfall process already observed by the upstream witness station is shifted in the positive direction to the target rain gauge station. when At that time, the witness station Located downstream of the propagation, usually has and Therefore, the time of the target rain gauge station can be deduced by using the same rainfall process observed later by the downstream witness station. The amount of rainfall.

[0083] Step S3.2: Generate forward and backward estimates.

[0084] The corrected rainfall estimates from the upstream witness stations are weighted and fused using the witness weights obtained in step S2 to obtain the positive rainfall estimates:

[0085] The same processing was applied to the corrected rainfall estimates at the downstream witness stations to obtain the reverse rainfall estimates:

[0086] in, Indicates the target rain gauge station At any moment Positive rainfall estimates; This indicates a reverse rainfall estimate; This represents a weighted median calculation, which arranges rainfall estimates from smallest to largest and outputs the rainfall estimate corresponding to the point where the cumulative weight first reaches half of the total weight. This formula integrates multiple witness stations to obtain a unidirectional rainfall estimate. The weighted median reduces the impact of outliers from individual witness stations on the directional estimation results. Forward estimation uses only rainfall from upstream witness stations, and backward estimation uses only rainfall from downstream witness stations. Neither directional estimate uses the current observed rainfall from the target rain gauge station. When... The positive rainfall estimate is invalid when the output is marked as invalid; when When the output is marked as invalid, the reverse rainfall estimate is invalid.

[0087] Further calculation of the robustness of discretization within both directions:

[0088]

[0089] in, This indicates the degree of dispersion of rainfall corrected by the upstream witness station; This indicates the dispersion of the corrected rainfall at downstream witness stations; 1.4826 is a constant coefficient used to convert the weighted median absolute deviation into a robust estimate of the standard deviation; the greater the dispersion, the greater the difference in estimates given by witness stations within the same direction to the target rainfall station; for the dispersion calculation, if the corresponding set is empty after filtering, then... or Marked as invalid; such as Figure 4 As shown, the rainfall sequence at the upstream witness station shifts positively based on the propagation arrival time difference.

[0090] Step S3.3: Construct the bidirectional propagation constraint interval.

[0091] When positive rainfall estimation and reverse rainfall estimation When both are valid, calculate the two-way fusion center value and the two-way estimation difference:

[0092]

[0093] in, Indicates the target rain gauge station At any moment The bidirectional fusion center value; This represents the absolute difference between the forward estimate and the reverse estimate.

[0094] By combining the internal dispersion in both directions and the difference in bidirectional estimates, the bidirectional combined dispersion is calculated:

[0095] in, Indicates a two-way merged discrete quantity; The expansion coefficient represents the difference between the two-way estimations and is used to control the impact of the difference between the forward and reverse estimations on the width of the constraint interval. This formula integrates the upstream discreteness, downstream discreteness, and the difference between the two-way estimations to obtain the merged discreteness, and the output is used to calculate the half width of the constraint interval.

[0096] Further determine the half-width of the bidirectional propagation constraint interval:

[0097] in, This represents the half-width of the bidirectional propagation constraint interval; Indicates the absolute uncertainty parameter; This represents the relative uncertainty parameter; This represents the interval width adjustment coefficient; the formula calculates the half-width of the constraint interval, and the interval is automatically widened when the witness station estimates a large discrepancy.

[0098] The lower and upper boundaries of the bidirectional propagation constraint interval are as follows:

[0099]

[0100] in, Indicates the lower boundary of the constraint interval; This represents the upper boundary of the constraint interval; the lower boundary is the larger of 0 and the calculated value, to ensure that the lower boundary is not less than 0, satisfying the physical condition that the rainfall value is non-negative.

[0101] When the estimates of the witness stations in both directions are relatively concentrated, the constraint interval is relatively narrow; when there is a large spatial difference between the witness stations, the constraint interval automatically widens; when the forward or reverse rainfall estimate is invalid, a complete bidirectional propagation constraint interval is not formed, and automatic repair is not initiated.

[0102] The processed output is a positive rainfall estimate. Reverse rainfall estimation Two-way fusion center value And a bidirectional propagation constraint interval, which is used for consistency judgment, anomaly identification and rainfall repair in step S4.

[0103] Step S4: Perform consistency gating identification and repair: Based on step S3, first determine whether the bidirectional estimates are consistent, and then identify the anomaly type of the target rain gauge station; for high-value anomaly candidates and low-value anomaly candidates, introduce a post-hoc confirmation of the propagation process, and use the temporal continuity of rainfall propagation to verify the rationality of the high value, so as to distinguish between real local heavy rainfall and sensor anomalies, and finally decide whether to perform automatic repair based on the confirmation results.

[0104] Step S4.1: Determine whether the two-way estimates are consistent.

[0105] Based on positive rainfall estimates Reverse rainfall estimation The difference between them is used to calculate the bidirectional consistency:

[0106] in, Indicates the target rain gauge station At any moment The bidirectional consistency is used for gating decisions; This represents the absolute difference between the forward rainfall estimate and the reverse rainfall estimate; Indicates the bidirectional fusion center value; Indicates the absolute uncertainty parameter; This represents the relative uncertainty parameter.

[0107] When the forward rainfall estimate and the reverse rainfall estimate are close Smaller, bidirectional consistency Approaching 1; as the difference between the estimates in the two directions increases, It decreases accordingly.

[0108] set up Indicates at time The number of upstream witness stations that can provide valid positive estimates, i.e., the number of upstream witness stations whose corrected rainfall estimates are valid in the set of upstream witness stations; This represents the number of downstream witness stations that can provide valid reverse estimates, i.e., the number of downstream witness stations in the set whose corrected rainfall estimates are valid; then the automatic repair threshold value is:

[0109] in, Indicates the automatic repair gate value; This indicates an indicator function; it takes the value 1 if all conditions within the parentheses are true, and 0 otherwise. This indicates the preset bidirectional consistency threshold; This indicates the minimum number of valid witness stations required for each direction of propagation.

[0110] In this embodiment, Pick , Take 2; when When the estimates of the two propagation directions support each other, the gating is enabled, allowing continued anomaly detection; when At that time, the gate is closed, and automatic replacement of target rain gauge data is not allowed.

[0111] Step S4.2: Identify the anomaly type of the target rain gauge.

[0112] Target rain gauge At any moment The original observed rainfall was recorded as According to the lower boundary of the constraint interval obtained in step S3 and upper boundary Determine the abnormal indication value:

[0113] in, Indicates the target rain gauge station At any moment The abnormal indicator value; when At that time, the target rain gauge observation value is outside the two-way propagation constraint range or there is a missing measurement; when At that time, the observation value of the target rain gauge was within the two-way propagation constraint interval.

[0114] when When =0, no automatic replacement is performed, the original observed rainfall is retained, and a verification mark is output.

[0115] exist Under the given conditions, the anomaly type shall be determined as follows: when If a test is missing, it is determined to be a test failure and proceeds directly to step S4.4, the repair branch. when When the value is high, it is determined to be a candidate for high-value anomaly and sent to step S4.3 for post-verification of the propagation process; when When the value is low, it is determined to be abnormal and sent to step S4.4 to perform the low value abnormality post-verification process. when If the data is normal, the original observed rainfall will be retained.

[0116] The determination of anomaly type is based on the premise of consistent bidirectional estimation. When Even if the target rain gauge observation deviates from the current constraint range, it is not directly identified as an anomaly that can be automatically repaired, but is marked as "two-way inconsistency, retain or review".

[0117] This processing method can distinguish between monitored anomalies and actual local rainfall. When a target rain gauge shows localized intensification, and the upstream forward estimate and the downstream reverse estimate cannot provide consistent evidence, the automatic repair gate remains closed, thereby preventing the replacement of the actual heavy rainfall with the estimate from a neighboring station; for example... Figure 7 As shown, the anomaly indication value is based on whether the original rainfall data observed at the target rain gauge station is missing.

[0118] Step S4.3: Perform post-propagation verification.

[0119] This step is used for bidirectional propagation posterior confirmation and intensity matching verification of high-value anomaly candidates; if it is a low-value anomaly candidate, this step is skipped and the process proceeds to step S4.4 to perform the posterior confirmation process for low-value anomalies.

[0120] For high-value anomaly candidates, this step determines the upstream prior window and downstream a posteriori window based on the direction of rainfall propagation and the time difference between stations. It calculates the rainfall enhancement at the upstream witness station before the high value occurs at the target station and at the downstream witness station after the high value occurs at the target station. This enhancement is then matched with the excess amount of the high value at the target station. Witness stations with excessive lateral offset are excluded from this matching judgment. Combining the time delay deviation and the duration of enhancement, it is determined whether the high value at the target station has evidence of the temporal evolution of actual rainfall propagation.

[0121] Step S4.3.1: Identify high-value anomaly candidates and high-value excesses.

[0122] when =1 and At time t, the target rain gauge station q is identified as a candidate for high-value anomaly; the excess amount of high-value anomalies is calculated. :

[0123] in, This represents the high value excess of the target rain gauge q at time t, i.e., the extent by which the observed value of the target station exceeds the upper boundary of the constraint interval; this formula calculates the high value excess magnitude as the benchmark for subsequent intensity matching.

[0124] Step S4.3.2: Determine the upstream prior window and the downstream a posteriori window.

[0125] For any upstream witness station i, i∈ Based on the propagation arrival time difference between the target rain gauge q and the upstream witness station i Determine the theoretical prior central moment of upstream witness station i. :

[0126] Where t represents the moment when a candidate for a high-value anomaly appears at the target rain gauge station; This represents the time difference of arrival of rainfall from witness station i to target rain gauge q, for upstream witness stations. >0; This formula calculates the central moment when the upstream witness station should theoretically observe the increase in rainfall, which is earlier than the high value moment of the target station.

[0127] Upstream Witness Station i's Pre-existing Window It is determined to be centered on the theoretical prior central time, with The time interval is half the window length:

[0128] Output the upstream prior window, which is used in step S4.3.3 to extract rainfall features within the window.

[0129] For downstream witness stations, <0, therefore the theoretical posterior central time for:

[0130] This time is later than the target station's high value time t, which is the expected time when the rainfall propagates to the downstream witness station after the target station's high value occurs.

[0131] The posterior window of downstream witness station i It is determined to be centered on the theoretical posterior central time, and with The time interval is half the window length:

[0132] Output the downstream posterior window, which is used in step S4.3.3 to extract rainfall features within the window; such as Figure 5 As shown, the candidate time t for high-value anomalies at the target rain gauge station is given.

[0133] Step S4.3.3: Calculate the propagation posterior enhancement.

[0134] For upstream witness station i, in its prior window Internal calculation of rainfall enhancement That is, the maximum rainfall within the window period:

[0135] Simultaneously record the time corresponding to this maximum value. :

[0136] The upstream prior rainfall exceeded the enhancement threshold. The number of consecutive sampling periods is denoted as the enhancement duration. :

[0137] in, This indicates an enhanced judgment threshold, with a recommended value of 2.0 mm / h to 5.0 mm / h. The duration of the enhancement of upstream witness station i within the prior window is represented by the number of consecutive sampling periods. This formula is used to determine the length of consecutive periods during which the rainfall exceeds the threshold within the statistical window, and to determine whether the enhancement process has sufficient duration to exclude false matches caused by instantaneous noise.

[0138] For downstream witness station i, in its posterior window Internal calculation of rainfall enhancement :

[0139] Simultaneously record the time corresponding to this maximum value. :

[0140] The downstream posterior window rainfall exceeded the enhancement threshold. The number of consecutive sampling periods is denoted as the enhancement duration. :

[0141] In the above formula, This represents the rainfall at witness station i during the k-th sampling period; and These represent the increase in rainfall at the upstream and downstream witness stations, respectively; and These represent the actual times when the enhancement occurs; these outputs are used for the intensity and delay matching judgment in step S4.3.4.

[0142] Step S4.3.4: Execution intensity and latency matching judgment.

[0143] Before performing a matching check, the witness stations participating in the check are first filtered: Filter 1: Lateral offset distance Exceeding the preset threshold The witness station does not participate in the upstream and downstream matching judgments below; this constraint ensures that only witness stations located on the same propagation path as the target station participate in the post-hoc confirmation, excluding interference from stations with excessive lateral offset; here The method used is compared to step S2.1. A stricter value, because step S2.3, to ensure the number of upstream and downstream witness stations, may result in a lateral offset close to... Candidate stations at the boundary are included in the witness set, so the post-validation stage adopts stricter lateral offset constraints to ensure that the stations participating in the verification of the propagation process are indeed located on the same rain cluster path in physical space. Filter 2: For witness stations participating in the matching, it is also necessary to determine whether their window data is complete; when the proportion of valid sampling points of the witness station within the corresponding window is not less than R... min If the data is complete, the window is considered complete; otherwise, it is considered missing data, and the witness station will not participate in the matching judgment. Valid sampling points refer to points that are not missing measurements and are within a reasonable rainfall range. min The minimum percentage of valid sampling points that represents the integrity of window data is recommended to be between 0.70 and 0.85, and 0.80 is used in this embodiment.

[0144] The upstream strength matching condition is that there exists at least one upstream witness station i that meets the lateral offset constraint, and the following three inequalities are satisfied simultaneously:

[0145]

[0146]

[0147] in, This represents the intensity matching ratio threshold, 0 < Coefficients ≤1; This represents the time delay deviation threshold, expressed in sampling intervals. The lower limit protection clause: to prevent damage caused by light rainfall or background noise. When the rainfall is very small, such as only 0.5 mm, even a trace amount of rainfall at the witness station can meet the proportionality condition and incorrectly pass the posterior confirmation; for low-value abnormal candidates, the above formula will be used to... Replace with the low value excess, which is the magnitude by which the target station's observed value is lower than the lower boundary of the constraint interval.

[0148] The first inequality determines whether the rainfall enhancement at the upstream witness station within the prior window reaches a certain proportion of the target station's high value, thus confirming whether the target station's high value may be caused by rain clusters propagating from upstream. The third inequality determines whether the actual time of the upstream enhancement is basically consistent with the theoretical propagation time, thus eliminating temporally irrelevant rainfall enhancement interference. The above conditions together determine whether the upstream direction provides effective post-propagation verification evidence.

[0149] The downstream strength matching condition is that there exists at least one downstream witness station i that meets the lateral offset constraint, and the following three inequalities are satisfied simultaneously:

[0150]

[0151]

[0152] The meanings of each symbol are consistent with those in the upstream intensity matching conditions. The first inequality determines whether the rainfall enhancement at the downstream witness station within the posterior window reaches a certain proportion of the target station's high value excess. The second inequality determines whether the actual time of downstream enhancement is basically consistent with the theoretical propagation time. These two conditions together determine whether the downstream direction provides effective post-propagation verification evidence.

[0153] Step S4.3.5: Determine the post-propagation confirmation result.

[0154] Define the propagation posterior confirmation indicator value It shall be determined according to the following rules: Determine whether the sum of the number of valid upstream witness stations and the number of valid downstream witness stations that meet the horizontal constraints and have complete window data is less than n. p If so, then it is deemed insufficient evidence. Setting the value to -1 disables automatic repair and marks the target rain gauge data as being in a verification state; where n p This represents the minimum number of valid witness stations required to determine that the evidence is sufficient; a value of 4 is recommended. If the total number of valid witness stations is not less than n p Further judgment: When, after filtering, at least one witness station with complete window data is found among the witness stations that meet the lateral constraints, but neither the upstream strength matching condition nor the downstream strength matching condition is met, then... A value of 0 indicates that the post-propagation confirmation failed. In this case, the high value of the target rain gauge is considered to lack corresponding temporal evolution evidence in the rainfall propagation chain, and is more likely to be an isolated anomaly caused by sensor failure, communication interference, or other non-rainfall factors, allowing automatic repair to be performed.

[0155] When the upstream strength matching condition is met or the downstream strength matching condition is met A value of 1 indicates successful post-propagation confirmation. At this point, the high value at the target rain gauge is considered to provide evidence of the temporal evolution of actual rainfall propagation. This high value represents a normal observation of genuine localized enhanced rainfall or moving rain clusters, and automatic repair should not be directly implemented. Figure 6 As shown, this refers to the proportion of rainfall increase within the upstream prior window that exceeds the high value.

[0156] Step S4.4: Generate the repair rainfall value.

[0157] Target rain gauge At any moment The amount of rainfall used for repair is recorded as , When low-value abnormal indicator variable When =1, the repair is not performed directly, but instead the process is switched to the low-value anomaly post-verification confirmation sub-process; it is necessary to determine whether there is a significant increase in rainfall at the corresponding propagation time upstream or downstream. If there is significant rainfall upstream and downstream, indicating an anomaly with the target station's value being lower than normal: Confirmation passed, bidirectional fusion of central values ​​adopted. Repair the data from the target rain gauge station; If there is no rainfall upstream or downstream, meaning the target station's rainfall has decreased: the confirmation fails, automatic repair is not performed, and the original rainfall observation at the target rain gauge station is retained. It also outputs a verification flag.

[0158] For high-value anomaly candidates, confirm the indication value based on the propagation posterior. Determine the handling method: when =0, meaning that automatic repair is performed when the post-propagation verification fails. = .

[0159] when =1, meaning that if the post-propagation verification is successful, automatic repair will not be performed, and the original observed rainfall of the target rain gauge will be retained. It also outputs a flag indicating suspicious high values ​​for verification.

[0160] when =-1, meaning that when the post-propagation confirmation evidence is insufficient, automatic repair is not performed, and the original observed rainfall of the target rain gauge is retained. It also outputs a verification flag.

[0161] As an optional implementation, for low-value anomalies, a similar posterior verification logic can also be selectively applied: calculating the amount of low value exceeding the limit. : ; To constrain the lower boundary of the interval, the upstream antecedent window and downstream posterior window are determined in a manner similar to that used for high-value anomalies, and the upstream enhancement is calculated. and downstream enhancement If there is a significant increase in rainfall upstream or downstream at the corresponding propagation time, and the time delay deviation meets the requirements and the enhancement duration reaches the target, while the observation value at the target station... If the value is significantly lower than normal, it is more likely that the sensor value is abnormally low and an automatic repair should be performed. If there is no corresponding increase in rainfall in either the upstream or downstream, it may be due to a decrease in actual rainfall or localized lack of rain. In this case, a forced repair should not be performed, and the original observation should be retained.

[0162] This invention supports two execution modes.

[0163] Offline rainfall event analysis mode: After all rainfall events have ended, the calculation is performed in a unified manner. The complete time series data of all stations are ready. There is no need to set a waiting window. The complete process of steps S1 to S4 is executed directly.

[0164] Real-time data stream processing mode: Set downstream confirmation wait window. Wait window duration. Set as the maximum propagation time difference between stations Times:

[0165] in, A value of 1.2 to 1.5 is recommended. For the set of downstream witness stations of target station q, Let q be the propagation arrival time difference between the target station q and the downstream witness station j.

[0166] exist During this period, the suspected abnormal data at the target station remained in a pending confirmation state. Even if the automatic repair gate is calculated =1, and the exception type determination and subsequent repair logic in step S4.2 are not executed; The use of one-way rainfall estimation to automatically replace the target station's observed rainfall is prohibited; Data awaiting confirmation will not proceed to the post-propagation confirmation step S4.3.

[0167] when After timeout: If the rainfall data corresponding to the downstream witness station has been received completely, the pending confirmation state will be exited, and bidirectional consistency judgment, anomaly identification and post-confirmation will be performed. If downstream data still hasn't arrived, it will be treated as insufficient evidence. =-1, do not perform automatic repair, retain the original observed rainfall and output the verification mark.

[0168] like Figure 8 As shown, in both offline rainfall event analysis mode and real-time data stream processing mode, through steps S1 to S4, the present invention finally outputs the original observed rainfall, anomaly type marker, propagation posterior confirmation status, repaired rainfall value and verification marker at each moment of the target rain gauge station, thus completing the identification and repair of rainfall anomalies based on bidirectional propagation constraints and propagation process posterior confirmation.

[0169] It should be noted that the parameters involved in this invention can be adjusted according to the density of rain gauges, temporal resolution, and regional rainfall characteristics. For example, a preset arrival ratio can be used. Recommended value range is 0.3 to 0.7; maximum distance along the propagation direction. Recommended range: 30 km to 150 km; maximum lateral offset distance Recommended distance is 10 km to 50 km; distance attenuation scale along the propagation direction. Recommended range: 20 km to 80 km; Lateral distance attenuation scale Recommended distance is 10 km to 40 km; maximum number of witness stations allowed per propagation direction. Recommended value: 3 to 5; minimum number of effective witness stations required for each propagation direction. Recommended value: 2; Absolute uncertainty parameter A value of 0.1 mm to 0.5 mm is recommended. Relative uncertainty parameter. Recommended values ​​are 0.05 to 0.15; the interval width adjustment coefficient λ is recommended to be 1.5 to 2.5; the two-way estimation difference expansion coefficient γ is recommended to be 0.1 to 0.4; and a preset two-way consistency threshold is recommended. Recommended value: 0.4 to 0.7; Upstream antecedent half-window length It is recommended to use 2 to 6 sampling intervals; downstream posterior half-window length It is recommended to use 2 to 6 sampling intervals; intensity matching ratio threshold. Recommended value: 0.3 to 0.8; Delay deviation threshold It is recommended to use 1 to 3 sampling intervals.

[0170] In a second embodiment, the present invention also discloses a computer-readable storage medium storing a computer program, which, when executed by a processor, causes the processor to perform the steps of the method described above.

[0171] In embodiment three, the present invention also discloses a computer device, including a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor performs the steps of the method described above.

[0172] In another embodiment provided in this application, a computer program product containing instructions is also provided, which, when run on a computer, causes the computer to perform any of the methods described in the above embodiments.

[0173] It is understood that the systems, devices, and storage media provided in the embodiments of the present invention correspond to the methods provided in the embodiments of the present invention, and the explanations, examples, and beneficial effects of the relevant content can be referred to the corresponding parts of the above methods.

[0174] In the above embodiments, implementation can be achieved entirely or partially through software, hardware, firmware, or any combination thereof. When implemented using software, it can be implemented entirely or partially in the form of a computer program product. The computer program product includes one or more computer instructions. When the computer program instructions are loaded and executed on a computer, all or part of the processes or functions described in the embodiments of this application are generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. The computer instructions can be stored in a computer-readable storage medium or transmitted from one computer-readable storage medium to another. For example, the computer instructions can be transmitted from one website, computer, server, or data center to another website, computer, server, or data center via wired or wireless means. The computer-readable storage medium can be any available medium accessible to a computer or a data storage device such as a server or data center that integrates one or more available media. The available medium can be a magnetic medium, an optical medium, or a semiconductor medium, etc.

[0175] It should be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element.

[0176] The various embodiments in this specification are described in a related manner. Similar or identical parts between embodiments can be referred to mutually. Each embodiment focuses on describing the differences from other embodiments. In particular, the system embodiments are basically similar to the method embodiments, so the description is relatively simple; relevant parts can be referred to the descriptions of the method embodiments.

[0177] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for identifying and repairing rainfall anomalies based on bidirectional propagation, characterized in that, include: Based on the rainfall data of each rain gauge station within a rainfall event, the arrival time of the event at each rain gauge station is determined. Plane fitting is performed based on the spatial location of the rain gauge station and the arrival time of the event to determine the direction of rainfall propagation. Based on the projected position of the rain gauge station in the direction of rainfall propagation, the propagation arrival sequence of the rain gauge stations is established. Based on the propagation arrival order of rain gauges, the candidate rain gauges are constrained by the propagation direction and the lateral offset distance to obtain the upstream witness station set and the downstream witness station set. The propagation arrival time difference is determined based on the arrival times of events at the witness station and the target rain gauge station. The phase correction is performed on the rainfall sequence of the witness station, and the forward rainfall estimate and the reverse rainfall estimate are generated by fusion respectively. The bidirectional propagation constraint interval of the target rain gauge station is constructed by combining the internal dispersion. Based on the bidirectional propagation constraint interval and forward and reverse rainfall estimates, the anomaly types of the target rain gauge are identified, distinguishing between missing anomalies, low-value anomalies, and high-value anomaly candidates. For missing anomalies, direct fusion and repair are performed. For low-value anomalies, repair or retention is performed after confirmation based on upstream and downstream rainfall increments. For high-value anomaly candidates, the a priori and posterior windows of the witness station are determined based on the propagation arrival time difference between the target rain gauge and the witness station. Within each window, the rainfall enhancement is calculated and intensity and time delay matching judgment is performed. The propagation posterior confirmation result is obtained based on the judgment result. Based on the propagation posterior confirmation result and the fusion value of forward and reverse rainfall estimates, the repaired rainfall value is generated.

2. The method for identifying and repairing rainfall anomalies based on bidirectional propagation according to claim 1, characterized in that, The process of determining the arrival time of an event at a rain gauge station includes: calculating the cumulative rainfall and the total rainfall of the rainfall event at each rain gauge station up to the current sampling time; determining the moment when the cumulative rainfall first reaches a preset proportion of the total rainfall of the rainfall event as the arrival time of the event; and when the total rainfall of the rainfall event is lower than a preset minimum rainfall threshold, the corresponding rain gauge station does not participate in the subsequent propagation direction calculation.

3. The method for identifying and repairing rainfall anomalies based on bidirectional propagation according to claim 1, characterized in that, The process of determining the direction of rainfall propagation in a rainfall event includes: using the least squares method to fit a propagation arrival time plane and calculating the fitting residual for the arrival time of each rain gauge station; identifying and excluding rain gauge stations that deviate from the propagation arrival time plane based on the absolute deviation of the median of the residuals, and then refitting; determining the direction of rainfall propagation based on the spatial gradient of the refitted plane, and arranging each rain gauge station in ascending order of its projection position in the propagation direction to obtain the order of propagation arrival.

4. The method for identifying and repairing rainfall anomalies based on bidirectional propagation according to claim 1, characterized in that, The process of obtaining the upstream witness station set and the downstream witness station set includes: Calculate the distance along the propagation direction and the lateral offset distance between the candidate rain gauge and the target rain gauge; Candidate rain gauges that meet the preset distance thresholds for both the distance along the propagation direction and the lateral offset distance are selected and used as the upstream candidate witness station set and the downstream candidate witness station set, respectively. The witness weight is calculated for the candidate witness stations according to the distance attenuation method. The candidate stations are sorted according to the witness weight and a preset number of upper limit is selected to obtain the upstream witness station set and the downstream witness station set respectively. If the number of either set is less than the preset minimum number, automatic repair is not performed.

5. The method for identifying and repairing rainfall anomalies based on bidirectional propagation according to claim 4, characterized in that, The process of performing phase correction on the rainfall sequence of the witness station includes: The propagation arrival time difference is calculated based on the arrival times of the events at the target rain gauge station and the witness station. The rainfall sequence at the witness station is shifted to the time phase of the target rain gauge station according to the propagation arrival time difference. The shifted rainfall is then corrected according to the ratio of the total rainfall of the events at the two stations or the ratio of the historical response. The fusion process for generating both forward and reverse rainfall estimates does not use the observed rainfall at the target rain gauge station at the current moment.

6. The method for identifying and repairing rainfall anomalies based on bidirectional propagation according to claim 4, characterized in that, The process of fusing and generating positive and negative rainfall estimates includes: When the witness station is located upstream of the propagation, the rainfall sequence of the upstream witness station is forward-shifted by the time difference of arrival to obtain the upstream corrected rainfall estimate; When the witness station is located downstream of the propagation, the propagation arrival time difference is back-calculated using the rainfall sequence of the downstream witness station to obtain the downstream corrected rainfall estimate; The upstream and downstream corrected rainfall estimates are fused using a weighted median fusion method based on witness weights to obtain forward and reverse rainfall estimates. The internal dispersion of each direction is determined based on the deviation between the internal corrected rainfall and the corresponding estimated value.

7. The method for identifying and repairing rainfall anomalies based on bidirectional propagation according to claim 1, characterized in that, The process of constructing the bidirectional propagation constraint interval of the target rain gauge includes: Calculate the two-way fusion center value and the two-way estimation difference between the forward and reverse rainfall estimates; The two-way combined dispersion is calculated by combining the upstream dispersion, downstream dispersion, and the two-way estimation difference; Based on the bidirectional merging discrete quantity, absolute uncertainty parameter, and relative uncertainty parameter, the half-width of the constraint interval is calculated, and the lower and upper boundaries of the bidirectional propagation constraint interval are obtained by combining the bidirectional fusion center value.

8. The method for identifying and repairing rainfall anomalies based on bidirectional propagation according to claim 1, characterized in that, The process of identifying the anomaly type of the target rain gauge includes: Based on forward and reverse rainfall estimates, the bidirectional consistency is calculated, and the automatic repair threshold is determined by combining the number of valid witness stations. When the automatic repair threshold value meets the preset conditions, it is determined whether the original observed rainfall of the target rain gauge exceeds the two-way propagation constraint range. When the original observed rainfall is missing, it is identified as a missing measurement anomaly; when it exceeds the upper boundary of the two-way propagation constraint interval, it is identified as a high-value anomaly candidate; when it is below the lower boundary of the two-way propagation constraint interval, it is identified as a low-value anomaly.

9. The method for identifying and repairing rainfall anomalies based on bidirectional propagation according to claim 1, characterized in that, The process of calculating the rainfall enhancement and performing intensity-time delay matching includes: The magnitude by which the target station's observation exceeds the upper boundary of the two-way propagation constraint interval is calculated to obtain the high-value excess. Based on the propagation arrival time difference between the target rain gauge station and the upstream witness station, the theoretical prior center time is determined, and the maximum rainfall within the upstream prior window is taken as the upstream rainfall enhancement amount. Based on the propagation arrival time difference between the target rain gauge station and the downstream witness station, the theoretical posterior center time is determined, and the maximum rainfall within the downstream posterior window is taken as the downstream rainfall enhancement amount. Witness stations whose lateral offset distance exceeds the preset offset threshold are filtered out. The proportion of upstream rainfall enhancement to the excess of high value and peak delay deviation are calculated. The proportion of downstream rainfall enhancement to the excess of high value and peak delay deviation are calculated. If either proportion or delay deviation meets the corresponding matching condition, the propagation post-hoc confirmation is determined to be successful and the original observation is retained. Otherwise, the propagation post-hoc confirmation is determined to be unsuccessful and repaired with the fused value.

10. A method for identifying and repairing rainfall anomalies based on bidirectional propagation according to claim 5, characterized in that, The process of generating the repaired rainfall value includes: For low-value anomalies, determine whether there is an increase in rainfall at the corresponding propagation time upstream or downstream. If there is an increase in rainfall upstream or downstream, the low value is confirmed as an anomaly and repaired using bidirectional fusion of central values. If there is no increase in rainfall upstream or downstream, the low value is confirmed as a reduction in actual rainfall and the original observed rainfall is retained. For high-value anomaly candidates, if the post-propagation confirmation result is not passed, the two-way fusion center value is used for repair; if the post-propagation confirmation result is passed or the evidence is insufficient, the original observed rainfall is retained and a verification mark is output. For missing test anomalies, the bidirectional fusion center value is used directly for repair.

Images