Disaster area unmanned aerial vehicle light source cooperative lighting advanced process control method

Abstract

The present application relates to the technical field of program control, and discloses an advanced process control method for disaster area unmanned aerial vehicle light source cooperative lighting. In view of the light intensity superposition distortion problem in the light field overlapping area of multiple machines, the unmanned aerial vehicle pose change sequence and the light source optical axis pointing angle change sequence under strong wind load are obtained, and a multivariate non-minimum phase transfer function matrix is constructed with the pose change as input and the ground grid node light intensity change as output. A diagonal matrix decoupling compensator is introduced to decouple the above-mentioned predicted output vector into independent single variable channels, and the compensation angular velocity control sequence of the light source holder is solved by optimizing the target of approximating the preset light intensity value with the predicted value. The sequence is reversely superimposed with the bottom attitude angular velocity to drive the holder to execute compensation. The scheme cuts off the cross-coupling path of the pose disturbance and the adjacent light beams, eliminates the light intensity flicker and overexposure in the overlapping area, and maintains the static consistency of the light distribution on the rescue operation surface.

Description

Technical Field

[0001] This invention relates to the field of program control technology, and discloses an advanced process control method for collaborative lighting using UAV light sources in disaster areas. Background Technology

[0002] In the field of collaborative lighting control using UAVs in disaster areas, existing conventional solutions typically employ a serial control architecture that separates flight control from gimbal stabilization. Specifically, the UAV flight control system receives attitude data from the inertial measurement unit and uses a proportional-integral-derivative (PID) control law to output rotation speed commands and adjust rotor thrust to maintain fuselage stability. Simultaneously, the gimbal, equipped with an independent gyroscope sensor, suppresses angular velocity disturbances transmitted from the fuselage through internal attitude calculation and isolation algorithms. In multi-UAV collaborative lighting scenarios, the ground control station allocates fixed lighting areas to each UAV based on a pre-defined geometric bounding box. Each UAV then executes the aforementioned independent attitude stabilization and gimbal isolation control according to the boundaries of its assigned area to maintain the projection direction of the light beam within the allocated area.

[0003] The core flaw of the aforementioned conventional technical solutions lies in the inherent response lag in the attitude correction of the flight control system, and the fact that the independent gimbal stabilization algorithm only suppresses disturbances on a single aircraft body, failing to incorporate the spatial cross-coupling effect formed by adjacent beams of multiple UAVs on the ground into the control loop. Under the continuous complex disturbances caused by strong winds in the disaster area, the beams of adjacent UAVs will produce synchronous swaying in the overlapping area. Due to the lack of coordinated decoupling control methods for multivariate illuminance distribution, the illuminance in the overlapping area will produce severe superimposed distortion and flicker, resulting in an unstable light field phenomenon with alternating local overexposed areas and blind spots, causing the illumination distribution on the rescue operation surface to fail to maintain static consistency. Summary of the Invention

[0004] The purpose of this invention is to provide an advanced process control method for drone-based coordinated lighting in disaster areas, which can solve the problems mentioned in the background art.

[0005] To achieve the above objectives, the technical solution adopted by the present invention is as follows: An advanced process control method for UAV-based coordinated lighting in disaster areas includes: acquiring a six-degree-of-freedom pose change sequence of a UAV under strong wind load in the disaster area and a corresponding light source optical axis pointing angle change sequence; and constructing a multivariable non-minimum phase transfer function matrix with the pose change sequence as input and the illuminance change sequence of a preset ground grid node as output through a system identification algorithm. Within the dynamic matrix control architecture, the step response model of the multivariable non-minimum phase transfer function matrix is ​​extracted, and the illuminance prediction output vector in the current time and future prediction time domain is calculated based on the step response model. A diagonal matrix decoupling compensator is introduced to decouple the illuminance prediction output vector into mutually independent univariate illuminance prediction channels; In each control cycle, the optimization objective is to approximate the preset illuminance value with the predicted value of each decoupled illuminance prediction channel. The compensation angular velocity control sequence of the light source gimbal in pitch and yaw degrees of freedom is calculated using a non-parametric model optimization algorithm. The compensation angular velocity control sequence is superimposed in reverse with the underlying attitude angular velocity of the UAV, driving the light source gimbal to perform reverse angular velocity compensation.

[0006] Preferably, a multivariable non-minimum phase transfer function matrix is ​​constructed by a system identification algorithm, with the pose change sequence as input and the illuminance change sequence of the ground preset grid nodes as output, including: injecting a pseudo-random binary sequence superimposed with multiple frequency bands into the UAV flight control system as an excitation signal, and simultaneously collecting the pose response data of the UAV under the action of the excitation signal and the illuminance response data of the ground grid nodes; Using the pose response data as the input vector and the illumination response data as the output vector, a model to be identified in the form of a difference equation is constructed. The parameter matrix in the model to be identified is solved iteratively using the recursive least squares method with a forgetting factor. The iteration stops when the parameter change is lower than the preset convergence threshold for three consecutive sampling periods. The solved parameter matrix is ​​then converted into the multivariable non-minimum phase transfer function matrix.

[0007] Preferably, the calculation of the illuminance prediction output vector in the current time and future prediction time domain based on the step response model includes: extracting the measured step response coefficients within a historical time window of a set length before the current time from the step response model, and constructing the basic prediction vector for dynamic matrix control. Obtain the illuminance deviation vector between the actual illuminance measurement value of the ground preset grid node at the current moment and the illuminance prediction value calculated in the previous control cycle; The illuminance deviation vector is weighted and corrected using a pre-defined error correction coefficient matrix. The corrected illuminance deviation vector is then added to the basic prediction vector to generate the illuminance prediction output vector, which includes error feedback correction features.

[0008] Preferably, a diagonal matrix decoupling compensator is introduced to decouple the illuminance prediction output vector into mutually independent univariate illuminance prediction channels, including: performing an inversion operation on the multivariate non-minimum phase transfer function matrix within a preset frequency band to obtain a generalized inverse matrix; Extract the off-diagonal elements of the generalized inverse matrix to construct the feedforward decoupling compensation matrix; The feedforward decoupling compensation matrix is ​​multiplied with the illuminance prediction output vector to generate the univariate illuminance prediction channel that contains only the channel information corresponding to the diagonal elements.

