Particle size pollution diffusion simulation method and system based on sailing reconstruction

Abstract

The application discloses a particle size pollution diffusion simulation method and system based on sailing reconstruction, and belongs to the technical field of atmospheric environment monitoring and pollution diffusion simulation. The particle size detection, positioning, inertia, meteorological and sampling chain working condition data are acquired, the inlet sampling efficiency and the discrete delay-broadening transfer kernel are determined according to the particle size, and the particle size space-time observation zone is constructed. The reverse Lagrange calculation is performed in the reconstructed wind field to form a particle size source-receptor matrix, the particle size source emission rate is inversed under the non-negative, spatial sparsity and cross-particle diameter common source activation constraints, and the concentration field and the settling simulation result are obtained by the forward particle diffusion. Through closed-loop calibration, residual discrimination, local update and rollback, the position deviation, concentration tailing and space-time aliasing are reduced, and the stability of the particle size composition, the settling area and the pollution boundary simulation is improved.

Description

Technical Field

[0001] This invention relates to the field of atmospheric environmental monitoring and pollution diffusion simulation technology, and in particular to a method and system for simulating particulate matter size pollution diffusion based on mobile reconstruction. Background Technology

[0002] Mobile particulate matter monitoring can simultaneously acquire multi-size particle concentration information at different locations during movement, and combine it with location and meteorological data to describe pollution distribution. Current processing typically uses the timestamp and vehicle location provided by the detection equipment directly as the sampling time and location, and then conducts pollution source location or diffusion simulation based on the concurrent wind field.

[0003] However, the sampling efficiency, transmission loss, response delay, and broadening degree of particles of different sizes vary at the sampling inlet, sampling pipeline, and detection equipment; the sampling frequency and time delay of particle size detection equipment, positioning equipment, inertial equipment, and meteorological sensors are also inconsistent. Direct registration will cause the detected values ​​to mismatch with the actual sampling time, actual sampling location, and corresponding wind field, resulting in positional offset, concentration tailing, and spatiotemporal aliasing.

[0004] Furthermore, particles of different sizes exhibit differences in advection, turbulent diffusion, gravity settling, and surface deposition. If source intensity inversion and forward simulation are not coupled with the sampling chain, migration trajectory, wind field, and particle size transport using the same observation relationship, the pollution center, diffusion direction, particle size composition, settling range, and pollution boundary are prone to distortion. Indiscriminately updating the sampling chain parameters based on the model residuals may also erase the true simultaneous pollution peaks of multiple particle sizes.

[0005] Therefore, there is a need for a method and system for simulating particulate matter size pollution diffusion based on mobile reconstruction that can overcome the shortcomings of the existing technologies. Summary of the Invention

[0006] The purpose of this invention is to provide a method and system for simulating particulate matter size-based pollution diffusion based on mobile reconstruction. This method uses an inlet sampling efficiency determined by particle size and a discrete delay-broadening transfer kernel to reconstruct the detected values ​​into a spatiotemporal observation band of particle size distributed across multiple historical actual sampling times, locations, and corresponding wind fields. A particle size source-receptor matrix is ​​established through inverse Lagrange calculations incorporating advection, turbulent random diffusion, gravity settling, and surface deposition. The particle size source emission rate is inverted under constraints of non-negativity, spatial sparsity, and cross-size common source activation, thereby forming a concentration field and sedimentation simulation results from forward particle diffusion. Simultaneously, a closed loop is constructed using co-pulse calibration, particle size-wise residual discrimination, local updates, and acceptance rollback.

[0007] This invention provides a method for simulating particulate matter particle size pollution diffusion based on mobile reconstruction, comprising:

[0008] S1. Acquire particle size detection data, positioning data, inertial data, meteorological data and sampling chain operating data of the mobile carrier, determine the inlet sampling efficiency and the discrete delay-broadening transfer kernel formed by the combination of sampling inlet, sampling pipeline and detection equipment according to particle size, and generate the reconstructed wind field of the mobile area according to the meteorological data, positioning data and inertial data.

[0009] S2. Based on the inlet sampling efficiency and discrete delay-broadening transfer kernel, each particle size detection value is characterized as the convolution result of the environmental concentration at multiple historical actual sampling times weighted by kernel weights. Each kernel delay unit is mapped to the corresponding historical actual sampling time's navigation position, heading, and local wind field in the reconstructed wind field to obtain the particle size spatiotemporal observation band.

[0010] S3. In the reconstructed wind field, perform inverse Lagrange calculations for each particle size, including advection, turbulent random diffusion, gravity settling and surface deposition. Incorporate the kernel weights of the particle size spatiotemporal observation zone into the transmission contribution from the candidate source grid to the particle size detection value to obtain the particle size source-receptor matrix.

[0011] S4. Based on the particle size source-acceptor matrix and the detection values ​​of each particle size, estimate the particle size source emission rate vector of each candidate source grid under non-negativity constraints, spatial sparsity constraints and cross-particle size common source activation constraints.

[0012] S5. Perform forward particle diffusion including advection, turbulent random diffusion, gravity settling and surface deposition based on the particle size source emission rate vector to obtain the particle size concentration field and particulate matter settling simulation results.

[0013] Optionally, S1 includes:

[0014] A low-hysteresis bypass, a main sampling path, a switching component for switching between the low-hysteresis bypass and the main sampling path, and a multi-particle size calibration pulse component are set at the determined location of the discrete delay-broadening transfer core.

[0015] The same calibration pulse is applied to the low-hysteresis bypass and the main sampling path by the multi-particle size calibration pulse component. The low-hysteresis bypass output sequence is used as the input reference and the main sampling path output sequence is used as the response output. The kernel mean and kernel covariance of the discrete delay-broadening transfer kernel are identified according to the particle size.

[0016] The inlet sampling efficiency is determined based on the ratio of the integral response of the main sampling path output sequence to the integral response of the low-hysteresis bypass output sequence, and the kernel mean, kernel covariance, inlet sampling efficiency, and corresponding sampling chain operating conditions are written into the particle size calibration record.

[0017] Optionally, S2 includes:

[0018] The detection time is determined by using the timestamp of the particle size detection data as the detection time, and the corresponding historical actual sampling time is determined according to the delay amount of each delay unit in the kernel mean value;

[0019] Based on the timestamps of the positioning data, inertial data, and meteorological data, time interpolation is performed on the data before and after the actual historical sampling time to obtain the corresponding navigation position, heading, and local wind field;

[0020] The historical actual sampling time, navigation position, heading, local wind field, particle size identifier, and kernel weight are combined to form a particle size spatiotemporal observation band, and the observation uncertainty of each observation item in the particle size spatiotemporal observation band is determined by the kernel covariance.

[0021] Optionally, S3 includes:

[0022] Release the opposite particle of the corresponding particle size from each travel position in the particle size spatiotemporal observation zone, and update the horizontal position of the opposite particle in the reconstructed wind field according to the opposite direction of the local wind speed.

[0023] Random displacements are generated based on the turbulence statistics of the reconstructed wind field. Gravity settling velocity is determined based on particle size, particle density, and air properties. The vertical position of the reverse particles is updated in the opposite direction of gravity settling according to the gravity settling velocity.

[0024] The retention weight of reverse particles is determined based on the deposition rate corresponding to the surface type.

[0025] The retention weights of the reverse particles entering each candidate source grid are accumulated with the kernel weights of the particle size spatiotemporal observation band to form the matrix elements of the particle size source-receptor matrix.

[0026] Optionally, S4 includes:

[0027] A data fitting term is constructed by the weighted residual between each particle size detection value and the product of the particle size source-receptor matrix and the particle size source emission rate vector, wherein the weight of each observation term is determined by the inverse value of the corresponding observation uncertainty;

[0028] Spatial sparsity penalty is applied to the emission rates of particle size sources in adjacent candidate source grids for the same particle size, and common source group penalty is applied to the emission rates of multiple particle size sources in the same candidate source grid. The sum of data fitting term, spatial sparsity penalty and common source group penalty is solved under the condition that the particle size source emission rate is non-negative, and the particle size source emission rate vector is obtained.

[0029] Optionally, S5 includes:

[0030] Forward particles are released from the corresponding candidate source grid according to the particle size source emission rate vector, and the particle size, position, mass and deposition state of each forward particle are recorded.

[0031] The mass of the undeposited forward particles is accumulated according to the spatial grid to obtain the concentration field of each particle size. The total particulate matter concentration field is obtained by summing the concentration fields of each particle size. The pollution center is determined by the mass-weighted position of the total particulate matter concentration field. The diffusion direction is determined by the main axis direction of the concentration distribution in the neighborhood of the pollution center. The particle size composition is determined by the proportion of each particle size concentration to the total particulate matter concentration.

[0032] The mass of the forward-deposited particles is accumulated according to the surface grid and time interval to obtain the sedimentation flux and sedimentation range. The pollution boundary is determined by the outer envelope of the spatial grid when the total particulate matter concentration reaches the boundary threshold. The boundary threshold is the preset quantile of the concentration sample in the background area of ​​the mobile survey.

[0033] Furthermore, it also includes: based on the particle size spatiotemporal observation zone, reading the environmental prediction concentration at each historical actual sampling time and navigation position in the particle size spatiotemporal observation zone from the particle size concentration field obtained by forward particle diffusion, and weighting the environmental prediction concentration with the corresponding kernel weight and correcting it with the inlet sampling efficiency to obtain the prediction detection response, and calculating the particle size residual based on each particle size detection value and the corresponding prediction detection response.

[0034] Within a discrimination window whose length is not less than the kernel support time of the corresponding discrete delay-broadened transfer kernel, the residual sign ratio is calculated, the correlation coefficient between the residual of the target particle size and the operating parameters of each sampling chain is calculated one by one, and one of the absolute values ​​of each correlation coefficient with a value not less than other absolute values ​​is selected as the operating condition correlation, and the multi-particle-size synchronous pollution peak identifier is determined at the same time.

[0035] When the proportion of residuals with the same sign of the target particle size reaches the first quantile threshold determined by the uncontaminated calibration residual sample, the working condition correlation reaches the second quantile threshold determined by the working condition correlation of the stable working condition calibration sample, and the multi-particle-size synchronous contamination peak identifier is not specified, a sampling chain drift identifier for the target particle size is generated.

[0036] The background concentration quantile thresholds for each particle size are determined by the preset quantiles of the corresponding particle size detection values ​​in the mobile background area. The equivalent actual sampling time of each particle size detection peak is determined by weighting the relevant historical actual sampling times with the kernel weight of the corresponding particle size spatiotemporal observation zone. The multi-particle-size synchronous pollution peak identifier is determined based on the condition that the detection peaks of at least two particle sizes are higher than their respective background concentration quantile thresholds and the difference between the equivalent actual sampling times is not greater than the sampling period of the positioning data.

[0037] Furthermore, the processing of the sampling chain drift marker includes: for each particle size residual sample within the discrimination window, the local eastward wind speed component and the local northward wind speed component at each historical actual sampling time are weighted by the kernel weight of the particle size spatiotemporal observation zone corresponding to the residual sample, so as to obtain the equivalent eastward wind speed component and the equivalent northward wind speed component.

[0038] Calculate the correlation coefficients between the residuals of each particle size and the corresponding equivalent eastward wind speed component and equivalent northward wind speed component, and select the one whose value is not less than the other absolute value as the residual-wind field correlation of the corresponding particle size.

[0039] When the target particle size generates a sampling chain drift flag, a low-hysteresis bypass is connected through the switching component and a calibration pulse is applied by the multi-particle size calibration pulse component. The discrete delay-broadening transfer kernel and kernel covariance of the target particle size are updated according to the new input reference and response output, while the discrete delay-broadening transfer kernel of the other particle sizes remains unchanged.

[0040] When no sampling chain drift identifier is generated and the multi-size synchronous pollution peak identifier is set to yes, the discrete delay-broadening transfer kernel of each particle size remains unchanged and the particle size source emission rate vector is re-estimated.