[0009] Preferably, the compensation angular velocity control sequence of the light source gimbal in pitch and yaw degrees of freedom is calculated using a nonparametric model optimization algorithm, including: constructing a quadratic objective function with the sum of squared deviations between the predicted value of the univariate illuminance prediction channel and the preset illuminance value as the optimized variable, and introducing a control increment weighting matrix for the compensation angular velocity control sequence and an output error weighting matrix for the sum of squared deviations into the quadratic objective function; Under the constraints of the control increment weighting matrix and the output error weighting matrix, the partial derivatives of the quadratic objective function with respect to the compensation angular velocity control sequence are calculated and set to 0. The compensation angular velocity control sequence is then calculated by solving the inverse of the Hessian matrix.

[0010] Preferably, the compensation angular velocity control sequence is superimposed in reverse with the underlying attitude angular velocity of the UAV to drive the light source gimbal to perform reverse angular velocity compensation, including: extracting the compensation angular velocity component of the first control cycle in the compensation angular velocity control sequence and taking the negative value of the compensation angular velocity component; The negative value of the compensated angular velocity component is used as the feedforward compensation amount, and algebraically added to the underlying attitude angular velocity calculated by the UAV inertial measurement unit to generate a comprehensive angular velocity command. The integrated angular velocity command is input to the motor driver of the light source gimbal. The motor driver outputs the corresponding three-phase current according to the integrated angular velocity command to drive the gimbal motor to operate in order to eliminate the influence of pose disturbance on the light field distribution.

[0011] Preferably, before injecting the multi-band superimposed pseudo-random binary sequence as an excitation signal into the UAV flight control system, the method further includes: performing wavelet packet decomposition on the collected original pose change sequence of the UAV under strong wind load, selecting the Daubechies wavelet basis function to perform multi-scale decomposition on the original pose change sequence, and obtaining wavelet packet decomposition coefficients of multiple frequency bands. Calculate the Shannon entropy value of the wavelet packet decomposition coefficients of each frequency band, and determine the frequency band with the Shannon entropy value greater than the preset entropy value threshold as the wind noise dominant frequency band, and determine the frequency band with the Shannon entropy value less than or equal to the preset entropy value threshold as the effective pose dominant frequency band. Set the wavelet packet decomposition coefficients corresponding to the dominant frequency band of wind load noise to zero, and use the wavelet packet decomposition coefficients corresponding to the dominant frequency band of effective pose to perform wavelet packet reconstruction, and output the pose change sequence that filters out the high-frequency flutter noise of wind load.

[0012] Preferably, the illuminance deviation vector is weighted and corrected using a pre-set error correction coefficient matrix, including: obtaining a time-varying softening coefficient sequence of the dynamic matrix control architecture in the future prediction time domain, and generating a diagonal softening matrix using the time-varying softening coefficient sequence; Calculate the matching degree between the current illuminance prediction time domain length and the dimension of the error correction coefficient matrix. When the matching degree indicates a change in the prediction step size, perform a shift and cyclic filling operation on the reference error correction coefficients according to the changed prediction step size to generate the error correction coefficient matrix. The error correction coefficient matrix and the softening matrix are subjected to Hadamard product operation, and the resulting matrix is ​​used to perform weighted correction on the illuminance deviation vector.

[0013] Preferably, the generalized inverse matrix is ​​obtained by performing an inversion operation on the multivariable non-minimum phase transfer function matrix within a preset frequency band, including: collecting the pitch angle change frequency and yaw angle change frequency of the UAV in the current control cycle, and determining the maximum value of the pitch angle change frequency and the yaw angle change frequency as the current main frequency disturbance frequency; Based on the current dominant frequency perturbation frequency, the adaptive forgetting factor value corresponding to the current dominant frequency perturbation frequency is determined by a lookup table method; The off-diagonal elements in the multivariate non-minimum phase transfer function matrix are attenuated and weighted using the adaptive forgetting factor value. The inversion operation is then performed on the attenuated and weighted multivariate non-minimum phase transfer function matrix to generate the generalized inverse matrix that weakens the off-diagonal coupling strength.

[0014] Preferably, after algebraically adding the negative value of the compensated angular velocity component as the feedforward compensation amount to the bottom attitude angular velocity calculated by the UAV inertial measurement unit, the method further includes: determining whether the absolute value of the comprehensive angular velocity command generated after algebraic addition is less than the static friction angular velocity threshold of the light source gimbal motor. When the absolute value of the comprehensive angular velocity command is less than the static friction angular velocity threshold, a dead zone compensation voltage pulse sequence consistent with the direction sign of the comprehensive angular velocity command is output to the motor driver. The amplitude of the dead zone compensation voltage pulse sequence is determined by linear interpolation based on the difference between the static friction angular velocity threshold and the comprehensive angular velocity command, so as to drive the light source gimbal motor to break through the static friction state and perform micro-motion compensation.

[0015] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. This invention establishes a multivariable non-minimum phase transfer function matrix with pose change as input and illuminance change of ground grid nodes as output. A diagonal matrix decoupling compensator is introduced into the dynamic matrix control architecture to decouple the illuminance prediction output vector into independent single-variable channels for optimization, thus cutting off the cross-coupling transmission path between the UAV's six-degree-of-freedom pose disturbance and adjacent beams from multiple UAVs. Combined with an illuminance prediction output vector calculation method with error feedback correction, the actual illuminance deviation of the ground grid nodes is directly used as the optimization target to solve the gimbal compensation angular velocity control sequence. This eliminates flickering and overexposure phenomena caused by illuminance superposition distortion in the overlapping area of ​​multiple UAV light fields under strong wind interference, maintaining a statically consistent illuminance distribution on the ground rescue operation surface.

[0016] 2. By performing wavelet packet decomposition and Shannon entropy calculation on the pose change sequence before system identification, high-frequency flutter noise from wind loads is filtered out, improving the modeling accuracy of the multivariable transfer function matrix. The adaptive forgetting factor value is determined by looking up the table based on the current main frequency disturbance frequency, and the off-diagonal elements are attenuated and weighted, weakening the residual coupling strength between channels under time-varying frequencies. When the comprehensive angular velocity command is in the static friction range of the motor, the dead zone compensation voltage pulse sequence is output based on the difference linear interpolation, overcoming the static friction hysteresis of the light source gimbal motor when performing micro-motion compensation, and ensuring the physical execution continuity of the reverse angular velocity compensation command under small disturbance conditions. Attached Figure Description

[0017] Figure 1 An overall flowchart of an advanced process control method for collaborative lighting of unmanned aerial vehicles (UAVs) in disaster areas, provided by an embodiment of the present invention; Figure 2 A flowchart for constructing a multivariable non-minimum phase transfer function matrix provided in an embodiment of the present invention; Figure 3 A flowchart for calculating the illuminance prediction output vector including error feedback correction is provided for embodiments of the present invention; Figure 4 This is a flowchart of channel decoupling using a diagonal matrix decoupling compensator provided in an embodiment of the present invention; Figure 5 This is a flowchart of solving the compensated angular velocity control sequence provided in an embodiment of the present invention; Figure 6 This is a flowchart illustrating the signal preprocessing and motor static friction compensation process provided in an embodiment of the present invention. Detailed Implementation

[0018] This invention relates to the field of drone lighting control technology, specifically to an advanced process control method for drone-based coordinated lighting in disaster areas. The invention will be further described in detail below with reference to embodiments. These specific embodiments are for illustrative purposes only and are not intended to limit the scope of the invention.

[0019] Please refer to Figure 1 This embodiment provides an advanced process control method for UAV-based coordinated lighting in disaster areas, acquiring the UAV's six-degree-of-freedom pose change sequence and the corresponding light source optical axis pointing angle change sequence under strong wind loads in disaster areas. Specifically, based on the sampling period of the UAV's onboard inertial measurement unit, an equally spaced sampling time sequence is set, with the sampling period consistent with the closed-loop control period of the UAV's flight control system. At each sampling time k, the UAV's six-degree-of-freedom pose data in the northeast-central geographic coordinate system is acquired through the combined positioning module of the onboard inertial measurement unit and the global satellite navigation system. The six-degree-of-freedom pose data includes three translational degrees of freedom position data along the X, Y, and Z axes, as well as the pitch angle. Yaw angle Roll angle Using the rotational degrees of freedom angle data, construct the pose vector X(k) at sampling time k, as shown in the formula: in, , , These represent the position coordinates of the UAV along the X, Y, and Z axes in the northeast-northeast coordinate system at sampling time k. , , These represent the pitch angle, yaw angle, and roll angle of the UAV at sampling time k, respectively, with the superscript T indicating matrix transpose.

[0020] Furthermore, the pose increments at adjacent sampling times are calculated to construct a pose change sequence. As shown in the formula: in, Let be the pose change vector at sampling time k. Let k-1 be the pose vector at the previous sampling time. The pose change sequence is composed of continuous... The pose change vectors at each sampling time point constitute the structure. The sequence length is given.

[0021] Synchronously, at each sampling time k, the pitch angle of the light source's optical axis relative to the UAV's body coordinate system is acquired through the angle measurement unit built into the light source gimbal. With yaw heading angle Construct the optical axis pointing angle vector Calculate the optical axis pointing angle increment at adjacent sampling times to construct the corresponding light source optical axis pointing angle change sequence. As shown in the formula: in, , These represent the changes in pitch and yaw pointing angles of the light source's optical axis at sampling time k, respectively.

[0022] In this embodiment, for a multi-drone collaborative lighting scenario, the number of drones involved in the collaborative operation is set to [number missing]. Each drone simultaneously collects its corresponding pose change sequence and optical axis pointing angle change sequence, constructing a multi-input sequence set. For the ground rescue operation surface, uniformly distributed grid nodes are pre-divided, with the number of grid nodes being [number missing]. Each grid node corresponds to a fixed ground coordinate. ,in z=0 is the plane where the ground is located. The actual illuminance values ​​of each grid node at each sampling time k are synchronously collected through an array of illuminance sensors deployed on the ground. Calculate the illuminance increment between adjacent sampling times. Construct a sequence of illuminance changes at pre-defined ground grid nodes. .

[0023] Furthermore, a multivariable non-minimum phase transfer function matrix is ​​constructed using a system identification algorithm, taking the pose change sequence as input and the ground preset grid node illumination change sequence as output. Specifically, the matrix is... The pose change sequence of the UAV and the optical axis pointing angle change sequence are merged into the system input vector. Input vector The dimension is The illuminance variation sequence of ground grid nodes is used as the system output vector. The dimension of the output vector Y(k) is A multi-input multi-output system model in the discrete-time domain is constructed, and the transfer function matrix of the system is solved using a system identification algorithm. The transfer function matrix The matrix is ​​a multivariable non-minimum phase transfer function with dimension . As shown in formula (4): in, Let z be the single-input single-output transfer function between the i-th output channel and the j-th input channel, and z be the z-transform operator in the discrete-time domain. The non-minimum phase characteristic of the transfer function matrix G(z) is reflected in a partial transfer function. The zero point is located outside the unit circle in the z-plane, corresponding to the inverse response characteristics of the UAV attitude change to the degree output under strong wind disturbance.

[0024] In this embodiment, within the dynamic matrix control architecture, the step response model of the multivariable non-minimum phase transfer function matrix is ​​extracted, and the illuminance prediction output vector in the current time and future prediction time domain is calculated based on the step response model. Specifically, for each single-input single-output transfer function in the transfer function matrix G(z)... Solve the unit step response sequence to obtain the step response coefficient sequence corresponding to each transfer function. Where L is the time-domain length of the model, and the value of L is greater than that of the transfer function. The number of sampling steps corresponding to the adjustment time. Based on the step response coefficient sequence of all channels, a dynamic matrix A of the dynamic matrix control architecture is constructed. The dimension of dynamic matrix A is... Where P is the prediction time domain length and M is the control time domain length. The element construction rule of dynamic matrix A is: for the j-th input channel and the i-th output channel, the element of dynamic matrix A is... line, number The elements of the column are The value is 0 if p ≥ m, otherwise it is 0, where p is the prediction step size index and m is the control step size index.