[0041] When no sampling chain drift identifier is generated, the multi-particle size synchronous pollution peak identifier is not generated, and the residual-wind field correlation of at least two particle sizes reaches the third quantile threshold determined by the stable wind field sample, the discrete delay-broadening transfer kernel of each particle size remains unchanged and the wind field is reconstructed according to the residual correction.

[0042] When no sampling chain drift identifier is generated, the multi-particle size synchronous pollution peak identifier is not present, and the particle size corresponding to the residual-wind field correlation that reaches the third quantile threshold is less than two, keep the discrete delay-broadening transfer kernel, particle size source emission rate vector and reconstructed wind field of each particle size unchanged, record the state to be determined and enter the next discrimination window.

[0043] Furthermore, after updating the discrete delay-broadening transfer kernel and kernel covariance of the target particle size, the spatiotemporal observation band of the target particle size is reconstructed using the discrete delay-broadening transfer kernel before and after the update, and inversion and forward particle diffusion are performed.

[0044] When the updated sum of squared residuals of the discrimination window is less than the sum of squared residuals of the discrimination window before the update, and the sum of the diagonal elements of the updated kernel covariance is not greater than the sum of the diagonal elements of the kernel covariance before the update, the updated discrete delay-widened propagation kernel and kernel covariance are received; otherwise, the process reverts to the previous discrete delay-widened propagation kernel and kernel covariance, and the values ​​of both before and after the update, as well as the receiving or reverting result, are written into the particle size calibration record.

[0045] On the other hand, the present invention also provides a particulate matter particle size pollution diffusion simulation system based on mobile reconstruction, comprising:

[0046] The sampling chain identification module acquires particle size detection data, positioning data, inertial data, meteorological data, and sampling chain operating condition data, and determines the inlet sampling efficiency and discrete delay-broadening transfer kernel according to particle size. The wind field reconstruction module generates a reconstructed wind field for the underway area based on the meteorological, positioning, and inertial data. The spatiotemporal observation zone construction module maps each kernel delay unit to the underway position, heading, and local wind field in the reconstructed wind field at the historical actual sampling time, obtaining the particle size spatiotemporal observation zone. The source-receptor matrix construction module performs inverse Lagrange multiplication to obtain the particle size source-receptor matrix. The source emission rate inversion module estimates the particle size source emission rate vector under non-negativity constraints, spatial sparsity constraints, and cross-particle size common source activation constraints. The diffusion simulation module... The first block is used to perform forward particle diffusion and obtain the concentration fields of each particle size, the total particulate matter concentration field, the pollution center, the diffusion direction, the particle size composition, the sedimentation flux, the sedimentation range, and the pollution boundary. The second block is a closed-loop calibration module, including a low-hysteresis bypass, a main sampling path, a switching component, and a multi-particle-size calibration pulse component. It is used to identify the discrete delay-broadening transfer kernel and the kernel covariance of each particle size. The forward particle diffusion results are processed by the inlet sampling efficiency, the discrete delay-broadening transfer kernel, and the particle size spatiotemporal observation band to obtain the predicted detection response. The sampling chain drift is judged according to the particle-size residual of the predicted detection response. The discrete delay-broadening transfer kernel of the target particle size, the particle size source emission rate vector, or the reconstructed wind field are updated according to the judgment results. The update is received or the rollback is performed according to the residual and kernel covariance before and after the update.

[0047] The beneficial effects of this invention are:

[0048] 1. By mapping historical trajectories and reconstructing wind fields item by item using particle size-related inlet sampling efficiency and discrete delay-broadened transfer kernel, the mismatch between detected values ​​and actual sampling time, location and wind field can be corrected, reducing location offset, concentration tail and spatiotemporal aliasing.

[0049] 2. By coupling reverse particle size transport, source strength constraint inversion, and forward particle diffusion under the same observation operator, the advection, turbulent random diffusion, gravity sedimentation, and surface deposition of various particle sizes are treated consistently, thereby improving the reliability of concentration field, particle size composition, sedimentation range, and pollution boundary.

[0050] 3. Differentiate between real pollution, sampling chain drift, and wind field deviation by using multi-particle-size synchronous pollution peaks, residual sign similarity, and sampling chain operating condition correlation. Simultaneously accept or backtrack the transfer kernel and kernel covariance to avoid erasing real pollution peaks and improve the stability and traceability of closed-loop updates. Attached Figure Description

[0051] The accompanying drawings are provided to further illustrate the invention and form part of the specification. They are used in conjunction with embodiments of the invention to explain the invention and do not constitute a limitation thereof. In the drawings:

[0052] Figure 1 This is a flowchart of the particulate matter size pollution diffusion simulation method based on mobile reconstruction according to the present invention.

[0053] Figure 2 This is a flowchart of the S4 multi-constraint source emission rate inversion of the present invention. Detailed Implementation

[0054] The present invention will now be described in further detail with reference to the accompanying drawings. These drawings are simplified schematic diagrams, illustrating only the basic structure of the invention, and therefore only show the components relevant to the invention.

[0055] refer to Figures 1-2 A method for simulating particulate matter particle size pollution diffusion based on mobile reconstruction includes:

[0056] S1. Acquire particle size detection data, positioning data, inertial data, meteorological data and sampling chain operating data of the mobile carrier, determine the inlet sampling efficiency and the discrete delay-broadening transfer kernel formed by the combination of sampling inlet, sampling pipeline and detection equipment according to particle size, and generate the reconstructed wind field of the mobile area according to the meteorological data, positioning data and inertial data.

[0057] S2. Based on the inlet sampling efficiency and discrete delay-broadening transfer kernel, each particle size detection value is characterized as the convolution result of the environmental concentration at multiple historical actual sampling times weighted by kernel weights. Each kernel delay unit is mapped to the corresponding historical actual sampling time's navigation position, heading, and local wind field in the reconstructed wind field to obtain the particle size spatiotemporal observation band.

[0058] S3. In the reconstructed wind field, perform inverse Lagrange calculations for each particle size, including advection, turbulent random diffusion, gravity settling and surface deposition. Incorporate the kernel weights of the particle size spatiotemporal observation zone into the transmission contribution from the candidate source grid to the particle size detection value to obtain the particle size source-receptor matrix.

[0059] S4. Based on the particle size source-acceptor matrix and the detection values ​​of each particle size, estimate the particle size source emission rate vector of each candidate source grid under non-negativity constraints, spatial sparsity constraints and cross-particle size common source activation constraints.

[0060] S5. Perform forward particle diffusion including advection, turbulent random diffusion, gravity settling and surface deposition based on the particle size source emission rate vector to obtain the particle size concentration field and particulate matter settling simulation results.

[0061] In this specific embodiment, S1 includes:

[0062] In this specific embodiment, the mobile monitoring vehicle is an electric monitoring vehicle equipped with a multi-channel optical particle size spectrometer, a dual-antenna satellite positioning unit, a six-axis inertial measurement unit, and a three-dimensional ultrasonic meteorological instrument. The three-dimensional ultrasonic meteorological instrument synchronously outputs three wind speed components in the instrument coordinate system. The particle size spectrometer divides the particle size into five grades: 0.3–0.5 μm, 0.5–1 μm, 1–2.5 μm, 2.5–5 μm, and 5–10 μm, with a detection cycle of 1 second. Each particle size record is recorded with optical equivalent diameter, standard polystyrene particle response curve, target particle refractive index, shape factor, density calibration batch, and diameter conversion version. The controller converts the optical equivalent diameter into the volume equivalent physical diameter for sedimentation calculation. The acquisition controller records particle size detection records, positioning records, inertial records, meteorological records, and sampling chain operating condition records with a uniform monotonic clock. The sampling chain operating condition records include at least the sampling pump flow rate, pipeline temperature, pipeline pressure difference, switching component position, and detection equipment operating status. All records are simultaneously recorded with mobile mission identifier, equipment identifier, timestamp, and quality mark.

[0063] The multi-particle-size calibration pulse assembly injects a stable, multi-disperse aerosol pulse lasting 2 seconds into a controllable calibration chamber upstream of the actual sampling inlet. The calibration chamber is equipped with a stirring fan and upstream / downstream concentration uniformity check ports. A low-hysteresis bypass directly samples from the calibration chamber, bypassing the actual sampling inlet by connecting to the reference particle size detection channel via a 0.25-meter-long conductive hose. The main sampling path passes sequentially from the same calibration chamber through the actual sampling inlet, the actual sampling pipeline, the drying section, and the particle size spectrometer. An electrically controlled three-way valve forms the switching component. The two detection channels, cross-calibrated with standard particles, simultaneously measure and convert the measured flow rate to the same volumetric flow rate. Before the start of the mobile mission, when the sampling pump flow rate changes by more than 5% of the calibration value, or when the pipeline pressure difference deviates from the calibration range for 30 consecutive seconds, the controller creates a calibration mission record, recording the triggering reason, target particle size set, current operating condition version, pulse duration, target flow rates for both paths, 120-second execution limit, and the current valid calibration version. The system then sequentially executes four steps: calibration chamber mixing, two-path flow stabilization, synchronous sampling, and restoration of the main sampling path. The system operates to obtain two comparable integral responses under the same upstream input. From the creation of the calibration task, the controller sets the data link status to calibration isolation and binds all particle size detection records, positioning records, meteorological records, and operating condition records generated during calibration chamber mixing, path switching, two-path stabilization, pulse injection, synchronous sampling, and restoration of the main sampling path to the same calibration task identifier and marks them as calibration-only. After the main sampling path is restored, the average value of the baseline response of the main sampling path 10 seconds before the pulse is used as the restoration baseline. The restoration tolerance is defined as three times the standard deviation of the 10-second response and 1% of the peak value of the current pulse, where the latter is not less than the other. Only when the absolute value of the difference between the baseline response of all target particle sizes and the restoration baseline is not higher than the restoration tolerance for 10 consecutive 1-second detection cycles, the calibration isolation ends and a new normal observation starting point is created. If the restoration conditions are not met within the 120-second execution limit, the isolation status is maintained, the release of normal simulation results for this cycle is stopped, and a manual review alarm is generated. Records after the timeout are not automatically converted into normal navigation records.

[0064] For each particle size class, the controller reads two raw sequences, channel sensitivity coefficients, and measured flow rates from the calibration task record. It then sequentially performs dark count subtraction, three-point median denoising, channel sensitivity conversion, volumetric flow rate conversion, and baseline zeroing. The rising edges of the same calibration pulse are aligned to a unified zero time. Subsequently, the processing rule version, input reference sequence, response output sequence, converted flow rate, alignment offset, and quality check results are written into the calibration data packet. The low-hysteresis bypass output serves as the input reference sequence upstream of the actual sampling inlet, and the main sampling path output serves as the input reference sequence after passing through the actual sampling inlet and... The response output sequence after the sampling chain; if the relative difference of the integral response of the two uniformity check ports of the calibration chamber is greater than 3%, the relative difference of the two converted volumetric flow rates is greater than 2%, the peak value of any sequence is lower than 2% of the full scale of the corresponding detection channel, or three consecutive saturated sampling points appear, the current calibration will be marked as invalid; if the particle size calibration record has a previous valid record, it will be used; if it does not exist, an empty version identifier and a calibration unavailable error code will be written, the particle size will be stopped from being transmitted to S2, the acquisition will be maintained, and recalibration will be performed when the next trigger condition is reached, without generating an alternative value for the inlet sampling efficiency or the transmission core;

[0065] The inlet sampling efficiency is calculated according to the formula. Sure, Indicates particle size Dimensionless inlet sampling efficiency Indicates the particle size index. Indicates the index of the sampling point in the calibration sequence. This represents the baseline response of the main sampling path at the corresponding sampling point, expressed in micrograms per cubic meter. This represents the baseline response of the low-hysteresis bypass at the corresponding sampling point, expressed in micrograms per cubic meter. The calibration sampling interval is 0.1 seconds in this specific embodiment. The integration window covers the period from 2 seconds before the pulse arrives to 5 seconds after the main sampling path response recovers to less than 1% of the peak value. The resulting inlet sampling efficiency is limited to between 0 and 1. If it exceeds this range, the calibration is deemed invalid.