[0025] Furthermore, measured step response coefficients within a historical time window of a set length prior to the current time k are extracted from the step response model to construct the basic prediction vector for dynamic matrix control. Basic prediction vector This is a sequence of predicted illuminance values ​​for P future sampling times without control increments, with dimension 1. Obtain the actual illuminance measurement value Y(k) of the ground preset grid node at the current time k and the illuminance prediction value calculated in the previous control cycle. Illuminance deviation vector between The illuminance deviation vector e(k) is weighted and corrected using a pre-defined error correction coefficient matrix H. The corrected illuminance deviation vector is then compared with the base prediction vector. The summation generates an illuminance prediction output vector that includes error feedback correction features. As shown in the formula: in, Let k be the output vector of predicted illuminance for the next P time domains calculated at the current time k, with dimension k. H is the error correction coefficient matrix with dimension 1. , Let be the illuminance deviation vector at the current moment, with dimension . .

[0026] In this embodiment, a diagonal matrix decoupling compensator is introduced to decouple the illuminance prediction output vector into mutually independent univariate illuminance prediction channels. Specifically, a compensation matrix D is constructed for the diagonal matrix decoupling compensator, and the dimension of the compensation matrix D is... The steady-state gain matrix G(z) of the multivariable non-minimum phase transfer function matrix is ​​solved within a preset operating frequency band to obtain the steady-state gain matrix G(1), where z=1 corresponds to the steady-state operating point of the discrete system. The generalized inverse operation of the steady-state gain matrix G(1) is then performed to obtain the generalized inverse matrix. Based on generalized inverse matrix Construct a feedforward decoupling compensation matrix, and then combine the feedforward decoupling compensation matrix with the illuminance prediction output vector. Matrix multiplication is performed to eliminate cross-coupling components between channels of varying optical axis pointing angles of different UAV light sources, generating a univariate illuminance prediction channel that contains only channel information corresponding to diagonal elements. As shown in the formula: in, is the set of decoupled univariate illuminance prediction channels. Each channel corresponds to the illuminance prediction sequence of a ground grid node. There are no cross-coupling components between channels. D is the compensation matrix of the diagonal matrix decoupling compensator. It is a diagonal matrix, and all off-diagonal elements are 0.

[0027] In this embodiment, in each control cycle, the optimization objective is to approximate the preset illuminance value with the predicted values ​​of each decoupled illuminance prediction channel. A non-parametric model optimization algorithm is used to calculate the compensation angular velocity control sequence for the light source gimbal in pitch and yaw degrees of freedom. Specifically, a corresponding preset illuminance value sequence is set for each decoupled univariate illuminance prediction channel. , Let P be the future predicted target illuminance value sequences in the time domain, with dimension P. A quadratic objective function is constructed with the sum of squared deviations between the predicted values ​​and preset illuminance values ​​from the univariate illuminance prediction channel as the optimized variable. A control increment weighting matrix for the compensated angular velocity control sequence and an output error weighting matrix for the sum of squared deviations are introduced into the objective function. By solving for the extrema of the objective function, the control increment sequence in the control time domain is obtained. The control increment sequence This includes the compensated angular velocity control sequence for each drone's light source gimbal in the pitch and yaw degrees of freedom, with dimensions of... .

[0028] Furthermore, the compensated angular velocity control sequence is superimposed in reverse with the underlying attitude angular velocity of the UAV, driving the light source gimbal to perform reverse angular velocity compensation. Specifically, the control increment sequence is extracted. The compensated angular velocity component of the first control step corresponding to the current control cycle, wherein the compensated angular velocity component includes the pitch degree of freedom compensated angular velocity. Angular velocity with yaw degree of freedom compensation By taking the negative value of the compensated angular velocity component, the reverse compensated angular velocity component is obtained. , The reverse compensation angular velocity component is used as the feedforward compensation amount and algebraically added to the bottom attitude angular velocity calculated by the UAV inertial measurement unit to generate a comprehensive angular velocity command. The comprehensive angular velocity command is input to the motor driver of the light source gimbal. The motor driver outputs the corresponding drive signal according to the comprehensive angular velocity command to drive the gimbal motor to operate, thereby offsetting the influence of UAV attitude disturbance on the optical axis pointing of the light source and maintaining the illumination stability of the ground grid nodes.

[0029] In this embodiment, the correspondence between the coordinates of the ground preset grid nodes and the preset illuminance values ​​is shown in Table 1. Table 1 is a table of preset illuminance parameters for the grid nodes of the ground rescue operation surface. The grid nodes are divided into a uniform grid of 5 rows and 5 columns, covering a square rescue operation surface. The grid nodes are evenly spaced, and the preset illuminance value of each node is set according to the rescue operation requirements. The preset illuminance value of the core operation area is higher than that of the edge area.

[0030] Table 1 Preset Illumination Parameters for Grid Nodes on Ground Rescue Operation Surface In this embodiment, by constructing a multivariable non-minimum phase transfer function matrix, a precise mapping relationship between UAV pose changes and ground illumination changes is established. Multi-step prediction and feedback correction of illumination are realized within the dynamic matrix control architecture. Cross-coupling between multiple UAV channels is eliminated through a diagonal matrix decoupling compensator. The gimbal compensation control sequence is solved directly with the prediction deviation of ground illumination as the optimization objective. This achieves coordinated lighting control of multiple UAV light sources under strong wind disturbances and maintains stable illumination distribution on the ground rescue operation surface.

[0031] refer to Figure 2 In one optional embodiment, the original pose change sequence of the acquired UAV under strong wind load is preprocessed. High-frequency flutter noise from the wind load is filtered out through wavelet packet decomposition and Shannon entropy calculation, outputting a denoised pose change sequence. Specifically, the Daubechies wavelet basis function is selected as the basis function for wavelet packet decomposition, and the original pose change sequence is processed. Perform multi-scale wavelet packet decomposition at scale j, where j is the decomposition level, to obtain the following results. The wavelet packet decomposition coefficients for each frequency band are given, and each frequency band corresponds to a wavelet packet subspace. The decomposition process is shown in formula (7). in, Let be the wavelet packet decomposition coefficients of the j-th scale and the l-th frequency band. h(n) are the low-pass decomposition filter coefficients of the Daubechies wavelet basis, g(n) are the high-pass decomposition filter coefficients, and n is the discrete time series index.

[0032] Furthermore, the Shannon entropy value of the wavelet packet decomposition coefficients for each frequency band is calculated, for the wavelet packet decomposition coefficient sequence of the j-th scale and the l-th frequency band. Calculate its normalized energy sequence ,in Given the total energy of this frequency band, calculate the Shannon entropy value of this frequency band based on the normalized energy sequence. As shown in the formula: Among them, when When, define To avoid meaningless logarithmic operations, a preset entropy threshold is set. Shannon entropy value Greater than the preset entropy threshold The frequency band was determined to be the dominant frequency band for wind noise, and the Shannon entropy value was used. Less than or equal to a preset entropy threshold The frequency band is determined as the dominant frequency band of the effective pose. All wavelet packet decomposition coefficients corresponding to the dominant wind noise frequency band are set to zero, while the wavelet packet decomposition coefficients corresponding to the dominant effective pose frequency band are retained. The processed wavelet packet decomposition coefficients are reconstructed using a wavelet packet reconstruction algorithm, outputting a pose change sequence filtered out of high-frequency wind flutter noise. .

[0033] In this embodiment, when constructing the multivariable non-minimum phase transfer function matrix using the system identification algorithm, a multi-band superimposed pseudo-random binary sequence is injected into the UAV flight control system as an excitation signal. Simultaneously, the pose response data of the UAV under the excitation signal and the illumination response data of the ground grid nodes are collected. Specifically, the pseudo-random binary sequence is generated using a linear feedback shift register. The multi-band superimposed pseudo-random binary sequence consists of three different frequency bands superimposed, corresponding to the low-frequency, mid-frequency, and high-frequency bands, respectively. The low-frequency band corresponds to the dominant frequency band of the UAV's rigid body motion, the mid-frequency band corresponds to the dominant frequency band of gimbal attitude adjustment, and the high-frequency band corresponds to the dominant frequency band of wind load disturbance. The parameters of the multi-band superimposed pseudo-random binary sequence are shown in Table 2.

[0034] Table 2 Excitation Parameters for Multi-Band Superimposed Pseudo-Random Binary Sequences The generated multi-band superimposed pseudo-random binary sequence is used as an excitation signal and injected into the attitude control loop of the UAV flight control system. The injection duration of the excitation signal covers the complete sequence period. During the entire process of excitation signal injection, the UAV's pose response data, gimbal optical axis pointing angle response data, and ground grid node illumination response data are collected synchronously at a fixed sampling period. The number of samples collected matches the total sequence length.

[0035] Furthermore, the pose response data and optical axis pointing angle response data are merged into an input vector U(k), and the illuminance response data is used as the output vector Y(k). A model to be identified in the form of a difference equation is constructed. The model to be identified adopts a controlled autoregressive moving average model structure, as shown in the formula: in The autoregressive order at the output end. The order of the moving average at the input end. for The output autoregressive parameter matrix has dimensions of 1. , for The input moving average parameter matrix has dimensions of 1. , The model residual vector has dimensions of . .

[0036] The recursive least squares method with a forgetting factor is used to iteratively solve the parameter matrix in the model to be identified. All parameter matrices to be identified are expanded into parameter vectors θ, and the corresponding regression vectors φ(k) are constructed. The iterative solution process is shown in the formula: Where K(k) is the iterative gain matrix, P(k) is the covariance matrix, λ is the forgetting factor with a value range of 0 < λ < 1, and I is the identity matrix. Let be the parameter vector estimate at time k, and the initial value for the iteration process be set to . P(0) is the product of the identity matrix and the preset initial coefficients.

[0037] During the iterative solution process, the parameter changes over three consecutive sampling periods are calculated in real time. When the parameter changes for three consecutive sampling periods are all below the preset convergence threshold Stop the iteration when the time is reached. Solve the resulting parameter matrix. , The transformation to a multivariable non-minimum phase transfer function matrix G(z) in the z-domain is achieved by performing a z-transform on the difference equation, resulting in the expression for the transfer function matrix: in, This is the shift operator.

[0038] In this embodiment, wavelet packet decomposition and Shannon entropy filtering are used to filter out high-frequency wind-borne flutter noise in the original pose sequence, improving the signal-to-noise ratio of the system's input data identification. By using a pseudo-random binary sequence excitation with multiple frequency bands superimposed, the entire operating frequency band of the UAV and gimbal system is covered. Combined with the recursive least squares method with a forgetting factor, the model parameters are quickly converged and accurately identified, ensuring the modeling accuracy of the multivariable non-minimum phase transfer function matrix.

[0039] refer to Figure 3 and Figure 4 In one optional embodiment, when calculating the illuminance prediction output vector for the current time and future prediction time domain based on the step response model, the measured step response coefficients within a historical time window of a set length prior to the current time k are extracted from the step response model to construct the basic prediction vector for dynamic matrix control. Specifically, the length of the historical time window is consistent with the model's time domain length L. Measured step response coefficients from L sampling times prior to the current time k are extracted and combined with the control increment sequence from M control cycles prior to the current time. Calculate the basic prediction vector without future control increments. Each element of the basic prediction vector corresponds to the illuminance prediction value at P future prediction times.

[0040] Obtain the actual illuminance measurement value Y(k) of the ground preset grid node at the current time k and the illuminance prediction value calculated in the previous control cycle. Illuminance deviation vector between The illuminance deviation vector is weighted and corrected using a pre-defined error correction coefficient matrix. Specifically, the time-varying softening coefficient sequence of the dynamic matrix control architecture in the future prediction time domain is obtained. The range of values ​​for the softening coefficient sequence is: , The softening coefficient sequence monotonically increases with the prediction step size. A diagonal softening matrix S is generated using the time-varying softening coefficient sequence. The dimensions of the diagonal softening matrix S are... The diagonal elements are repeating. The softening coefficient sequence of the next step .

[0041] The matching degree between the predicted illuminance temporal length P and the dimension of the error correction coefficient matrix H at the current moment is calculated. The matching degree is calculated as the ratio of the predicted temporal length P to the number of rows in the error correction coefficient matrix H. When the matching degree is not equal to 1, it is determined that the prediction step size has changed. Based on the changed prediction step size, a shift and cyclic filling operation is performed on the reference error correction coefficients to generate an error correction coefficient matrix H that matches the current predicted temporal length. The reference error correction coefficients are a one-dimensional sequence of length L. The elements of the baseline error correction coefficient sequence monotonically decrease as the sequence index increases. When the prediction step size changes to... At that time, the reference error correction coefficient sequence is cyclically shifted, and the first part is truncated. 1 element, repeated Next, the generation dimension is The error correction coefficient matrix H.