[0066] The discrete delay-stretched transfer kernel is obtained through non-negative convolution fitting. The controller divides the delay into 0.5-second units within a delay limit of 60 seconds and solves the problem. And satisfy , Indicates particle size In delay unit Dimensionless kernel weights Indicates the particle size index. Indicates the delay unit index, Indicates the index of the sampling point in the calibration sequence. This represents the baseline response of the main sampling path, expressed in micrograms per cubic meter. Indicates the input reference sequence at a relative delay The baseline response at each location is identical in units. This represents the dimensionless inlet sampling efficiency determined by the previous formula. This represents the maximum absolute value of the peak values ​​of the two baseline sequences, expressed in micrograms per cubic meter. The dimensionless smoothing coefficient representing the suppression of abrupt changes in adjacent nuclear weights is selected by the minimum leave-one-out verification error of 20 sets of effective calibration pulses; After normalization, both the data terms and the smoothing term are dimensionless. The solution is obtained using the projection acceleration gradient method, and the L2 norm of the kernel weight change is less than 1. Stop at this time;

[0067] The controller calculates the kernel mean and kernel covariance based on the kernel weights, and adopts... and , Indicates particle size The kernel mean, in seconds. Indicates particle size The one-dimensional kernel covariance, with units of square seconds. This represents the particle size level obtained by fitting the aforementioned non-negative convolution. In delay unit Dimensionless kernel weights Represents delay unit The delay amount, in seconds. and These represent the particle size index and the delay unit index, respectively. The kernel support interval is the smallest continuous delay interval corresponding to the cumulative kernel weight from 0.5% to 99.5%. The kernel weight outside the interval is reset to zero and then renormalized.

[0068] The particle size calibration record uses the particle size level and sampling chain operating condition version as a joint key. The fields include kernel weight array, delay array, kernel mean, kernel covariance, inlet sampling efficiency, sampling pump flow rate, pipeline temperature, pipeline pressure difference, equipment status, calibration time, valid flag, and check hash. In this specific embodiment, a record has a particle size of 1-2.5 micrometers, a flow rate of 16.7 liters per minute, a pipeline temperature of 28 degrees Celsius, a pressure difference of 1.2 kPa, a kernel mean of 8.5 seconds, a kernel covariance of 6.8 square seconds, and an inlet sampling efficiency of 0.84. The record is generated by the median of three valid pulse results. When the operating condition exceeds the valid range of the record, the record with the smallest Euclidean normalized distance is selected and an operating condition extrapolation flag is added.

[0069] The wind field reconstruction module first converts the satellite positioning latitude and longitude into local coordinates (East-North-Sky) with the southwest corner of the navigation area as the origin. It then uses the quaternions output by the inertial unit for rotation from the meteorological instrument coordinate system to the local coordinate system, and follows the... Get the first Environmental wind speed of a valid meteorological record This represents the environmental wind speed vector in the local coordinate system, with units of meters per second. This indicates a synchronized meteorological record index. Indicates the first The dimensionless attitude quaternion output by a synchronous inertial record. This represents the dimensionless rotation matrix determined by the attitude quaternion. This represents the relative wind speed vector measured by the weather instrument, with units of meters per second. This represents the carrier velocity vector obtained by differential positioning and five-point Savitzky-Golay filtering, with units of meters per second.

[0070] On a reconstruction grid with a horizontal dimension of 10 meters, a vertical dimension of 5 meters, and a time interval of 5 seconds, the wind field reconstruction module uses spatiotemporal Gaussian weights from each valid meteorological record for interpolation, with the weights satisfying the following conditions: , and according to Generate a reconstructed wind field. Representing records Position and time Dimensionless weights The first value obtained by the aforementioned attitude rotation and carrier velocity compensation is... A valid meteorological record of environmental wind speed vector, with units in meters per second. This indicates the center position of the grid to be reconstructed, in meters. Indicates the time to be reconstructed, with the unit being seconds. Representing records The location and the unit is meters. Representing records The timestamp is in seconds. The dimension represents a spatially relevant scale, and in this specific embodiment, it is 30 meters. The time-related scale is 15 seconds in this specific embodiment. This indicates the reconstructed wind speed, expressed in meters per second. This represents the index of valid meteorological records. For each reconstructed grid center, the wind field reconstruction module selects wind speed records after attitude rotation and carrier velocity compensation for 30 seconds before and after the center, with a horizontal distance not exceeding 30 meters. First, it performs first-order least squares detrending on the eastward, northward, and vertical components over time. Then, it removes outlier residuals by using three times the absolute deviation of the median residuals for each component, and uses the unbiased sample variance of the remaining residuals to obtain the turbulent velocity variance in the three directions. A window is only valid when there are at least 30 valid records and the valid time span is at least 20 seconds. The variance of each window is then... The spatiotemporal Gaussian weights are interpolated to the 5-second reconstruction time. If the current window is invalid, the time is backed up to the nearest valid variance no earlier than 5 minutes ago within the same navigation mission and a low confidence flag is written. If the valid variance does not exist, the turbulence variance status of the grid is written as unavailable and the grid is prohibited from entering the S2 observation zone and S3 random diffusion calculation. The turbulence velocity variance in the three directions, the start and end times of the window, the number of valid samples, the number of outliers removed, the quality flag, and the wind field version consistent with the mean wind field are written together into the mission data packet for S2, S3, and S5 to read according to the same version.

[0071] When the sum of the Gaussian weights of a certain reconstructed grid is less than 0.05, the wind field reconstruction module uses the most recent valid meteorological record and marks the grid as a low-confidence grid. When the satellite positioning quality flag is invalid, the position is extrapolated by inertial velocity within a maximum of 3 seconds and an extrapolation flag is added. If the positioning is not restored after more than 3 seconds, the generation of new mobile observation records is stopped. Finally, the quality-controlled multi-source time series, particle size calibration records, inlet sampling efficiency, discrete delay-broadening transfer kernel, and reconstructed wind field with confidence flags are written into the mission data package for S2 to construct the particle size spatiotemporal observation band.

[0072] In this specific embodiment, S2 includes:

[0073] The spatiotemporal observation band construction module reads the task data packet from S1, uses the timestamp of the particle size detection record as the detection time, and organizes the observation items according to particle size, detection record, and delay unit. Only delay units with valid particle size calibration records, inlet sampling efficiency greater than zero, kernel weight within the kernel support interval, and data link status of normal observation are included in subsequent calculations. The particle size detection value retains the original value, quality flag, and corresponding sampling link operating condition version. Records bound to calibration task identifiers or in calibration isolation state are not allowed to enter the particle size spatiotemporal observation band, S3 source-receptor matrix, S4 source strength inversion, S5 mobile background sample, detection peak determination, prediction detection response, or residual discrimination window in this step. After calibration isolation is completed, the module starts using the currently valid kernel version after receiving from the new normal observation starting point created in S1. Any discrimination window that crosses the calibration isolation zone at any start and end time is discarded and the valid sample count is cleared to zero. Only after accumulating no less than 30 normal observation samples is it allowed to generate drift identifiers or wind field corrections.

[0074] Each historical actual sampling time is according to Sure, Indicates particle size The Each detection record corresponds to a delay unit The actual historical sampling time, in seconds. This indicates the detection time of the particle size measurement record, and the unit is seconds. This represents the delay unit defined in S1, with the unit being seconds. , and These represent the particle size index, detection record index, and delay unit index, respectively. Observations whose actual historical sampling time exceeds the start time of this mobile mission are discarded and the remaining kernel weights are renormalized.

[0075] Piecewise linear interpolation is performed between two valid positioning records before and after the actual historical sampling time. The east, north and elevation coordinates use the same time scale. The heading is obtained by projecting the inertial quaternion onto the horizontal plane after performing spherical linear interpolation. If the time interval between adjacent positioning records exceeds 3 seconds or the time interval between inertial records exceeds 0.5 seconds, the observation is written into the positioning extrapolation or attitude low confidence flag, and the uncertainty of subsequent observations is increased to twice the original value.

[0076] The local wind field is read from the reconstructed wind field generated by S1. The spatiotemporal observation zone construction module first locates eight spatial grids surrounding the historical travel position, and then performs spatial trilinear and temporal linear interpolation on two adjacent reconstruction times to obtain the local wind speed components in the east, north and vertical directions, as well as the turbulent velocity variance. If the historical travel position is located at the outer boundary of the reconstructed wind field, the value of the nearest boundary grid is used and the boundary extrapolation flag is written. Observations with an extrapolation distance exceeding one horizontal grid width are not included in the source strength inversion.

[0077] The relationship between particle size detection value and environmental concentration is adopted Characterization, Indicates particle size The Each test value is expressed in micrograms per cubic meter. This represents the corresponding inlet sampling efficiency and is dimensionless. This represents the dimensionless kernel weight corresponding to the discrete delay unit. Indicates historical position and historical moments The environmental concentration is expressed in micrograms per cubic meter. This indicates the historical navigation position obtained through interpolation, in meters. This represents the detection noise and unmodeled error, expressed in micrograms per cubic meter; other indices follow the aforementioned definition.

[0078] Each observation item is recorded with the actual historical sampling time, eastward coordinates, northward coordinates, elevation, heading, local eastward wind speed, local northward wind speed, local vertical wind speed, turbulent velocity variance, particle size identifier, kernel weight, inlet sampling efficiency, detection record identifier, sampling chain operating condition version, and quality mark. All valid observation items of the same detection record form an ordered band structure through the detection record identifier, with the sorting keys being detection time, particle size, and delay amount in sequence.

[0079] Before calculating the observation uncertainty, the spatiotemporal observation band construction module performs non-negative Tikhonov deconvolution on the quality-compliant particle size detection sequence up to the current detection time using the S1 current effective inlet sampling efficiency and discrete delay-stretching transfer kernel, specifically solving for... ; Indicates the particle size range to be determined. Environmental concentration time series, with units in micrograms per cubic meter. Indicates the current effective kernel weight A dimensionless lower triangular Toeplitz convolution matrix constructed in the order of sampling time. This represents a sequence of corresponding quality compliance test values, with the unit being micrograms per cubic meter. This represents the current effective inlet sampling efficiency and is dimensionless. This represents the dimensionless matrix used to perform second-order differencing on the time series. Indicates particle size The dimensionless deconvolution regularization coefficient, This represents a particle size index; every fifth valid sample in time sequence is fixed as a validation sample, and the remaining samples are fitted samples. to The candidate with the smallest mean square error for verification was selected from 13 logarithmic candidates with a base-10 interval. When there are ties for the minimum, the larger candidate value is fixed, and the solution is recalculated using all valid samples; the non-negative projection gradient method is used, and the relative change of the objective function is no greater than [value missing]. Furthermore, the infinite norm of the projected gradient is no greater than the initial value. The process stops when the time limit is reached, with a maximum of 5000 iterations. If no candidate value converges, the deconvolution for that particle size is marked as failed. This yields environmental concentration prior samples with timestamps and navigation locations, and low-resolution prior fields are formed by inverse distance weighted interpolation using a 40m x 40m x 10m spatial grid and a 10-second time interval. The deconvolution target, candidate set, validation sample division, stopping threshold, detection data version, calibration version, and prior field version are all written into the prior field record. This prior field is generated within S2 and does not read the S5 concentration field that has not yet been generated in the current cycle. When there are at least 30 valid prior samples for a certain particle size, the prior field is centrally differencing. If there are insufficient samples, deconvolution fails, or convergence is not achieved, the previous valid prior field is read. If neither exists, the gradient propagation term is set to zero, and the square of the 95th percentile of the absolute value of the first-order difference of the qualified detection sequence is included in the additional concentration variance. If the quantile is also unavailable, the square of 5% of the detection channel range is used as a conservative additional variance.