[0042] Perform a Hadamard product operation on the error correction coefficient matrix H and the diagonal softening matrix S. The Hadamard product operation involves multiplying corresponding elements of two matrices with the same dimensions. The resulting matrix is ​​the corrected error correction matrix. , where ⊙ is the Hadamard product operator. The corrected error correction matrix is ​​used. The illuminance deviation vector e(k) is weighted and corrected, and the corrected illuminance deviation vector is compared with the base prediction vector. The summation generates an illuminance prediction output vector that includes error feedback correction features. .

[0043] In this embodiment, when a diagonal matrix decoupling compensator is introduced to decouple the illuminance prediction output vector into mutually independent univariate illuminance prediction channels, the multivariate non-minimum phase transfer function matrix G(z) is inverted within a preset frequency band to obtain the generalized inverse matrix. Specifically, the pitch angle change frequency of the UAV within the current control cycle is collected. With yaw angle change frequency Pitch angle change frequency The yaw angle variation frequency is obtained by solving the dominant frequency using a Fast Fourier Transform (FFT) on a sequence of pitch angle changes with a continuously set number of sampling periods. Solve in the same way to get the pitch angle change frequency. With yaw angle change frequency The maximum value in is determined as the current main frequency perturbation frequency. .

[0044] Based on the current main frequency disturbance frequency The adaptive forgetting factor value corresponding to the current dominant frequency perturbation frequency is determined by looking up a table. The adaptive forgetting factor values ​​corresponding to different main frequency perturbation frequencies are shown in Table 3. Table 3 is a parameter table corresponding to the main frequency perturbation frequency and the adaptive forgetting factor. The frequency intervals in the table are divided according to logarithmic intervals. The adaptive forgetting factor value decreases as the main frequency perturbation frequency increases.

[0045] Table 3. Parameters corresponding to the dominant frequency perturbation frequency and the adaptive forgetting factor. The adaptive forgetting factor value obtained by looking up a table For the multivariable non-minimum phase transfer function matrix G(z) at the current dominant frequency perturbation frequency Frequency response matrix at ... The off-diagonal elements in the matrix are attenuated and weighted, resulting in the frequency response matrix after attenuation weighting. The elements are as shown in the formula: in, The element in the i-th row and j-th column of the frequency response matrix after attenuation weighting is... , This represents the adaptive forgetting factor value.

[0046] Frequency response matrix after attenuation weighting Perform the inversion operation to generate a generalized inverse matrix that weakens the off-diagonal coupling strength. The generalized inverse matrix is ​​solved using the Moore-Penrose generalized inverse. Off-diagonal elements of the generalized inverse matrix are extracted to construct a feedforward decoupling compensation matrix. This matrix is ​​then combined with the illuminance prediction output vector. Matrix multiplication is performed to eliminate the cross-coupling components between different UAV light source optical axis pointing angle variation channels, generating a univariate illuminance prediction channel that only contains channel information corresponding to diagonal elements. Each univariate illuminance prediction channel corresponds to the illuminance prediction sequence of only one ground grid node and has no coupling relationship with other channels.

[0047] In this embodiment, adaptive error correction is achieved when the predicted time domain length changes by using time-varying softening coefficients and shift-cycle-filled error correction coefficient matrices, which improves the accuracy of the illuminance prediction output vector. The adaptive forgetting factor is determined by looking up the table at the current dominant frequency disturbance frequency, and the off-diagonal coupling elements of the transfer function matrix are attenuated and weighted, which weakens the residual coupling strength between channels under time-varying wind load disturbance frequency, improves the decoupling compensation effect, and ensures the independence of the univariate illuminance prediction channel.

[0048] refer to Figure 5In one optional embodiment, when calculating the compensated angular velocity control sequence of the light source gimbal in pitch and yaw degrees of freedom using a nonparametric model optimization algorithm, a quadratic objective function is constructed with the sum of squared deviations between the predicted values ​​of the decoupled univariate illuminance prediction channel and the preset illuminance value as the optimized variable. Specifically, the predicted value vector of the decoupled univariate illuminance prediction channel is... The corresponding preset illuminance value sequence is R(k), where R(k) is the target illuminance value sequence for the next P predicted time domains, with dimension P. Construct a quadratic objective function J, as shown in formula (12): in, This is a weighted L2 norm operation, where Q is the output error weighting matrix with dimension 1. Let be a diagonal matrix, where the diagonal elements are the illumination error weights for each grid node at each prediction time. The weights of grid nodes in the core operating area are higher than those in the edge area. R is the control increment weighting matrix with dimension . , is a diagonal matrix, with diagonal elements being the weight values ​​for each control increment channel, used to limit the variation of the control increment and avoid overloading the gimbal motor.

[0049] Under the constraints of the control increment weighting matrix R and the output error weighting matrix Q, the quadratic objective function J is subjected to control sequence with respect to the compensated angular velocity. By calculating the partial derivatives and setting them to zero, the necessary condition for the objective function to reach its minimum is obtained. The compensation angular velocity control sequence is then calculated by solving for the inverse of the Hessian matrix. The analytical solution is shown in the formula: Where A is the dynamic matrix, This is the base prediction vector after decoupling. To control the compensated angular velocity control sequence in the time domain M, it includes the compensated angular velocity increments of each UAV light source gimbal in the pitch and yaw degrees of freedom. The compensated angular velocity corresponding to each control step is the sum of all control increments before that step.

[0050] refer to Figure 6 Furthermore, the compensated angular velocity control sequence is superimposed in reverse with the underlying attitude angular velocity of the UAV, driving the light source gimbal to perform reverse angular velocity compensation. Specifically, the compensated angular velocity control sequence is extracted. The compensated angular velocity component of the first control cycle, wherein the compensated angular velocity component is the compensated angular velocity of the gimbal pitch degree of freedom in the current control cycle. Compensated angular velocity with yaw degree of freedom By taking the negative value of the compensated angular velocity component, the reverse compensated angular velocity component is obtained. , The UAV's three-axis angular velocity data are collected by its onboard inertial measurement unit. Attitude calculations are then performed to obtain the UAV's underlying attitude angular velocity, which includes the fuselage pitch angular velocity. yaw rate Roll angular velocity The negative value of the compensated angular velocity component is used as the feedforward compensation amount, and algebraically added to the corresponding pitch and yaw angular velocities in the UAV's underlying attitude angular velocities to generate a comprehensive angular velocity command, as shown in the formula: in, For the combined angular velocity command of the pitch degree of freedom, This is the combined angular velocity command for the yaw degree of freedom.

[0051] The integrated angular velocity command is input to the motor driver of the light source gimbal. Based on the integrated angular velocity command, the motor driver outputs the corresponding drive current through the space vector pulse width modulation algorithm to drive the gimbal motor to run, perform reverse angular velocity compensation, and counteract the influence of UAV body attitude disturbance on the optical axis pointing of the light source, thus maintaining the illumination stability of the ground grid nodes.

[0052] In this embodiment, after algebraically adding the negative value of the compensated angular velocity component to the underlying attitude angular velocity of the UAV as the feedforward compensation amount, dead-zone compensation processing is also performed on the comprehensive angular velocity command. Specifically, the static friction angular velocity threshold of the light source gimbal motor is obtained in advance through gimbal motor test bench calibration. Static friction angular velocity threshold To determine the angular velocity value corresponding to the minimum angular velocity command required for the motor to transition from static friction to dynamic friction, the absolute value of the composite angular velocity command generated by algebraic summation is determined. Is it less than the static friction angular velocity threshold? ,in This is a combined angular velocity command for pitch or yaw degrees of freedom.

[0053] When the absolute value of the combined angular velocity command Less than the static friction angular velocity threshold At that time, based on the direction sign of the comprehensive angular velocity command, a dead-zone compensation voltage pulse sequence with the same direction sign is output to the motor driver. The amplitude of the dead-zone compensation voltage pulse sequence is based on the static friction angular velocity threshold. Combined angular velocity command The difference is determined by linear interpolation, as shown in the formula: in, The amplitude of the dead-zone compensation voltage pulse sequence, The minimum driving voltage amplitude required for the motor to overcome static friction is given by sign(·), where sign(·) is the sign function. hour, ,when hour, ,when hour, .

[0054] The pulse width of the dead-zone compensation voltage pulse sequence is set to one control cycle. When the comprehensive angular velocity command remains within the static friction range, the dead-zone compensation voltage pulse sequence is output once per control cycle. This drives the light source gimbal motor to break through the static friction state and perform micro-motion compensation, eliminating the influence of the motor's static friction hysteresis on the micro-compensation command. The corresponding parameters of the static friction angular velocity threshold and the dead-zone compensation voltage of the light source gimbal motor are shown in Table 4.

[0055] Table 4 Static Friction Characteristics and Dead Zone Compensation Parameters of the Light Source Gimbal Motor In this embodiment, the gimbal compensation angular velocity control sequence that satisfies the constraints is obtained by analytically solving the quadratic objective function. By superimposing the compensation angular velocity with the UAV's underlying attitude angular velocity in reverse, feedforward compensation for the body attitude disturbance is achieved. By judging the static friction range and linearly interpolating the dead zone compensation voltage pulse sequence, the static friction hysteresis of the light source gimbal motor when performing micro-motion compensation is overcome, ensuring the physical execution continuity of the reverse angular velocity compensation command under small disturbance conditions and improving the lighting control accuracy under weak disturbance conditions.

Claims

1. An advanced process control method for coordinated lighting using drone light sources in disaster areas, characterized in that, include: The six-degree-of-freedom pose change sequence of the UAV under strong wind load in the disaster area and the corresponding light source optical axis pointing angle change sequence are obtained. A multivariable non-minimum phase transfer function matrix is ​​constructed with the pose change sequence as input and the ground preset grid node illumination change sequence as output through the system identification algorithm. Within the dynamic matrix control architecture, the step response model of the multivariable non-minimum phase transfer function matrix is ​​extracted, and the illuminance prediction output vector in the current time and future prediction time domain is calculated based on the step response model. A diagonal matrix decoupling compensator is introduced to decouple the illuminance prediction output vector into mutually independent univariate illuminance prediction channels; In each control cycle, the optimization objective is to approximate the preset illuminance value with the predicted value of each decoupled illuminance prediction channel. The compensation angular velocity control sequence of the light source gimbal in pitch and yaw degrees of freedom is calculated using a non-parametric model optimization algorithm. The compensation angular velocity control sequence is superimposed in reverse with the underlying attitude angular velocity of the UAV, driving the light source gimbal to perform reverse angular velocity compensation.

2. The advanced process control method for coordinated lighting of UAV light sources in disaster areas according to claim 1, characterized in that, The system identification algorithm constructs a multivariable non-minimum phase transfer function matrix with the pose change sequence as input and the illuminance change sequence of the ground preset grid nodes as output. This includes: injecting a pseudo-random binary sequence with multiple frequency bands superimposed as an excitation signal into the UAV flight control system, and simultaneously collecting the pose response data of the UAV under the action of the excitation signal and the illuminance response data of the ground grid nodes. Using the pose response data as the input vector and the illumination response data as the output vector, a model to be identified in the form of a difference equation is constructed. The parameter matrix in the model to be identified is solved iteratively using the recursive least squares method with a forgetting factor. The iteration stops when the parameter change is lower than the preset convergence threshold for three consecutive sampling periods. The solved parameter matrix is ​​then converted into the multivariable non-minimum phase transfer function matrix.

3. The advanced process control method for coordinated lighting of UAV light sources in disaster areas according to claim 1, characterized in that, The calculation of the illuminance prediction output vector in the current time and future prediction time domain based on the step response model includes: extracting the measured step response coefficients within a historical time window of a set length before the current time from the step response model, and constructing the basic prediction vector for dynamic matrix control. Obtain the illuminance deviation vector between the actual illuminance measurement value of the ground preset grid node at the current moment and the illuminance prediction value calculated in the previous control cycle; The illuminance deviation vector is weighted and corrected using a pre-defined error correction coefficient matrix. The corrected illuminance deviation vector is then added to the basic prediction vector to generate the illuminance prediction output vector, which includes error feedback correction features.