[0080] Observation uncertainty according to estimate, Indicates particle size Test records and delay unit The concentration variance corresponding to the observed item, with units of square micrograms per square cubic meter. , and These represent the particle size index, the detection record index, and the delay unit index, respectively. This indicates the detection variance obtained from repeated measurements of zero gas and standard aerosol in the particle size detection channel, with the same units. The four-dimensional spatiotemporal gradient column vector of the aforementioned low-resolution prior field has the first three components being the eastward, northward, and vertical concentration gradients, respectively, with units of micrograms per fourth square meter. The fourth component is the temporal concentration gradient, with units of micrograms per cubic meter per second. Each component is estimated using the width of two adjacent prior space grids and the central difference between two adjacent prior time points. The delay Jacobian vector is a 4x1 matrix. The first three elements are the derivatives of the historical position with respect to the delay, in meters per second, and are the negatives of the eastward, northward, and vertical velocity components of the mobile carrier at the corresponding historical sampling time. The fourth element is the derivative of the historical time with respect to the delay, negative one, and dimensionless. The mobile carrier velocity is obtained from the qualified positioning trajectory of S1. When there are valid positioning records within 2 seconds before and after the corresponding historical sampling time, the center difference is used by dividing the difference between the two record positions by the time difference. When there are two consecutive valid positioning records on only one side, the same-side difference is used. The resulting velocity is then filtered using the positioning velocity filtering parameters of S1. If there are fewer than two valid positioning records within the aforementioned 4-second range, the observation item is deleted and the reason for the unobtainable carrier velocity is written. The local wind speed is not used to replace the carrier velocity; the local wind speed is only used for wind field interpolation and particle transport. Therefore... The unit is micrograms per cubic meter per second; This represents the kernel covariance obtained from S1, expressed in square seconds. This represents the variance of wind field quality and prior missing concentrations. Records with high confidence and valid prior information are represented as zero. For low confidence grids, the variance is the square of 25% of the absolute value of the median concentration in the neighborhood of the observation. For extrapolated records, the variance is the square of 50% of the absolute value of the median. When both low confidence and extrapolation are present, the extrapolated value is used. For prior missing information, the aforementioned conservative additional variance is applied. The neighborhood median is determined from valid samples taken within a 30-second radius of the observation, excluding those taken 30 seconds before and after the observation. (Transpose symbol) This indicates matrix transpose, where all three terms retain the units of concentration variance.

[0081] The spatiotemporal observation band construction module sums the squared variances of all observations in the same detection record according to their kernel weights to obtain the observation uncertainty of that detection record. Before normalization, it sums the variances of observations with kernel weights not less than 0.005 and effective spatiotemporal mapping to obtain the original retained kernel quality, which is then written into the quality record. If the original retained kernel quality is lower than 0.95, the effective proportion of historical positions is lower than 80%, or the detection record has a saturation mark, the entire band structure corresponding to that detection record is deleted, and the reason for deletion is written. Only when the original retained kernel quality reaches 0.95 are the remaining observations pruned and the retained kernel weights re-normalized. At the same time, the normalized checksum is written into the quality record, and its deviation from 1 is required to be no more than 1. ;

[0082] When a certain particle size class has no usable banded structure within a continuous nuclear support duration, the spatiotemporal observation band construction module marks that particle size class as non-reversible for this period and does not transmit the empty observation set to S3. When all particle size classes are non-reversible, the data conditions for this period are deemed infeasible, the source strength inversion for this period is terminated, and if a previous valid simulation result exists, that result is retained; otherwise, an empty version identifier, a non-reversible error code for this period, and a no-simulation-result status are written, and S3 to S5 are stopped. Acquisition continues until at least one particle size re-forms a valid banded structure. Otherwise, the spatiotemporal observation bands of the particle sizes that have passed the quality check, the detection value vector, the observation uncertainty vector, and the particle size availability flag are written into the spatiotemporal observation band file for joint reading by the inverse Lagrange calculation of S3 and the predictive detection response back calculation of S5.

[0083] In this specific embodiment, S3 includes:

[0084] The source-receptor matrix construction module uses the particle size spatiotemporal observation band and available particle size flags of S2 as inputs to the receptor end. It only reads observations and constructs matrices for particles marked as invertible in this period. For non-invertible particles, it writes particle size identifiers, unestimateable states, and cause codes. It does not create empty observation matrices or zero-response matrices. For available observations, it first sorts them according to the detection record identifier, particle size level, and delay amount verification order and verification hash, deletes duplicate observations, rejects records with missing coordinates, kernel weights, or observation uncertainties, and sets kernel weights and deviations exceeding 1. The strip structure was renormalized, and a 20m x 20m x 5m candidate source grid was established in the area 500m beyond the cover-cleaning trajectory. The horizontal range of each candidate source grid was taken from the corresponding 20m x 20m grid boundary, and the vertical range was taken from a 5m thick layer centered on the effective release height and trimmed by the surface elevation and the top of the simulation domain. The remaining space after deducting the underground portion and the portion occupied by buildings was defined as the effective release volume of the grid. The grid record includes the grid identifier, horizontal boundary, effective release height, clipped vertical boundary, effective release volume, surface type, and source base version. The volume normalization in the subsequent reverse Green's contribution and the S5 forward release are based on this uniform volume source base. For each available particle size observation, 2000 equally weighted reverse particles are released from its historical travel position. The reverse integration time is 30 minutes and the integration step size is 1 second.

[0085] The displacement of the opposite particle according to renew, Represents the opposite particle In the integration step The position vector and its unit is meters. Indicates the inverse particle index. Indicates the index of the reverse integration step. This represents the local wind speed obtained from the reconstructed wind field interpolation, with units of meters per second. This indicates the integration time, expressed in seconds. This represents the inverse integration step size, expressed in seconds. This represents a diagonal matrix of turbulent diffusion coefficients, with units of square meters per second. This represents a dimensionless random vector whose components follow a standard normal distribution. Represents an upward unit vector. Indicates particle size The gravitational settling velocity is expressed in meters per second.

[0086] The turbulent diffusion coefficient is obtained by multiplying the turbulent velocity variance saved by the S1 reconstructed wind field with the Lagrange integration time. In this specific implementation, the horizontal integration time is 30 seconds and the vertical integration time is 10 seconds. The diffusion coefficient is limited to 0.1 to 100 square meters per second. The random vector is generated by a fixed seed jointly generated by the mobile task identifier, detection record identifier, particle size level and reverse particle index to ensure that the same particle size source-receptor matrix is ​​obtained by repeated runs.

[0087] Gravity settling velocity is corrected for low Reynolds number spherical particles. , Indicates particle size The gravity settling velocity, calculated using the actual particle density, volume equivalent physical diameter, gravitational acceleration, Cunningham slip correction factor, and aerodynamic viscosity, is expressed in meters per second and is taken as a positive value. In the reverse integral, it acts in the opposite direction of gravity settling, and in the forward integral, it acts in the vertically downward direction. Indicates particle size The actual particle density used, and in this specific embodiment, is 1500 kg per cubic meter. S1 represents the geometric mean of the volumetric equivalent physical diameter of this particle size class, calculated based on optical response, refractive index, shape factor, and density calibration batch, with units in meters. It does not use the aerodynamic diameter already converted with reference density. This represents the acceleration due to gravity, taken as 9.80665 meters per second squared. Represents the equivalent physical diameter based on the volume. The dimensionless Cunningham slip correction factor was calculated. Expresses the aerodynamic viscosity converted from meteorological temperature, with the unit being Pascal-second (Pa·s). This indicates the particle size index; when material parameters are missing, the default combination of density 1500 kg / m³ and shape factor 1.2 is used, along with a low confidence flag and the applicable particle size range.

[0088] When the reverse particles approach the surface, the deposition rate is read according to the surface type. The surface deposition parameters are recorded with the surface type and grain size as the joint key. The fields include deposition rate, unit, effective humidity range, calibration batch and version number. In this specific embodiment, one record is low vegetation, grain size 1-2.5 micrometers, deposition rate 0.008 meters per second, effective relative humidity range 30%-80% and version 2026A. The record is obtained by fitting the wind tunnel deposition test of the same batch of standard particles and the boundary value of the adjacent interval is used when the surface humidity exceeds the effective range.

[0089] Retention weight according to renew, Represents the opposite particle In the integration step Dimensionless retention weights that have not yet been deposited Indicates the inverse particle index. Indicates the index of the reverse integration step. Indicates particle size In surface type The deposition rate on the surface is expressed in meters per second. Indicates the particle size index. Indicates the surface type index. This indicates the reverse integration step size used in the aforementioned reverse particle displacement formula, with the unit being seconds. The height of the near-surface grid is 5 meters in this specific embodiment. When the particle position is below the surface, it is mirrored back to the near-surface and the calculation continues according to the updated retention weight.

[0090] When a reverse particle crosses the boundary of the reconstructed wind field, it continues to integrate using the nearest boundary wind speed and multiplies its contribution by 0.5. When it crosses the candidate source region by more than one horizontal grid width or the retention weight is less than 0.001, the particle is terminated. When it encounters a building grid, it is randomly displaced by mirror reflection according to the surface normal and the advection displacement is retained. All termination reasons, final positions and final retention weights are written into the reverse particle state record.

[0091] For each receptor observation and candidate source grid, the inverse Green's contribution is first calculated. , Indicates particle size Test records Delay unit With candidate source grid The transmission contribution between them, expressed in seconds per cubic meter. , , and These represent the particle size index, detection record index, delay cell index, and candidate source grid index, respectively. This represents the inverse integration step size, expressed in seconds. This represents the total number of antiparticles released for each observation. Represents candidate source grid The volume, and the unit is cubic meters. Represents the opposite particle In the integration step Dimensionless retention weight, This represents the position vector of the inverse particle at this integration step, in meters. Indicates the inverse particle index. This represents the index of the reverse integration step from the release step to the particle's termination step. The outer layer sums over all valid integration steps. Indicates the position of the reverse particle falling into the grid. A dimensionless indicator function that takes the value 1 if the condition is met and 0 otherwise.

[0092] Particle size source-receptor matrix elements according to form, Represents candidate source grid Particle size Unit emission rate for the first The response coefficient of each detection value is given in seconds per cubic meter. It represents the inlet sampling efficiency and is dimensionless. Indicates particle size In delay unit The kernel weight is dimensionless. Indicates the delay cell index within the core support interval. This represents the particle size range calculated using the aforementioned reverse Green's contribution formula. Test records Delay unit and candidate source grid The transmission contribution between them, expressed in seconds per cubic meter. , and These represent the particle size index, detection record index, and candidate source grid index, respectively, thereby directly incorporating the weights of the discrete delay-broadening transfer kernel into the observation relationship between the candidate source grid and the particle size detection value;

[0093] The source-receptor matrix construction module calculates the relative change of non-zero matrix elements every 500 inverse particles added. When the 95th percentile of two consecutive relative changes is less than 5%, the particle release for that observation term is stopped early. If convergence is not achieved, the number of inverse particles is increased to a maximum of 8000. Finally, the non-zero matrix elements, Monte Carlo standard error, boundary extrapolation ratio, and matrix version are written into a compressed sparse line file using particle size, detection record identifier, and candidate source grid identifier as keys, for S4 to perform particle size source emission rate inversion.

[0094] In this specific embodiment, S4 includes:

[0095] The source emission rate inversion module reads the particle size source-receptor matrix of S3 and the detection value vector, observation uncertainty vector, and particle size availability flag of S2. It only arranges the detection values ​​of invertible particle sizes in this period as receptor observation vectors and arranges the emission rates to be estimated only by candidate source grids and invertible particle sizes. For non-invertible particle sizes, no optimization variables are established, and they are written into the non-estimateable state instead of the zero emission state. The unit of the source-receptor matrix coefficients is per cubic meter per second, and the unit of the emission rate is micrograms per second. The product of the two and the detection value are both micrograms per cubic meter.