4. The advanced process control method for coordinated lighting of UAV light sources in disaster areas according to claim 1, characterized in that, A diagonal matrix decoupling compensator is introduced to decouple the illuminance prediction output vector into mutually independent univariate illuminance prediction channels, including: performing an inversion operation on the multivariate non-minimum phase transfer function matrix in a preset frequency band to obtain a generalized inverse matrix; Extract the off-diagonal elements of the generalized inverse matrix to construct the feedforward decoupling compensation matrix; The feedforward decoupling compensation matrix is ​​multiplied with the illuminance prediction output vector to generate the univariate illuminance prediction channel that contains only the channel information corresponding to the diagonal elements.

5. The advanced process control method for coordinated lighting of UAV light sources in disaster areas according to claim 1, characterized in that, The compensation angular velocity control sequence of the light source gimbal in pitch and yaw degrees of freedom is calculated using a nonparametric model optimization algorithm, including: constructing a quadratic objective function with the sum of squared deviations between the predicted value of the univariate illuminance prediction channel and the preset illuminance value as the optimized variable, and introducing a control increment weighting matrix for the compensation angular velocity control sequence and an output error weighting matrix for the sum of squared deviations into the quadratic objective function; Under the constraints of the control increment weighting matrix and the output error weighting matrix, the partial derivatives of the quadratic objective function with respect to the compensation angular velocity control sequence are calculated and set to 0. The compensation angular velocity control sequence is then calculated by solving the inverse of the Hessian matrix.

6. The advanced process control method for coordinated lighting of UAV light sources in disaster areas according to claim 1, characterized in that, The compensation angular velocity control sequence is superimposed in reverse with the underlying attitude angular velocity of the UAV to drive the light source gimbal to perform reverse angular velocity compensation, including: extracting the compensation angular velocity component of the first control cycle in the compensation angular velocity control sequence and taking the negative value of the compensation angular velocity component; The negative value of the compensated angular velocity component is used as the feedforward compensation amount, and algebraically added to the underlying attitude angular velocity calculated by the UAV inertial measurement unit to generate a comprehensive angular velocity command. The integrated angular velocity command is input to the motor driver of the light source gimbal. The motor driver outputs the corresponding three-phase current according to the integrated angular velocity command to drive the gimbal motor to operate in order to eliminate the influence of pose disturbance on the light field distribution.

7. The advanced process control method for coordinated lighting of UAV light sources in disaster areas according to claim 2, characterized in that, Before injecting a multi-band superimposed pseudo-random binary sequence as an excitation signal into the UAV flight control system, the process further includes: performing wavelet packet decomposition on the original pose change sequence of the UAV under strong wind load, selecting the Daubechies wavelet basis function to perform multi-scale decomposition on the original pose change sequence, and obtaining wavelet packet decomposition coefficients of multiple frequency bands. Calculate the Shannon entropy value of the wavelet packet decomposition coefficients of each frequency band, and determine the frequency band with the Shannon entropy value greater than the preset entropy value threshold as the wind noise dominant frequency band, and determine the frequency band with the Shannon entropy value less than or equal to the preset entropy value threshold as the effective pose dominant frequency band. Set the wavelet packet decomposition coefficients corresponding to the dominant frequency band of wind load noise to zero, and use the wavelet packet decomposition coefficients corresponding to the dominant frequency band of effective pose to perform wavelet packet reconstruction, and output the pose change sequence that filters out the high-frequency flutter noise of wind load.

8. The advanced process control method for coordinated lighting of UAV light sources in disaster areas according to claim 3, characterized in that, The illuminance deviation vector is weighted and corrected using a pre-set error correction coefficient matrix, including: obtaining a time-varying softening coefficient sequence of the dynamic matrix control architecture in the future prediction time domain, and generating a diagonal softening matrix using the time-varying softening coefficient sequence; Calculate the matching degree between the current illuminance prediction time domain length and the dimension of the error correction coefficient matrix. When the matching degree indicates a change in the prediction step size, perform a shift and cyclic filling operation on the reference error correction coefficients according to the changed prediction step size to generate the error correction coefficient matrix. The error correction coefficient matrix and the softening matrix are subjected to Hadamard product operation, and the resulting matrix is ​​used to perform weighted correction on the illuminance deviation vector.

9. The advanced process control method for coordinated lighting of UAV light sources in disaster areas according to claim 4, characterized in that, The generalized inverse matrix is ​​obtained by performing an inversion operation on the multivariable non-minimum phase transfer function matrix in a preset frequency band, including: collecting the pitch angle change frequency and yaw angle change frequency of the UAV in the current control cycle, and determining the maximum value of the pitch angle change frequency and the yaw angle change frequency as the current main frequency disturbance frequency; Based on the current dominant frequency perturbation frequency, the adaptive forgetting factor value corresponding to the current dominant frequency perturbation frequency is determined by a lookup table method; The off-diagonal elements in the multivariate non-minimum phase transfer function matrix are attenuated and weighted using the adaptive forgetting factor value. The inversion operation is then performed on the attenuated and weighted multivariate non-minimum phase transfer function matrix to generate the generalized inverse matrix that weakens the off-diagonal coupling strength.

10. The advanced process control method for coordinated lighting of UAV light sources in disaster areas according to claim 6, characterized in that, After algebraically adding the negative value of the compensated angular velocity component as the feedforward compensation amount to the bottom attitude angular velocity calculated by the UAV inertial measurement unit, the method further includes: determining whether the absolute value of the comprehensive angular velocity command generated after algebraic addition is less than the static friction angular velocity threshold of the light source gimbal motor. When the absolute value of the comprehensive angular velocity command is less than the static friction angular velocity threshold, a dead zone compensation voltage pulse sequence consistent with the direction sign of the comprehensive angular velocity command is output to the motor driver. The amplitude of the dead zone compensation voltage pulse sequence is determined by linear interpolation based on the difference between the static friction angular velocity threshold and the comprehensive angular velocity command, so as to drive the light source gimbal motor to break through the static friction state and perform micro-motion compensation.

Images