[0096] The diagonal elements of the observation weight matrix are the reciprocals of the corresponding observation standard uncertainties, specifically: , Indicates particle size The The observation weight of each detection record is expressed in cubic meters per microgram. This represents the standard uncertainty of the observation obtained by summing up S2, with the unit being micrograms per cubic meter. This represents the lower limit of uncertainty in the same unit, determined by the standard deviation of repeated zero-gas measurements. and These represent the particle size index and the detection record index, respectively, thus making the weighted residual dimensionless and preventing a single minimal uncertainty from dominating the solution;

[0097] The particle size source emission rate vector is solved by get, This represents the emission rate vector for all candidate source grids and particle size levels, expressed in micrograms per second. Represents a diagonal matrix consisting of observation weights, such that the weighted residuals are dimensionless. This represents a vector of detected values ​​in micrograms per cubic meter. This represents the particle size source-acceptor matrix, with units of seconds per cubic meter. This represents the spatial sparsity penalty coefficient, with units of micrograms per second. This represents the dimensionless first-order difference matrix constructed from horizontally adjacent candidate source grids. Indicates particle size The emission rate subvector across all candidate source grids, in micrograms per second. This represents the group penalty coefficient for cross-particle size common source activation, expressed in micrograms per second. Represents candidate source grid The multi-particle-size emission rate subvector is expressed in micrograms per second. and These represent the granularity-level index and the candidate source grid index, respectively; therefore, the data item and the two penalty terms are dimensionless.

[0098] To ensure the three target items have comparable scales, the source emission rate inversion module first normalizes the source-receptor matrix columns of each particle size according to their L2 norm and restores the emission rate scale after solving. At the same time, it converts the penalty coefficients back to physical units of per second per microgram according to the corresponding column scale. The spatial sparsity penalty coefficient and the group penalty coefficient are selected from the logarithmically equidistant candidate set. The candidate ranges are 0.001 to 0.1 times and 0.002 to 0.2 times the maximum absolute value of the initial gradient of the data fitting term with respect to the physical emission rate, respectively. The gradient unit is per second per microgram. The weighted prediction error is calculated using the validation subset extracted every five mobile observation points, and the coefficient combination with a validation error not higher than the minimum value plus one standard error and fewer non-zero grids is selected.

[0099] The solution employs the alternating direction multiplier method, performed on the dimensionless variables after column normalization. It separates the updates of non-negative emission rates, spatial difference soft thresholds, and cross-size group soft thresholds. The spatial difference operator and the cross-size group identity operator are stacked row-wise to form a dimensionless splitting operator. In the Rotational original residuals Dual residuals ,in Indicates the first Round normalized nonnegative emission rate variable, This represents the normalized auxiliary variable that supports the results of spatial difference and cross-particle size grouping. Indicates the first Dimensionless penalty parameter of the wheel, and Let them represent the dimensionless original residual vector and the dual residual vector, respectively. Indicates the iteration round index. This represents the matrix transpose; in each round, the non-negative least squares subproblem with a quadratic penalty is solved first, then element-wise soft thresholding is performed on the differences between adjacent grids, and grouping is performed on the multi-granular radial quantities of the same grid according to their second norm. Finally, the auxiliary variables and scaling dual variables are updated; the penalty parameter is initially set to 1, when... The penalty parameter for the next round is twice the current value. The penalty parameter for the next round is half of the current value, while all other cases remain unchanged, and the penalty parameter is limited to... to When the penalty parameter changes, the scaling dual variable is multiplied by the ratio of the penalty parameter before adjustment to the penalty parameter after adjustment to maintain the continuity of the unscaled dual variable.

[0100] Dimensional adaptive stopping tolerance according to and calculate, Indicates the first Dimensionless original residual tolerance of the wheel Indicates the first dimensionless dual residual tolerance of the wheel Represents the splitting operator the number of rows, Represents the normalized emission rate variable The total number of elements, Indicates absolute tolerance and takes a fixed value. , Indicates relative tolerance and is fixed. , Indicates the first The dimensionless, unscaled dual variable corresponding to the wheel and split constraints. , , and Continuing with the previous definition; when the relative change in the objective function between two adjacent rounds is less than... , and The solution is stopped when the maximum number of rounds is 2000. If the solution fails to converge after reaching the limit, the spatial sparsity penalty coefficient and the group penalty coefficient are increased by 20% and the solution is restarted from the current non-negative solution. If the solution fails to converge again after the second round, the solution that most recently satisfies the non-negative constraint is output and written into the non-convergence quality flag. The ordinary least squares result after truncation of negative values ​​is not output.

[0101] In the solution results, grids with emission rates less than 0.5% of the maximum emission rate for the same particle size are set to zero. Adjacent, non-zero candidate source grids are merged into source regions based on four-neighbor connectivity. The source region record includes the source region identifier, grid members, emission rates for each particle size, total emission rate, number of active particle sizes, fitting residuals, penalty coefficient version, and quality flag. When the optimized non-negative particle size source emission rate vector is all zero and the weighted residual sum of squares is higher than the 95th percentile of the background survey data, no zero-emission-rate grid is marked as a source region. Write the error codes for no available source terms and model mismatch, generate empty source term results for the binding matrix version, observation band version, residual value and threshold version, and stop the current cycle S5. Retain the previous valid source term result but do not spoof it as the output of the current cycle. S3 and S4 are re-executed after 30 new mobile observations are added. When the weighted residual square sum of the all-zero vector is not higher than the quantile, write the current cycle's no active source state and zero source term result. S5 only retains the background baseline result and does not release local source forward particles.

[0102] The zeroing ratio of 0.5% is determined by controlled release data with known locations and known multi-size emission rates. During calibration, the selection criteria are that the retention rate of the real source grid is not less than 95% and the false activation rate of the blank grid is not higher than 5%. The predicted value of the validation subset is obtained by multiplying the source-receptor matrix and the candidate emission rate vector and compared with the reserved measured detection values ​​one by one. If all combinations of penalty coefficients make the validation residual more than twice the unconstrained baseline, the current observation matrix is ​​determined to be unidentifiable. If the previous valid source term result exists, its version is retained. If it does not exist, the empty source term version, the source term unidentifiable error code and the status of the no-source term result are output and the current cycle S5 is stopped. After accumulating observation data, S3 and S4 are re-executed.

[0103] The source emission rate inversion module writes the non-negative source emission rate subvectors of invertible particle sizes, the list of non-invertible particle sizes, source region records, and inversion quality records into the source term results file. The candidate source grid order and invertible particle size order in the file are consistent with the particle size source-receptor matrix of S3. Non-invertible particle sizes only save the unestimateable state and must not be filled with zero values. The matrix version, observation zone version, particle size availability flag version, and penalty coefficient version are written into the file header for S5 to decide whether to publish the complete simulation results.

[0104] In this specific embodiment, S5 includes:

[0105] The diffusion simulation module reads the particle size source emission rate vector, non-reversible particle size list, and source region record from S4. When the non-reversible particle size list is not empty, the S4 result for this cycle is marked as a partial inversion result for diagnostic purposes only. Forward particle release and full result publication for this cycle are stopped. It is prohibited to mark the partial summation of available particle sizes as total particulate matter concentration, complete particle size composition, pollution center, pollution boundary, sedimentation range, or as a threshold or closed-loop determination for requiring full particle size distribution. If a previous complete and valid simulation version exists, that version is retained and its timestamp and non-current cycle status are clearly indicated. If it does not exist, a missing particle size list, a status of no complete simulation results, and an empty version identifier are written. Acquisition continues until all particle sizes are restored to reversible state before executing S3 to S5. Forward particles are released by particle size class from candidate source grids with non-zero emission rates in each 10-second release cycle, with at least 500 particles released per particle size class per grid. The initial position of the forward particles is the same effective release volume recorded in the corresponding candidate source grid in S3. The data is generated according to a uniform volume distribution. East and north coordinates are uniformly sampled within the recorded horizontal boundaries, and vertical coordinates are uniformly sampled within the clipped vertical boundaries. Samples falling into the building's occupied portion are sampled using the same random sequence until they fall into the effective volume. The random seed is obtained by calculating the SHA-256 of an ordered string containing the mission identifier, release cycle start time, candidate source grid identifier, particle size, and forward particle index, and then taking the lower 64 bits, ensuring that the same version of input produces the same release coordinates. This uniform volume source base is divided by the reverse Green's contribution of S3. The source base used is consistent, and the source base version, effective release volume boundary, coordinate generation rules, and random seed are written into the forward particle record; the mass of a single forward particle is determined according to... Sure, Indicates particle size From candidate source grid The mass of the released single particles, measured in micrograms. This indicates the emission rate of the corresponding particle size source, expressed in micrograms per second. This indicates the release period, which is 10 seconds in this specific embodiment. This indicates the number of forward particles released at the corresponding grid and particle size level within that period. and These represent the particle size index and the candidate source grid index, respectively.

[0106] Each forward particle records the particle size, current location, release location, current mass, release time, sedimentation state, and random seed. Location updates use the same reconstructed wind field, turbulent diffusion coefficient, gravity settling velocity, and surface sedimentation parameters as in S3, but the advection term is along the local wind speed direction, and the gravity settling term is along the vertically downward direction. The forward integration step size is... The time is 1 second; when a particle escapes the simulation domain, the escape state is written, and it enters a height of 1 second. When near the ground, according to Determine the deposition probability in this step. Indicates particle size In surface type The dimensionless deposition probability on the surface This indicates the deposition rate of the version with the same reverse retention weight as S3, expressed in meters per second. This represents the forward integration step size, in seconds. This indicates the near-surface grid height as 5 meters. and These represent the particle size index and the surface type index, respectively. The controller generates uniform random numbers from 0 to 1 using a fixed random seed for the particles. When the probability is less than the specified probability, the particle is marked as deposited and the position update is stopped. When crossing the surface, the particle is first mirrored back to the near-surface layer and then a probability determination is performed. Undeposited particles can re-enter the near-surface layer and be determined again, so that the forward non-deposited probability and the reverse retention weight of S3 use the same exponential deposition operator and parameter version.

[0107] The mass of undeposited forward-facing particles is accumulated using a 20m x 20m x 5m spatial grid to obtain the concentration field for each particle size. The particle size concentration is calculated according to... calculate, Indicates time Lower particle size In spatial grid The concentration is expressed in micrograms per cubic meter. Indicates the particle size index. Indicates the forward simulation space grid index. This indicates the forward output time, measured in seconds, starting from the task's inception. This indicates that the current moment is within the grid. And undeposited grain size Forward particle assembly, Represents forward-moving particles The current mass, in micrograms. Representing spatial grid The volume is expressed in cubic meters. Indicates the forward particle index;

[0108] The mass of deposited particles was accumulated using a 20m x 20m surface grid and 60-second time intervals. Settlement flux was calculated using... , Indicates particle size In the surface grid and time period Settlement flux within the area, expressed in micrograms per square meter per second. Indicates the particle size index. Indicates the surface grid index. This represents the index of the cumulative settlement period starting from the beginning of the task. This indicates the grain size that transitioned to a deposited state within that grid and time period. Forward particle assembly, Represents forward-moving particles The current mass, in micrograms. Indicates the forward particle index. This represents the area of ​​the land surface grid, expressed in square meters. This indicates the cumulative settlement time interval, which is 60 seconds in this specific embodiment;

[0109] The total particulate matter concentration field is obtained by summing the particle size concentrations of the same spatiotemporal grid. The pollution center takes the total particulate matter concentration as the mass weight of the horizontal coordinate weighted mean. The diffusion direction calculates the concentration weighted horizontal coordinate covariance matrix within a 300-meter neighborhood around the pollution center and takes the unit eigenvector corresponding to the largest eigenvalue. If the ratio of the two eigenvalues ​​is less than 1.1, the diffusion direction is marked as uncertain. The particle size composition is written into each grid according to the proportion of each particle size concentration to the total particulate matter concentration, and all proportions are set to zero when the total concentration is zero.

[0110] The background area for mobile monitoring consists of valid detection records for 5 consecutive minutes before the start of the mission, which are more than 500 meters away from any inversion source area. For each particle size, the diffusion simulation module directly takes the 95th percentile of the qualified original detection value sample in the background area as the particle size detection background threshold, and binds the particle size, inlet sampling efficiency version, transfer kernel version, and background sample version for subsequent determination of synchronous pollution peaks. When there are fewer than 60 valid detection samples for a certain particle size, saturation occurs, or the background area does not exist, the detection background threshold and synchronous pollution peak identifier for that particle size are set as undeterminable and enter the pending determination path, without borrowing the total environmental concentration threshold.

[0111] The diffusion simulation module then performs non-negative regularized deconvolution on the background detection sequence according to the entry sampling efficiency and discrete delay-broadening propagation kernel bound to each record. This deconvolution calls the same non-negative Tikhonov objective function, second-order difference operator, 13 candidate regularization coefficients, validation partitioning every fifth sample, the rule of taking the larger coefficient in parallel, the projection gradient stopping threshold, and the upper limit of 5000 iterations defined in S2 for each particle size. It also writes the detection sequence version, kernel version, entry sampling efficiency version, regularization coefficient, and convergence status into the background deconvolution record, thereby obtaining the environmental background concentration sequence for each particle size consistent with the simulation field. If the deconvolution of any particle size fails or does not converge, that particle size is marked as zero. The conversion method does not use truncated negative values ​​or direct division by the inlet sampling efficiency as a sequence replacement; for all successfully converted particle sizes, synchronous pairing within a 1-second tolerance is performed according to the unified historical actual sampling time, and the total particulate matter environmental background concentration sample is formed by summing time by time. The pollution boundary threshold is directly taken as the 95th percentile of the synchronous total concentration sample; when there are fewer than 60 effective synchronous samples or any particle size cannot be converted, the pollution boundary is marked as undetermined. If there is a previous effective threshold, its version record is used. If there is no previous effective threshold, an empty threshold version is written and the pollution boundary is not extracted; the diffusion simulation module performs eight-neighbor connectivity on adjacent spatial grids with total concentration not lower than the effective environmental concentration threshold and extracts the outer envelope as the pollution boundary.

[0112] The background deposition flux baseline uses the aforementioned environmental background concentration sequences for each particle size, the concurrently reconstructed wind field, and the same deposition parameters as the formal simulation. The background forward model is run without setting local inversion sources. For each 60-second interval, the upwind boundary is determined based on the horizontal average wind direction. The environmental background concentration sequence is linearly interpolated to this interval, and background particles are continuously released at each 20m x 20m x 5m grid along the upwind boundary in a uniform horizontal distribution and an exponentially decaying vertical distribution from the ground upwards. The scale height of the exponential distribution is 20 meters, calculated from the normalized mean square error of the concentration within 0 to 100 meters in the vertical profile calibration of the background in the underway area. The minimum candidate scale is determined, and the upper limit of the release height is taken as the value that is no greater than either 100 meters or the top height of the simulation domain. The first vertical grid is calculated from the ground. For each 5-meter vertical grid, the definite integral of the exponential function between the upper and lower boundaries of that grid is used as the original weight. Then, the weights are normalized by the sum of the original weights of all effective vertical grids, so that the sum of the mass shares of each height grid equals 1. When the background vertical profile calibration is invalid, a conservative scale height of 10 meters is used and a low confidence flag is added. When there are fewer than two effective vertical grids within the release limit, that period is marked as unusable. The total particle mass is calculated based on the environmental background concentration and the normal wind speed at the inflow boundary. The product of the boundary grid area and the 60-second duration is used to determine the boundary, which is then distributed to each height grid according to the aforementioned normalized mass share. 500 equal-mass particles are released per particle size per boundary grid. When the normal wind speed is less than 0.1 m / s and a previous effective upwind boundary exists, that boundary is used but marked as low confidence. When the normal wind speed is less than 0.1 m / s and no effective upwind boundary has ever been formed before, that 60-second period is marked as unusable, no background particles are released, and the reasons for low wind speed and no historical boundary are recorded. An upwind boundary is only established after the first effective average wind direction is formed and the normal wind speed is not less than 0.1 m / s. The aforementioned invalid periods without historical boundaries are not... Background sedimentary flux quantile samples are included. When the wind direction reverses, the boundary is switched in the next time period. If the environmental sequence gap exceeds 10 seconds, the time period is invalid. The background model generates background sedimentary flux samples of each grain size using the same 20m x 20m surface grid and 60-second time periods. The threshold of each surface grid is taken as the 95th percentile of the corresponding sample. For grids with formal deposition flux higher than the threshold of the same grain size and grid, the outer envelope is extracted as the deposition range. If each grid has fewer than 5 valid background time periods, the wind field version is inconsistent, or the background forward model fails, the grid and the summarized deposition range are marked as undeterminable, and the air concentration threshold is not used to replace the sedimentary flux threshold.

[0113] To calculate the visible response of the detection equipment, the diffusion simulation module performs spatiotemporal linear interpolation on the forward concentration field based on the particle size spatiotemporal observation band of S2 at each historical actual sampling time and travel position, and according to... Calculate the predicted detection response. Indicates particle size The A predictive detection response, with units of micrograms per cubic meter. It represents the inlet sampling efficiency and is dimensionless. Represents delay unit The kernel weight is dimensionless. This represents the predicted environmental concentration at the corresponding historical actual sampling time and navigation location, expressed in micrograms per cubic meter. and These represent the detection record index and the delay unit index, respectively.

[0114] Particle size residual according to calculate, Indicates particle size The Each residual is expressed in micrograms per cubic meter. This indicates the corresponding particle size detection value, with the unit being micrograms per cubic meter. This represents the particle size level calculated from the current forward concentration field and the spatiotemporal observation zonation of particle size according to the aforementioned predictive detection response formula. No. A predictive detection response, with units of micrograms per cubic meter. and These represent the particle size index and the detection record index, respectively. The runtime source for the predicted detection response is the current forward concentration field, the spatiotemporal observation zone of the particle size, the inlet sampling efficiency, and the kernel weight. This back calculation does not perform candidate actions or resource budget screening. When any source version is inconsistent, the historical position exceeds the simulation domain, or the environmental predicted concentration is missing, the corresponding sample is marked as uncalculation and excluded. The residual, predicted detection response, detection value, observation zone version, source term version, and wind field version are all written into the closed-loop discrimination cache.

[0115] The discrimination window length for each particle size is taken as the duration of the corresponding discrete delay-broadened propagation kernel support time and 120 seconds, where the value is not less than the other. The window slides in 10-second increments, and the residual sign ratio is adopted. , Indicates particle size In the discrimination window The proportion of dimensionless residuals with the same sign within the range Indicates the index of the discrimination window. Indicates the number of positive residuals. Indicates the number of negative residuals. This represents the total number of valid residuals. The first quantile threshold is determined by the 95th percentile of the residuals with the same sign proportion in the uncontaminated calibration task. If there are fewer than 30 valid residuals in the window, a sampling chain drift flag is not generated and the insufficient sample status is recorded.

[0116] The sampling chain operating parameters include pump flow rate, pipeline differential pressure, pipeline temperature, and sheath gas flow rate of the detection equipment. The closed-loop calibration module first normalizes each parameter according to the mean and standard deviation of the stable operating condition samples, and then... Calculate the correlation coefficient. Indicates particle size Residuals and operating parameters The dimensionless Pearson correlation coefficient, Indicates the index of operating condition parameters. Indicates the first The normalized chemical condition parameters corresponding to each detection record This represents the mean of the normalized chemical condition parameters within the window. This represents the mean residual particle size within the window. This represents the particle size range calculated using the aforementioned residual formula. No. Each residual is expressed in micrograms per cubic meter. and These represent the particle size index and the detection record index within the window, respectively. The correlation coefficient is calculated only when the number of effective paired samples is not less than 30 and the sample variance of the residual and the corresponding operating condition parameter is greater than 1% of their respective stable operating condition variance. The one with the largest absolute value is taken as the operating condition correlation, and the corresponding parameter name and correlation coefficient symbol are retained. The second quantile threshold is taken as the 95th quantile of the operating condition correlation in the stable operating condition calibration samples. Uncalculable parameters do not participate in the maximum value selection. When all parameters are uncalculable, the operating condition correlation is written as unavailable, no sampling chain drift flag is generated, and the reason for zero variance or insufficient samples is recorded.

[0117] Each particle size detection peak is established when the detected value is higher than the corresponding background concentration threshold and is a local maximum within two consecutive detection points. The equivalent actual sampling time is calculated according to... calculate, Indicates particle size The The equivalent actual sampling time for each detection peak, in seconds. The kernel weight is dimensionless. This indicates the corresponding historical actual sampling time, with the unit being seconds. Indicates the particle size index. This indicates the index of the detection peak records. This indicates the index of the delay unit within the current effective core support interval and sums up all delay units in that interval; when at least two particle size levels have detection peaks higher than their respective background thresholds and the difference between the equivalent actual sampling times is no greater than 1 second in the positioning data sampling period, the multi-particle size synchronous contamination peak flag is set to yes;

[0118] When the proportion of the same sign of the residual of the target particle size reaches the first quantile threshold, the correlation of the operating condition reaches the second quantile threshold, and the synchronous pollution peak of multiple particle sizes is marked as no, the closed-loop calibration module generates a sampling chain drift identifier for the target particle size. The drift record includes the particle size level, the start and end time of the discrimination window, the proportion of the same sign, the correlation of the operating condition, the dominant operating condition parameters, the three quantile threshold versions, and the synchronous pollution peak identifier. If any of the three conditions is not met, the drift identifier is not generated.

[0119] For each particle size residual sample within the discrimination window, the closed-loop calibration module uses the kernel weights of its particle size spatiotemporal observations to weight the local eastward and northward wind speed components at the actual historical sampling times, forming equivalent eastward and northward wind speed sequences. For stable wind field survey samples, historical discrimination windows without sampling chain drift markers, multi-particle-size synchronous pollution peak markers, consistent wind field versions, and where both residuals and equivalent wind speeds are calculable are selected. These windows form variance benchmark sets for each particle size residual sample and variance benchmark sets for eastward and northward equivalent wind speed samples. The residual variance threshold is the 1st quantile of the variance benchmark set for the same particle size residual sample, with units of square micrograms per square cubic meter. The directional wind speed variance threshold is the 1st quantile of the variance benchmark set for the corresponding directional equivalent wind speed sample, with units of square meters per square second. Only one sample is used for each direction. The Pearson correlation coefficient is calculated when the number of effective paired samples is no less than 30, the sample variance of the residual sequence is strictly greater than the residual variance threshold of the same particle size, and the sample variance of the wind speed sequence in that direction is strictly greater than the wind speed variance threshold of the corresponding direction. If these conditions are not met, the direction is marked as unavailable in terms of correlation. If the effective historical window of any benchmark set is less than 30, the corresponding 1% quantile threshold is not generated, the affected direction is marked as unavailable in terms of correlation, and the reason for insufficient stable variance benchmark is written. The correlation coefficient with an absolute value not less than that of the other one is taken from the available directions as the residual-wind field correlation of that particle size and the corresponding direction is recorded. When the two absolute values ​​are equal, the eastward wind speed direction is fixed. When both directions are unavailable, the particle size is written as unavailable in terms of wind field correlation. The third quantile threshold is taken as the 95th quantile of the residual-wind field correlation in the stable wind field sailing sample.

[0120] When multiple particle size generation sampling chain drift markers exist within a single discrimination window, the controller first sorts them from highest to lowest based on the absolute value of the correlation between operating conditions. If the values ​​are the same, they are sorted from smallest to largest based on the particle size level number, forming a calibration queue. A mutual exclusion lock is set between the shared switching component and the multi-particle size calibration pulse component, allowing only the head of the queue to hold the lock and execute at any given time. The controller generates a calibration action record containing the action type, target particle size level, queue number, execution time, current core version, candidate core version, and timeout limit. It then drives the switching component to connect to a low-hysteresis bypass and applies a single pulse calibration. Simultaneously, it immediately enters the calibration isolation state defined in S1, generating candidate discrete delay-broadening transfer kernels and candidate kernels for the target particle size according to the fitting process in S1. Covariance and other particle size records remain unchanged. Candidate records use incremental version numbers and do not overwrite the current valid records before the 120-second timeout acceptance. Isolation ends, the mutex lock is released, and a new normal observation starting point is created only after the candidate reception or rollback is completed and the main sampling path recovery conditions specified in S1 are met. All records and discrimination windows during the isolation period and those crossing the isolation area are excluded according to the S2 rule. If any recovery or acceptance action times out, the current valid kernel version and isolation status are maintained and a manual review alarm is generated. After the current target completes the aforementioned reception, rollback, and recovery, the controller recalculates the drift conditions of the remaining targets with the new current version and the normal observation samples re-accumulated after isolation. Those that no longer meet the conditions are removed from the queue.

[0121] If no sampling chain drift identifier is generated and the multi-particle-size synchronous pollution peak identifier is "yes", the closed-loop calibration module keeps all particle size calibration records unchanged, calls S3 and S4 to reconstruct the source-receptor matrix and particle size source emission rate vector of the corresponding window with the current particle size spatiotemporal observation band, and continues to execute S5 to generate the concentration field, sedimentation results, predicted detection response and residuals of the same version chain; after the matrix, source terms and forward results are all successfully completed, they are submitted atomically with a unified transaction identifier. Before completion, the old source terms and old forward results are kept as the current valid versions and the window is marked as updating. Mixing the old and new versions is prohibited.

[0122] If no sampling chain drift identifier is generated, no multi-particle-size synchronous pollution peak identifier is generated, and the residual-wind field correlation of at least two particle sizes reaches the third quantile threshold, keep the particle size calibration record and the current effective source emission rate vector unchanged, apply small perturbations of ±0.05 m / s to the eastward and northward components of the reconstructed wind field respectively, and rerun the prediction detection response on the same source term version, the same discrimination window, and the same candidate effective record set; first press Calculate the wind speed sensitivity for each record. Indicates particle size The Each detection record is in the direction The wind speed sensitivity is measured in micrograms per second per cubic meter. and They represent directions respectively. Particle size distribution after applying positive and negative perturbations No. A predictive detection response, with units of micrograms per cubic meter. This indicates the magnitude of the wind speed disturbance on one side, which is 0.05 meters per second in this specific embodiment. , and These represent the particle size index, the detection record index, and the eastward or northward wind speed direction index, respectively. For the current particle size and direction, the set of record indexes that are all calculable and have the same version, including positive disturbance, negative disturbance, baseline prediction response, residual, and observation variance, is denoted as... Records lacking any input are removed from the set, and dimensionless weights are obtained by normalizing the inverse of the observed variance of each record within the set. And the sum of the weights is 1; when the set contains no fewer than 30 records, according to , and The particle size wind speed sensitivity, the average residual of the particle size window, and the reciprocal of the average observation variance of the window are obtained respectively in this direction; Indicates particle size In direction The wind speed sensitivity is measured in micrograms per second per cubic meter. Indicates the direction of the current situation Granularity obtained using the same set of records and the same weights The window average residual, expressed in micrograms per cubic meter. This indicates the aforementioned particle size residual, with the unit being micrograms per cubic meter. Indicates the direction of the current situation The reciprocal of the window mean observation variance obtained from the same set of records, in units of sixth-power meters per square microgram. Indicates the particle size distribution obtained by S2 No. The variance of each detection record is expressed in square micrograms per square cubic meter. , , , and All are defined directly according to this section; if there are fewer than 30 records in the set, the particle size in that direction is marked as unusable for sensitivity.

[0123] For each direction, only the set The available particle size is defined as a particle size with at least 30 valid records and a residual-wind field correlation reaching the third quantile threshold. The closed-loop calibration module is configured according to... Synthetic unique wind field correction amount Indicates direction The wind speed correction is given in meters per second. This indicates that the preceding part refers to the current direction. Grain size obtained from the same set of records The reciprocal of the window mean observation variance, expressed in sixth-power meters per square microgram. This indicates the particle size-level wind speed sensitivity obtained by weighting the sensitivity of each record, with the unit being micrograms per second per fourth cubic meter. This represents the average residual of the particle size window obtained from the previous segment using the same set of records and the same weights, with the unit being micrograms per cubic meter. Indicates the available particle size index in the current direction. This indicates the wind speed direction index; therefore, the numerator unit is seconds per meter, the denominator unit is square seconds per square meter, and the ratio is meters per second; the nonlinear error of positive and negative disturbances is calculated according to... calculate, Indicates direction Dimensionless nonlinear error, Represents the set of all available particle sizes in the current direction. The set obtained by merging according to particle size-detection record index This means to combine the sets The inverse of the observed variance used is again in the set Internally normalized dimensionless weights that sum to 1. and They represent directions respectively. Particle size distribution after applying perturbations of positive 0.05 m / s and negative 0.05 m / s No. A predictive detection response, with units of micrograms per cubic meter. This represents the baseline predicted response with the same particle size and detection record when no perturbation is applied, and the units are the same. This represents the maximum digital resolution of the particle size detection channels contained in the set, expressed in micrograms per cubic meter. The expression represents the resolution of a number in squares, with the unit being square micrograms per square cubic meter. Both the numerator and denominator of the candidate expression are in square micrograms per square cubic meter. , and These represent the particle size index, the detection record index, and the eastward or northward wind speed direction index, respectively; the current direction can use fewer than two particle sizes, or a set. Empty, any disturbance branch calculation failed, any absolute value of sensitivity is lower than the 5th percentile lower limit obtained from the stable wind field disturbance test, Alternatively, if the denominator of the wind field correction formula is less than the 1st quantile of the corresponding stable wind field sample, the direction is not updated and the process enters the pending determination stage; otherwise, the direction is... After limiting the amplitude to ±0.2 m / s, the reconstructed wind field is added, bound to the new wind field version, and the S2 observation zone is reconstructed using the new wind field. The current effective source emission rate is then used to perform S5 forward diffusion and predicted response backcalculation, without proceeding to S4 estimation. If only one direction meets all update conditions, only that direction is updated; the other direction remains unchanged and is recorded as pending judgment. If neither direction meets the conditions, the reconstructed wind field remains unchanged, and the entire window is set to pending judgment. If S3 or S4 is run only for diagnostic purposes, its temporary results must not be submitted or replace the current source item, and the sensitivity and set of records are recorded. and Sample identifiers, aggregation weights, , The orientation status, fresh air field version, and fixed source item version are bound and written into the closed-loop record;

[0124] When no sampling chain drift identifier is generated, the multi-particle size synchronous pollution peak identifier is not set, and there are fewer than two particle sizes that reach the third quantile threshold, the closed-loop calibration module keeps the discrete delay-broadening transfer kernel, particle size source emission rate vector, and reconstructed wind field unchanged, writes the window status as pending judgment, and enters the next discrimination window. When three consecutive windows are pending judgment, only a manual review prompt is added, without changing the current simulation results.

[0125] After generating the target particle size candidate calibration record, the closed-loop calibration module writes the current kernel version, candidate kernel version, current observation band version, and effective detection record set into the acceptance status. It defines receiving the candidate version and reverting to the current version as mutually exclusive actions. Within the same 10-second closed-loop cycle, it reconstructs the spatiotemporal observation band of the target particle size using both the pre-update and candidate discrete delay-broadening transfer kernels. It then sequentially executes the inverse Lagrange calculation in S3, the source emission rate inversion in S4, and the forward particle diffusion in this step, calculating within the same discrimination window and the same effective detection record set. , , and , This represents the sum of squared residuals before the update, expressed in square micrograms per square cubic meter. This indicates the sum of squared residuals after the update, with the same units. This represents the residual calculated using the aforementioned particle size residual formula from the pre-update kernel version, pre-update observation zone version, pre-update source term version, and pre-update wind field version, with units in micrograms per cubic meter. This represents the residual calculated from the candidate kernel version, candidate observation band version, candidate source term version, and candidate forward result using the aforementioned particle size residual formula, with units in micrograms per cubic meter. This represents the sum of the diagonal elements of the kernel covariance before the update, expressed in square seconds. This represents the sum of the diagonal elements of the updated kernel covariance, with all elements having the same unit. This represents the kernel covariance calculated from the unupdated discrete delay-stretched transfer kernel using the aforementioned kernel covariance formula, with units of square seconds. This represents the candidate kernel covariance calculated by the candidate discrete delay-stretched transfer kernel according to the aforementioned kernel covariance formula, with units of square seconds. Indicates the target granularity level index. This represents the index of the detection record within the same set of valid detection records. The summation in all four equations iterates through this same set. It represents matrix trace operations and is not applicable to physical units as an operator; the acceptance objective is to reduce the fitting residuals without increasing the kernel delay uncertainty. Candidates are only accepted when the sum of squared residuals after the update is strictly less than the value before the update and the sum of the diagonal elements of the candidate kernel covariance is not greater than the value before the update. The candidate kernels, kernel covariances, candidate observation bands, candidate source-receptor matrices, candidate source emission rate vectors, and candidate forward results are changed from temporary to current valid state in one go with a unified transaction identifier. If any object fails to be committed, the entire transaction is rolled back and all versions before the update are retained. It is forbidden to commit only candidate kernels.

[0126] If both acceptance criteria are not met simultaneously, the closed-loop calibration module generates a rollback action record containing the target particle size level, rollback kernel version, restored observation zone version, restored source term version, restored wind field version, and execution time. It rolls back to the discrete delay-broadened transfer kernel and kernel covariance before the update, and restores the particle size spatiotemporal observation zone, source-acceptor matrix, particle size source emission rate vector, and forward diffusion results to their pre-update versions. If 120 seconds have elapsed since the candidate record was generated, or if any recalculation stage from S2 to S5 fails, candidate acceptance is immediately terminated, the candidate record is marked as timed out and invalid, or the calculation failed, and no part of the recalculation results are used. The system maintains the current valid core, observation zone, source-receptor matrix, source terms, wind field, and forward results unchanged, and records the failure stage, elapsed time, and recovery version. Each candidate can be automatically retried at most once, and a manual verification prompt is added if it still fails. The particle size calibration record saves the values ​​before and after the update, the sum of squares of the two sets of residuals, the sum of the diagonal elements of the two sets of core covariances, received, rolled back, timeout invalid or calculation failure results, execution time, and associated mobile task identifiers. Finally, it outputs the concentration fields of each particle size, the total particulate matter concentration field, the pollution center, the diffusion direction, the particle size composition, the sedimentation flux, the sedimentation range, the pollution boundary, and the closed-loop calibration record.

[0127] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.

[0128] This invention integrates sampling chain error, travel trajectory, reconstructed wind field, and particle size transport into the particle size observation relationship, so that the detection response used for source strength inversion is consistent with the predicted detection response of forward diffusion back-calculation, thereby forming a continuous technology chain from calibration, observation reconstruction, inversion to diffusion simulation.

[0129] The closed-loop calibration of this invention selects the local update object of the transfer nucleus, source term or reconstructed wind field based on the discrimination condition, and performs update acceptance and rollback record with residuals and nucleus covariance, thereby maintaining the stability of particle size diffusion results while retaining the true pollution peak.

Claims

1. A method for simulating particulate matter particle size pollution diffusion based on mobile reconstruction, characterized in that, include: S1. Acquire particle size detection data, positioning data, inertial data, meteorological data and sampling chain operating data of the mobile carrier, determine the inlet sampling efficiency and the discrete delay-broadening transfer kernel formed by the combination of sampling inlet, sampling pipeline and detection equipment according to particle size, and generate the reconstructed wind field of the mobile area based on meteorological data, positioning data and inertial data. S2. Based on the inlet sampling efficiency and the discrete delay-broadening transfer kernel, the detection value of each particle size is characterized as the convolution result of the environmental concentration at multiple historical actual sampling times weighted by the kernel weight. Each kernel delay unit is mapped to the local wind field in the corresponding historical actual sampling time, the navigation position, the heading, and the reconstructed wind field to obtain the particle size spatiotemporal observation band. S3. In the reconstructed wind field, reverse Lagrange calculations including advection, turbulent random diffusion, gravity settling and surface deposition are performed for each particle size. The kernel weight of the spatiotemporal observation zone of particle size is incorporated into the transmission contribution from the candidate source grid to the particle size detection value to obtain the particle size source-receptor matrix. S4. Based on the particle size source-acceptor matrix and the detection values ​​of each particle size, estimate the particle size source emission rate vector of each candidate source grid under non-negativity constraints, spatial sparsity constraints and cross-particle size common source activation constraints. S5. Perform forward particle diffusion including advection, turbulent random diffusion, gravity settling, and surface deposition based on the particle size source emission rate vector to obtain the particle size concentration field and particulate matter settling simulation results.

2. The particulate matter particle size pollution diffusion simulation method based on mobile reconstruction according to claim 1, characterized in that, S1 includes: A low-hysteresis bypass, a main sampling path, a switching component for switching between the low-hysteresis bypass and the main sampling path, and a multi-particle size calibration pulse component are set at the determined location of the discrete delay-broadening transfer core. The same calibration pulse is applied to the low-hysteresis bypass and the main sampling path by the multi-particle size calibration pulse component. The low-hysteresis bypass output sequence is used as the input reference and the main sampling path output sequence is used as the response output. The kernel mean and kernel covariance of the discrete delay-broadening transfer kernel are identified according to the particle size. The inlet sampling efficiency is determined based on the ratio of the integral response of the main sampling path output sequence to the integral response of the low-hysteresis bypass output sequence, and the kernel mean, kernel covariance, inlet sampling efficiency, and corresponding sampling chain operating conditions are written into the particle size calibration record.

3. The particulate matter particle size pollution diffusion simulation method based on mobile reconstruction according to claim 2, characterized in that, S2 include: The detection time is determined by using the timestamp of the particle size detection data as the detection time, and the corresponding historical actual sampling time is determined according to the delay amount of each delay unit in the kernel mean value; Based on the timestamps of the positioning data, inertial data, and meteorological data, time interpolation is performed on the data before and after the historical actual sampling time to obtain the corresponding navigation position, heading, and local wind field. The historical actual sampling time, navigation position, heading, local wind field, particle size identifier, and kernel weight are combined to form a particle size spatiotemporal observation band, and the observation uncertainty of each observation item in the particle size spatiotemporal observation band is determined by the kernel covariance.

4. The particulate matter particle size pollution diffusion simulation method based on mobile reconstruction according to claim 3, characterized in that, S3 includes: releasing reverse particles of corresponding particle sizes from each travel position in the particle size spatiotemporal observation band; updating the horizontal position of the reverse particles in the reconstructed wind field according to the opposite direction of the local wind speed; generating random displacements based on the turbulence statistics of the reconstructed wind field; determining the gravity settling velocity based on particle size, particle density, and air properties; updating the vertical position of the reverse particles in the opposite direction of gravity settling according to the gravity settling velocity; determining the retention weight of the reverse particles based on the deposition velocity corresponding to the surface type; and accumulating the retention weight of the reverse particles entering each candidate source grid with the kernel weight of the particle size spatiotemporal observation band to form the matrix elements of the particle size source-receptor matrix.

5. The particulate matter particle size pollution diffusion simulation method based on mobile reconstruction according to claim 4, characterized in that, S4 includes: A data fitting term is constructed using the weighted residual between the detection values ​​of each particle size and the product of the particle size source-receptor matrix and the particle size source emission rate vector. The weight of each observation term is determined by the inverse value of the corresponding observation uncertainty. A spatial sparsity penalty is set for the particle size source emission rate in adjacent candidate source grids for the same particle size, and a common source group penalty is set for the emission rates of multiple particle sizes in the same candidate source grid. The sum of the data fitting term, spatial sparsity penalty, and common source group penalty is solved under the condition that the particle size source emission rate is non-negative to obtain the particle size source emission rate vector.

6. The particulate matter particle size pollution diffusion simulation method based on mobile reconstruction according to claim 5, characterized in that, S5 includes: Forward particles are released from the corresponding candidate source grid according to the particle size emission rate vector, and the particle size, position, mass, and deposition state of each forward particle are recorded. The mass of undeposited forward particles is accumulated according to the spatial grid to obtain the particle size concentration field. The total particulate matter concentration field is obtained by summing the particle size concentration fields. The pollution center is determined by the mass-weighted position of the total particulate matter concentration field. The diffusion direction is determined by the main axis direction of the concentration distribution in the neighborhood of the pollution center. The particle size composition is determined by the proportion of each particle size concentration to the total particulate matter concentration. The mass of deposited forward particles is accumulated according to the surface grid and time interval to obtain the deposition flux and deposition range. The pollution boundary is determined by the spatial grid outer envelope when the total particulate matter concentration reaches the boundary threshold. The boundary threshold is the preset quantile of the concentration sample in the background area of ​​the mobile survey.

7. The particulate matter particle size pollution diffusion simulation method based on mobile reconstruction according to claim 6, characterized in that, Step S5 further includes: based on the particle size spatiotemporal observation band, reading the environmental prediction concentration at each historical actual sampling time and navigation position in the particle size spatiotemporal observation band from the particle size concentration field obtained by forward particle diffusion, and weighting the environmental prediction concentration with the corresponding kernel weight and correcting it with the inlet sampling efficiency to obtain the prediction detection response; calculating the particle size residual based on the detection value of each particle size and the corresponding prediction detection response; calculating the residual sign ratio within a discrimination window with a length not less than the kernel support time of the corresponding discrete delay-broadening transfer kernel; calculating the correlation coefficient between the residual of the target particle size and the operating parameters of each sampling chain one by one; and selecting one with a value not less than other absolute values ​​from the absolute values ​​of each correlation coefficient as the operating condition correlation, while determining the multi-particle-size synchronous pollution peak identifier. When the proportion of residuals with the same sign of the target particle size reaches the first quantile threshold determined by the uncontaminated calibration residual sample, the working condition correlation reaches the second quantile threshold determined by the working condition correlation of the stable working condition calibration sample, and the multi-particle-size synchronous pollution peak identifier is not specified, a sampling chain drift identifier for the target particle size is generated. The background concentration quantile thresholds for each particle size are determined by the preset quantiles of the corresponding particle size detection value samples in the mobile background area. The equivalent actual sampling time of each particle size detection peak is determined by weighting the relevant historical actual sampling time by the kernel weight of the corresponding particle size spatiotemporal observation zone. The multi-particle-size synchronous pollution peak identifier is determined based on the condition that the detection peaks of at least two particle sizes are higher than their respective background concentration quantile thresholds and the difference between the equivalent actual sampling times is not greater than the sampling period of the positioning data.

8. The particulate matter particle size pollution diffusion simulation method based on mobile reconstruction according to claim 7, characterized in that, The processing of the sampling chain drift flag includes: for each particle size residual sample within the discrimination window, the local eastward and northward wind speed components at each historical actual sampling time are weighted by the kernel weight of the corresponding particle size spatiotemporal observation zone to obtain the equivalent eastward and equivalent northward wind speed components; the correlation coefficient between each particle size residual and the corresponding equivalent eastward and equivalent northward wind speed components is calculated, and the one with a value not less than the other absolute value is selected as the residual-wind field correlation of the corresponding particle size; when the target particle size generates a sampling chain drift flag, a low-hysteresis bypass is connected through the switching component and a calibration pulse is applied by the multi-particle size calibration pulse component; the discrete delay-broadening transfer kernel and kernel covariance of the target particle size are updated according to the new input reference and response output, while the others are maintained. The discrete delay-broadening transfer kernel for each particle size remains unchanged. When no sampling chain drift indicator is generated and the multi-particle-size synchronous pollution peak indicator is positive, the discrete delay-broadening transfer kernel for each particle size remains unchanged, and the particle size source emission rate vector is re-estimated. When no sampling chain drift indicator is generated, the multi-particle-size synchronous pollution peak indicator is negative, and the residual-wind field correlation of at least two particle sizes reaches the third quantile threshold determined by stable wind field samples, the discrete delay-broadening transfer kernel for each particle size remains unchanged, and the wind field is reconstructed based on the residual correction. When no sampling chain drift indicator is generated, the multi-particle-size synchronous pollution peak indicator is negative, and the particle size corresponding to the residual-wind field correlation that reaches the third quantile threshold is less than two, the discrete delay-broadening transfer kernel for each particle size, the particle size source emission rate vector, and the reconstructed wind field remain unchanged, the pending determination state is recorded, and the process proceeds to the next discrimination window.

9. The particulate matter particle size pollution diffusion simulation method based on mobile reconstruction according to claim 8, characterized in that, After updating the discrete delay-broadened transfer kernel and kernel covariance of the target particle size, the spatiotemporal observation band of the target particle size is reconstructed using the discrete delay-broadened transfer kernel before and after the update, and inversion and forward particle diffusion are performed. When the sum of squares of the residuals of the updated discrimination window is less than the sum of squares of the residuals of the discrimination window before the update, and the sum of the diagonal elements of the updated kernel covariance is not greater than the sum of the diagonal elements of the kernel covariance before the update, the updated discrete delay-broadened transfer kernel and kernel covariance are received; otherwise, the process reverts to the discrete delay-broadened transfer kernel and kernel covariance before the update, and the values ​​before and after the update, as well as the received or reverted results, are written into the particle size calibration record.

10. A particulate matter size pollution diffusion simulation system based on mobile reconstruction, used to execute the particulate matter size pollution diffusion simulation method based on mobile reconstruction as described in any one of claims 1 to 9, characterized in that, include: The sampling chain identification module is used to acquire particle size detection data, positioning data, inertial data, meteorological data, and sampling chain operating condition data, and determine the inlet sampling efficiency and discrete delay-broadening transfer kernel according to particle size; the wind field reconstruction module is used to generate the reconstructed wind field of the underway area based on the meteorological data, positioning data, and inertial data. The spatiotemporal observation band construction module is used to map each core delay unit to the navigation position, heading and local wind field in the reconstructed wind field at the historical actual sampling time, so as to obtain the particle size spatiotemporal observation band; The module constructs a source-receptor matrix to perform inverse Lagrange multiplication and obtain the particle size source-receptor matrix. It also includes a source emission rate inversion module to estimate the particle size source emission rate vector under non-negativity constraints, spatial sparsity constraints, and cross-particle size common source activation constraints. A diffusion simulation module performs forward particle diffusion and obtains the concentration fields for each particle size, the total particulate matter concentration field, the pollution center, diffusion direction, particle size composition, sedimentation flux, sedimentation range, and pollution boundary. A closed-loop calibration module, including a low-hysteresis bypass, a main sampling path, a switching component, and a multi-particle size calibration pulse component, identifies the discrete delay-broadening transfer kernel and kernel covariance for each particle size. It processes the forward particle diffusion results through inlet sampling efficiency, discrete delay-broadening transfer kernel, and particle size spatiotemporal observation bands to obtain the predicted detection response. Based on the particle-by-particle size residuals of the predicted detection response, it determines sampling chain drift and updates the discrete delay-broadening transfer kernel and particle size source emission rate vector or reconstructs the wind field for the target particle size based on the determination results. Finally, it receives updates or performs rollbacks based on the residuals and kernel covariance before and after the update.

Images