An AIS signal constraint least square angle estimation method, device, equipment and storage medium

Abstract

The application provides an AIS signal constraint least square angle estimation method, device and equipment and a storage medium, three kinds of TDOA observation data in time domain, frequency domain and phase domain are collected, and each observation accuracy is respectively evaluated based on signal-to-noise ratio, goodness of fit and phase stability; then, multi-source observation values are mapped into a unified nonlinear observation equation, and a weighted matrix is constructed according to the observation accuracy, so that a high-precision channel obtains greater weight; then, combined with the physical incident angle range constraint of a double / three antenna array, the range constraint of the time difference of arrival determined by the baseline length and the propagation speed, and the three-antenna closed TDOA geometric constraint, an equation and inequality constraint set is formed; finally, an optimization target is constructed in the form of a constraint least square, an optimal angle estimation value is output under the premise of meeting all physical and geometric constraints by using a sequential quadratic programming (SQP) algorithm iteration, and high-precision estimation of the AIS signal angle of arrival in a complex electromagnetic environment is realized.

Description

Technical Field

[0001] This invention relates to the field of AIS signal angle estimation, and particularly to an AIS signal constrained least squares angle estimation method, apparatus, device, and storage medium. Background Technology

[0002] In modern maritime traffic management and ship navigation systems, Automatic Identification Systems (AIS) have become an important technical means to ensure maritime safety and traffic supervision. Angle estimation of AIS signals is a key technological foundation for realizing functions such as ship positioning, track tracking, and collision warning, and its estimation accuracy directly affects the reliability and safety of maritime traffic management systems.

[0003] Traditional AIS signal angle estimation methods are mainly based on single-point estimation or simple peak detection algorithms. These methods have the following technical limitations: First, they lack statistical optimality guarantees and are not designed based on rigorous statistical theory, so they cannot guarantee the optimality of the estimation results. Second, their uncertainty quantification capabilities are insufficient, and they cannot provide reliable confidence intervals and accuracy assessments, making it difficult to meet the needs of safety-critical applications. Third, they do not fully utilize constraints and do not effectively utilize physical and geometric constraints to improve estimation performance and eliminate ill-conditioned solutions. Finally, the multi-observation information fusion methods are simple, often using simple averaging or weighting methods, and lack scientific fusion strategies that consider differences in observation quality.

[0004] In view of the above, this application is hereby submitted. Summary of the Invention

[0005] This invention discloses an angle estimation method, apparatus, device, and storage medium for constrained least squares of AIS signals, aiming to achieve high-precision estimation of the arrival angle of AIS signals in complex electromagnetic environments.

[0006] The first embodiment of the present invention provides an angle estimation method based on AIS signal-constrained least squares, comprising: Receive AIS signal time difference of arrival observation data from multiple observation domains, and perform quality assessment on each observation data to obtain the measurement accuracy of each observation data; A unified observation equation is established based on the time difference of arrival observation data, and an observation weight matrix is ​​constructed based on the measurement accuracy of each observation data. The observation equation is used to represent the functional relationship between the observed value and the angle of arrival of the signal to be estimated. Based on the physical characteristics and geometric configuration of the antenna array, constraints are established, including physical constraints that limit the physically feasible range of the estimated angle and geometric constraints that ensure the consistency of multi-antenna measurements. By combining the observation equation, the observation weight matrix, and the constraints, a constrained least squares optimization problem is constructed. The constrained optimization algorithm is used to solve the constructed optimization problem, and the angle estimate that satisfies all constraints and minimizes the objective function is obtained through iterative calculation.

[0007] Preferably, the step of establishing a unified observation equation based on the arrival time difference observation data specifically involves: The m observation values ​​are combined to form an observation vector z = [z1, z2, ..., z]. m ] T Each observation zᵢ represents the arrival time difference τ obtained from time-domain cross-correlation. i The arrival time difference τ obtained by fitting the time in the time domain or frequency domain is _time or frequency domain fitting. i The arrival time difference τ obtained from _freq or carrier phase i _phase; Establish the theoretical model function h(θ) = [h1(θ), h2(θ), ..., h m (θ)] T , where θ is the arrival angle of the signal to be estimated; For a dual-antenna system, the theoretical model function is represented as h. i (θ)=(d / c)×sin(θ), where d is the antenna baseline length and c is the signal propagation speed; For a three-antenna system, the theoretical model function is represented as h1(θ) = (d 12 / c)×sin(θ),h2(θ)=(d 13 / c)×sin(θ-α 13 h3(θ)=(d 23 / c)×sin(θ-α 23 ), where d ij Let α be the baseline length between the i-th antenna and the j-th antenna. ij This is the angle between the corresponding baseline and the reference direction; Establish the observation equation z = h(θ) + ε, where ε is the observation error vector.

[0008] Preferably, the step of constructing an observation weight matrix based on the measurement accuracy of each observation data is as follows: Quality assessments were performed on observation data from the time, frequency, and phase domains respectively: Specifically, the measurement accuracy σ of the time-domain observation data was assessed. i The measurement accuracy σ of the frequency domain observation data is determined based on the signal-to-noise ratio (SNR) and interpolation quality. i The measurement accuracy σ of the phase domain observation data is determined based on the goodness of fit and coherence. i Based on phase stability and ambiguity confidence, the inverse square of the precision of each observation data is used to construct the initial weight matrix: W=diag([1 / σ1 2 ,1 / σ22 ,...,1 / σ m 2 When there is correlation between observations, the weight matrix is ​​corrected to W=R. -1 , where R is the observation covariance matrix; During the iteration process, the weights w are adaptively adjusted using a robust weight function based on the magnitude of the residual ᵢ. i =w i0 ×robust_weight(residual i ).

[0009] Preferably, the constraints include: Angle range constraint: -90°≤θ≤90° for a dual-antenna system, and 0°≤θ<360° for a three-antenna system; Time difference of arrival range constraint: The absolute value of the time difference of arrival between any two antennas shall not exceed the ratio of the baseline length to the signal propagation speed, expressed as |τ|. ij |≤d ij / c; Three-antenna closure constraint: The time difference between the first antenna and the third antenna is equal to the time difference between the first antenna and the second antenna plus the time difference between the second antenna and the third antenna, expressed as τ. 13 =τ 12 +τ 23 ; The inequality constraint is expressed as g(θ)≤0, where g(θ)=[θ-θ max ,θ min -θ,|h(θ)|-τ max ] T The equality constraint is expressed as c(θ)=0. For a three-antenna system, c(θ)=h3(θ)-h1(θ)-h2(θ).

[0010] Preferably, the construction of the constrained least squares optimization problem includes: Define the weighted least squares objective function as f(θ) = (zh(θ)). T W(zh(θ)) represents the weighted sum of squares of the differences between the observed values ​​and the theoretical model values; Construct the Lagrangian function L(θ,λ,μ)=f(θ)+λ T c(θ)+μ T g(θ), where λ is the Lagrange multiplier vector of equality constraints and μ is the Lagrange multiplier vector of inequality constraints; Establish the optimality condition of KKT: gradient condition ∇θL=∇f+A T λ+B Tμ=0, equality constraint c(θ)=0, inequality constraint g(θ)≤0, nonnegativity constraint μ≥0, complementary relaxation condition μ T g(θ)=0, where A=∂c / ∂θ is the equality-constrained Jacobian matrix, and B=∂g / ∂θ is the inequality-constrained Jacobian matrix.

[0011] Preferably, the optimization problem constructed by employing a constrained optimization algorithm is solved iteratively to obtain an angle estimate that satisfies all constraints and minimizes the objective function, specifically as follows: In the k-th iteration, θ is estimated based on the current angle. k and its multiplier λ k μ k Constructing a quadratic subproblem: min(1 / 2)d T H k d+∇f k T d, the constraint is A k d+c(θ k )=0 and B k d+g(θ k )≤0, where d=θ k+1 -θ k H is the search direction. k Let the Lagrangian function L(θ,λ,μ) be at θ k The approximate matrix of the Hessian matrix at point ∇f k The gradient of the objective function; Update the Hessian matrix using the BFGS formula: H k+1 =H k +(y k y k T ) / (y k T s k )-(H k s k s k T H k ) / (s k T H k s k ), where s k =θ k+1 -θ k y k =∇L(θ k+1 ,λ k+1 ,μ k+1 )-∇L(θ k ,λ k ,μ k ) The step size α is determined using a line search based on the Armijo criterion. k This makes merit_function(θ) k +α k d k )≤merit_function(θ k )+c1α k ∇merit_function k T d k Where c1 is the line search parameter, and merit(·) is the comprehensive evaluation function that combines the objective function and the penalty term; The iteration terminates and the optimal estimate is output when the following convergence conditions are met simultaneously: Gradient norm ||∇L||<ε0; Parameter increment ||θ k+1 -θ k ||<ε θ ; Change of objective function |f(θ) k+1 )-f(θ k )|<ε f ; The constraint satisfies ||c(θ) k )||<ε c And max(g(θ) k ))<ε g .

[0012] Among them, ε0, ε θ ε f ε c and ε g These are the preset convergence thresholds.

[0013] A second embodiment of the present invention provides an angle estimation device for AIS signal-constrained least squares, comprising: The quality assessment unit is used to receive AIS signal time difference observation data from multiple observation domains and to perform quality assessment on each observation data to obtain the measurement accuracy of each observation data. The observation model construction unit is used to establish a unified observation equation based on the time difference of arrival observation data, and to construct an observation weight matrix based on the measurement accuracy of each observation data. The observation equation is used to represent the functional relationship between the observed value and the arrival angle of the signal to be estimated. The constraint establishment unit is used to establish constraints based on the physical characteristics and geometric configuration of the antenna array. The constraints include physical constraints that limit the physically feasible range of the estimated angle and geometric constraints that ensure the consistency of multi-antenna measurements. The optimization problem construction unit is used to construct a constrained least squares optimization problem by combining the observation equation, the observation weight matrix, and the constraints. The angle estimation calculation unit is used to solve the constructed optimization problem using a constrained optimization algorithm, and iteratively calculates the angle estimate that satisfies all constraints and minimizes the objective function.

[0014] The third embodiment of the present invention provides an AIS signal constrained least squares angle estimation device, including a memory and a processor. The memory stores a computer program, which can be executed by the processor to implement an AIS signal constrained least squares angle estimation method as described in any of the above embodiments.

[0015] The fourth embodiment of the present invention provides a computer-readable storage medium, characterized in that it stores a computer program, which can be executed by the processor of the device in which the computer-readable storage medium is located, to implement an AIS signal constrained least squares angle estimation method as described in any of the above claims.

[0016] Based on the constrained least squares angle estimation method, apparatus, device, and storage medium for AIS signals provided by this invention, three types of TDOA observation data—time domain, frequency domain, and phase domain—are collected, and the observation accuracy of each channel is evaluated based on signal-to-noise ratio, goodness of fit, and phase stability. Then, the multi-source observations are mapped to a unified nonlinear observation equation, and a weighted matrix is ​​constructed according to each observation accuracy, giving higher-precision channels greater weight. Next, the physical incident angle range constraint of dual / triple antenna arrays, the arrival time difference range constraint determined by baseline length and propagation speed, and the geometric constraint of three-antenna closed TDOA are combined to form a set of equality and inequality constraints. Finally, an optimization objective is constructed in the form of constrained least squares, and iteratively solved using a sequential quadratic programming (SQP) algorithm. Under the premise of satisfying all physical and geometric constraints, the optimal angle estimate is output, achieving high-precision estimation of the AIS signal arrival angle in complex electromagnetic environments. Attached Figure Description

[0017] Figure 1 This is a flowchart illustrating an angle estimation method based on AIS signal-constrained least squares provided in the first embodiment of the present invention. Figure 2 This is a schematic diagram of a module of an AIS signal constrained least squares angle estimation device provided in the second embodiment of the present invention. Detailed Implementation

[0018] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0019] To better understand the technical solution of the present invention, the embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0020] It should be understood that the described embodiments are merely some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0021] The terminology used in the embodiments of this invention is for the purpose of describing particular embodiments only and is not intended to limit the invention. The singular forms “a,” “the,” and “the” as used in the embodiments of this invention and the appended claims are also intended to include the plural forms unless the context clearly indicates otherwise.

[0022] It should be understood that the term "and / or" used in this article is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A existing alone, A and B existing simultaneously, and B existing alone. Additionally, the character " / " in this article generally indicates that the preceding and following related objects have an "or" relationship.

[0023] Depending on the context, the word "if" as used here can be interpreted as "when," "when," "in response to determination," or "in response to detection." Similarly, depending on the context, the phrase "if determination" or "if detection (of the stated condition or event)" can be interpreted as "when determination," "in response to determination," "when detection (of the stated condition or event)," or "in response to detection (of the stated condition or event)."

[0024] The terms "first" and "second" used in the embodiments are merely to distinguish similar objects and do not represent a specific ordering of objects. It is understood that "first" and "second" can be interchanged in a specific order or sequence where permissible. It should be understood that the objects distinguished by "first" and "second" can be interchanged where appropriate so that the embodiments described herein can be implemented in orders other than those illustrated or described herein.

[0025] The specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0026] This invention discloses an angle estimation method, apparatus, device, and storage medium for constrained least squares of AIS signals, aiming to achieve high-precision estimation of the arrival angle of AIS signals in complex electromagnetic environments.

[0027] The first embodiment of the present invention provides an AIS signal-constrained least squares angle estimation method, which can be executed by an estimation device, specifically by one or more processors within the estimation device, to at least implement the following steps: S101 receives AIS signal time difference observation data from multiple observation domains, performs quality assessment on each observation data, and obtains the measurement accuracy of each observation data. In this embodiment, the estimation device (or system) can be a workstation, desktop computer, laptop computer, server, or other terminal with data processing capabilities. The estimation device can be equipped with a corresponding operating system and application software, and the functions required in this embodiment can be realized through the combination of the operating system and application software.

[0028] In this embodiment, the estimation device is equipped with three receiving antennas, forming three baselines for measuring the time difference of arrival (TDOA) of the AIS signal. When an AIS signal is received from the ship, the estimation device simultaneously extracts TDOA observation data from three different processing domains: the time domain, the frequency domain, and the phase domain.

[0029] Specifically, the system first receives AIS signals through an antenna array. The signal received by each antenna is processed by the radio frequency front end and then digitally sampled. For time-domain processing, the system uses a cross-correlation algorithm to calculate the time delay between the signals received by different antennas, and determines the arrival time difference τ by finding the peak position of the cross-correlation function. i _time. During this process, the system records the sharpness of the cross-correlation peaks and the signal-to-noise ratio (SNR). A higher SNR and sharper cross-correlation peaks indicate better quality time-domain observations. The system comprehensively evaluates the measurement accuracy σ of the time-domain observations based on the SNR value and the accuracy of the interpolation algorithm. i _time, for example, when the SNR is 20dB, combined with the accuracy of parabolic interpolation, the measurement accuracy can be obtained as approximately 0.1 sampling periods.

[0030] For frequency domain processing, the system converts the received signal to the frequency domain using a Fast Fourier Transform (FFT) and uses the signal's phase spectrum information for time delay estimation. In the frequency domain, the cross-power spectrum phase of the signals received by the two antennas exhibits a linear relationship, with its slope corresponding to the time difference of arrival τ. i The system extracts this slope by fitting a line using least squares, and simultaneously calculates the goodness-of-fit Rfreq. 2 And spectral coherence index. When the goodness of fit R 2A value close to 1 and high coherence indicate good frequency domain observation quality. The system determines the measurement accuracy σ of the frequency domain observation based on the fitting residuals and the coherence coefficient. i _freq typically achieves an accuracy of 0.05 sampling periods under high coherence conditions.

[0031] For phase domain processing, the system extracts the carrier phase information of the AIS signal and obtains a high-precision time difference of arrival τ through phase difference measurement. i _phase. Due to the 2π ambiguity problem in phase measurements, the system needs to first perform ambiguity resolution using a coarse estimate in the time or frequency domain. Under conditions of continuous and stable phase, the phase domain provides the highest measurement accuracy. The system evaluates the measurement accuracy σ of phase domain observations by monitoring phase stability indices, such as the phase change rate and Allan variance, as well as the confidence level of the ambiguity resolution. i _phase. Under ideal conditions, the phase domain measurement accuracy can reach 0.01 sampling periods.

[0032] In practical implementation, the system organizes all observation data into an observation vector z=[τ 12 _time,τ 13 _time,τ 23 _time,τ 12 _freq,τ 13 _freq,τ 23 _freq,τ 12 _phase,τ 13 _phase,τ 23 _phase] T The subscript indicates an antenna pair, containing a total of 9 observations. Each observation has a corresponding measurement accuracy assessment result. Furthermore, the system implements a real-time monitoring mechanism for observation quality. When the quality index of a certain observation domain falls below a preset threshold, the system automatically reduces the weight of that observation or removes it, ensuring that only reliable observation data participates in the final angle estimation. For example, if the SNR of the time domain cross-correlation is less than 10dB, the system considers the observation unreliable, sets its weight to 0, or uses a robust estimation method to reduce its impact.

[0033] S102, A unified observation equation is established based on the time difference of arrival observation data, and an observation weight matrix is ​​constructed based on the measurement accuracy of each observation data, wherein the observation equation is used to represent the functional relationship between the observed value and the signal arrival angle to be estimated; It should be noted that the system organizes all time difference of arrival (TDOA) measurements from different processing domains and antenna pairs into a unified observation vector. In specific implementation, assuming the system is configured with three receiving antennas, forming three baselines (1-2, 1-3, 2-3), and each baseline provides TDOA measurements in all three processing domains (time domain, frequency domain, and phase domain), then the observation vector contains nine elements, represented as z = [z1, z2, ..., z9]. T Where z1 to z3 correspond to the time-domain arrival time differences τ of the three baselines, respectively. 12 _time,τ 13 _time,τ 23 _time, the frequency domain measurement value τ corresponding to z4 to z6 12 _freq,τ 13 _freq,τ 23 _freq, the phase domain measurement value τ corresponding to z7 to z9 12 _phase,τ 13 _phase,τ 23 _phase. Vectorized representation enables the fusion of observational information from different domains within the same mathematical framework.

[0034] To establish the mathematical relationship between the observed values ​​and the angle to be estimated, the system constructs a theoretical model function h(θ). This function, based on the physical laws of electromagnetic wave propagation, describes the theoretical arrival time difference corresponding to each antenna when the signal is incident from the direction of angle θ. In the three-antenna system of this embodiment, antenna 1 is set as the reference antenna, located at the origin, while antennas 2 and 3 are located at different positions, forming a baseline length d. 12 =10 meters, d 13 =10 meters, d 23 An equilateral triangle array with a length of 10 meters. Baseline 1-2 is along the positive x-axis, and baseline 1-3 makes an angle α with the x-axis. 13 =60°, the angle α between baseline 2-3 and the x-axis 23 =120°.

[0035] The specific form of the theoretical model function is: for the theoretical time difference h between baselines 1 and 2 i (θ)=(d / c)×sin(θ), where c=3×10 8 meters per second is the propagation speed of electromagnetic waves; for baseline 1-3, due to its 60° angle with the reference direction, the theoretical time difference h2(θ) = (d 13 / c)×sin(θ-α 13 For baseline 2-3, the theoretical time difference h3(θ) = (d 23 / c)×sin(θ-α 23Since the theoretical values ​​for each baseline are identical across the three processing domains, the complete theoretical model vector is: h(θ)=[h1(θ),h2(θ),h3(θ),h1(θ),h2(θ),h3(θ),h1(θ),h2(θ),h3(θ)] T .

[0036] The observation equation z = h(θ) + ε establishes the relationship between actual observations and the theoretical model, where ε represents the observation bias caused by various error factors, including thermal noise, multipath effects, and systematic errors. Ideally, if the true angle of incidence is known and there are no errors, the observed value should be exactly equal to the theoretical value. However, in reality, biases always exist, and it is by minimizing these biases that the most probable angle of incidence is estimated.

[0037] When constructing the observation weight matrix, the system fully utilizes the aforementioned quality assessment results. Taking a specific measurement as an example, when the time-domain measurement of baseline 1-2 shows an SNR of 25dB and good interpolation quality, the system assesses its measurement accuracy as σ1 = 0.08 sampling periods; while the frequency-domain measurement of the same baseline, due to its goodness of fit R... 2 =0.98, coherence coefficient reaches 0.95, and higher accuracy assessment σ4=0.05 sampling periods is obtained; under the condition of phase stability and correct ambiguity resolution, the accuracy of phase domain measurement can reach σ7=0.01 sampling periods.

[0038] Based on these accuracy assessments, the system constructs a diagonal weight matrix W=diag([1 / 0.08 2 ,1 / 0.10 2 ,1 / 0.09 2 ,1 / 0.05 2 ,1 / 0.06 2 ,1 / 0.05 2 ,1 / 0.01 2 ,1 / 0.015 2 ,1 / 0.012 2 This weighting ensures that high-precision phase domain observations dominate the estimation, while lower-precision time domain observations still provide auxiliary information, achieving optimized fusion of multi-domain information.

[0039] In practical processing, the system also considers potential correlations between observations. For example, when using the same signal segment for time-domain and frequency-domain processing, the errors of the two may be correlated. In this case, the system obtains the observation covariance matrix R by analyzing historical data or theoretical derivation, and corrects the weight matrix to W=R. -1 This ensures that the correlations are handled correctly.

[0040] More importantly, the system implements an adaptive weight adjustment mechanism. During the iterative solution process, after obtaining a new angle estimate each time, the system calculates the residual *residualᵢ* = *zᵢ* - *hᵢ*(θ) for each observation. When the residual of an observation is unusually large, it may indicate that the observation has been disturbed or contains gross errors. The system adopts a Huber-type robust weight function, where *residualᵢ* = *zᵢ* - *hᵢ*(θ). i |≤1.345σ i When |residual i |>1.345σ i At that time, the weight is based on w i =w i0 ×1.345σ i / |residual i |Weigh down the data. This adaptive mechanism enables the algorithm to automatically identify and suppress anomalous observations.

[0041] S103, Based on the physical characteristics and geometric configuration of the antenna array, establish constraints, including physical constraints that limit the physically feasible range of the estimated angle and geometric constraints that ensure the consistency of multi-antenna measurements. In this embodiment, the system establishes a constraint system based on the actual configuration of the antenna array and the propagation characteristics of the AIS signal. For the configured three-antenna equilateral triangular array, the physically feasible range of the signal incident angle is 0° to 360° omnidirectional coverage. This is because the three-antenna system can eliminate the inherent ambiguity problem of the two-antenna system. The system expresses this angle range constraint as two inequalities: θ-360°≤0 and 0°-θ≤0, ensuring that the estimated angle is always within the effective range.

[0042] The physical constraint of the time difference of arrival stems from a fundamental fact: electromagnetic waves propagate at a finite speed. Therefore, the time difference between any two antennas receiving the same signal cannot exceed the time required for the signal to travel the distance between the two antennas. In this embodiment, the length of each of the three baselines is 10 meters, and the speed of light c = 3 × 10⁻⁶. 8 m / s, therefore the maximum possible time difference of arrival is 10 / (3×10 m / s). 8 The time difference of arrival for each observation is approximately 33.3 nanoseconds. The system imposes a constraint on the time difference of arrival for each observation: |τ| ≈ 33.3 nanoseconds. 12 |≤33.3ns、|τ 13 |≤33.3ns、|τ 23 |≤33.3ns. These constraints apply not only to actual observations but also indirectly to the angle estimates through the theoretical model function h(θ).

[0043] The unique feature of a three-antenna system lies in its geometric closure characteristic. When a signal is incident from a certain direction, the three arrival time differences must satisfy a vector closure relationship. Specifically, the time difference from antenna 1 to antenna 3 must be equal to the time difference from antenna 1 to antenna 2 plus the time difference from antenna 2 to antenna 3, i.e., τ. 13 =τ 12 +τ 23 This relationship should hold true for measurements in the time, frequency, and phase domains. Ideally, this equation holds strictly; however, in actual measurements, slight deviations may occur due to various errors. The system incorporates this relationship as an equality constraint into the optimization framework, forcing the estimation results to satisfy this geometric consistency.

[0044] To handle these constraints uniformly within the mathematical optimization framework, the system organizes all inequality constraints into a vector function g(θ). For the angular range constraint, g1(θ) = θ - 360° and g2(θ) = 0° - θ; for the time difference of arrival constraint, since the theoretical model h(θ) calculates the theoretical time difference between each baseline, the constraint is expressed as g3(θ) = |h1(θ)| - 33.3 × 10⁻⁶. -9 g4(θ) = |h2(θ)| - 33.3 × 10 -9 g5(θ) = |h3(θ)| - 33.3 × 10 -9 The complete inequality constraint vector is g(θ) = [g1(θ), g2(θ), g3(θ), g4(θ), g5(θ)]. T The requirement is that g(θ) ≤ 0.

[0045] The equality constraint is expressed by the function c(θ). In a three-antenna system, the geometric closure constraint requires h3(θ) - h1(θ) - h2(θ) = 0. Expanded, this becomes: (d 23 / c)×sin(θ-120°)-(d 12 / c)×sin(θ)-(d 13 / c)×sin(θ-60°)=0. Since the three baselines are of equal length in this embodiment, the constraint simplifies to: sin(θ-120°)-sin(θ)-sin(θ-60°)=0. This trigonometric equation ensures the geometric consistency between the estimated angle and the three time difference measurements.

[0046] In practical implementation, the system also considers the accuracy of numerical calculations. Due to rounding errors in floating-point operations, the strict equality constraint c(θ)=0 may never be precisely satisfied numerically. Therefore, the system sets a very small tolerance ε_c=10. -8 An equality constraint is considered satisfied when |c(θ)| < ε_c. Similarly, for inequality constraints, when g_i(θ) < ε_g = 10...-8 The boundary where the constraint is considered active is then considered to be in a state of equilibrium. It should be noted that the introduction of constraints improves the reliability of angle estimation. For example, when an observation is severely disturbed, leading to abnormal measurements, without constraints, least squares estimation might yield physically impossible results, such as angles exceeding 360° or time differences exceeding physical limits. With these constraints, the optimization algorithm automatically seeks solutions that both fit most normal observations well and satisfy all physical and geometric constraints, thus providing more robust estimation results.

[0047] S104. Combining the observation equation, the observation weight matrix, and the constraints, a constrained least squares optimization problem is constructed. It should be noted that after establishing the observation equations, constructing the weight matrix, and setting the constraints, all the aforementioned elements are integrated into a unified mathematical optimization framework to form a solvable constrained least squares problem. Specifically, a weighted least squares objective function f(θ) is first constructed. This function quantifies the weighted deviation between the actual observed value and the theoretical predicted value at a given angle θ. In practice, the system calculates the observation residual vector zh(θ), where z contains the arrival time differences of nine actual measurements, and h(θ) is the corresponding theoretical calculation value. For example, when θ = 45°, the theoretical time difference between baselines 1 and 2 is h1(45°) = (10 / 3 × 10) / 2. 8 )×sin(45°)≈23.6 nanoseconds. If the actual measured value z1=24.1 nanoseconds, then the residual of this observation is 0.5 nanoseconds.

[0048] Objective function f(θ)=(zh(θ)) T W(zh(θ)) performs a weighted sum of squares on each residual using a weight matrix W. In this embodiment, assuming the weight of phase domain observation z7 is w7=10000 (corresponding to high precision σ7=0.01), and the weight of time domain observation z1 is w1=156.25 (corresponding to lower precision σ1=0.08), then the contribution of the 0.1 nanosecond phase domain residual to the objective function is comparable to the contribution of the 0.8 nanosecond time domain residual. This weighting mechanism ensures that high-precision observations play a dominant role in angle estimation.

[0049] However, simply minimizing the objective function f(θ) can lead to unreasonable results. For example, if all observations are affected by systematic biases in the same direction, unconstrained optimization may yield angle estimates that are beyond physical limits. Therefore, the system incorporates constraints into the optimization framework using the Lagrangian method. The Lagrangian function L(θ,λ,μ)=f(θ)+λ T c(θ)+μ T g(θ) cleverly transforms the constrained optimization problem into finding the stationary point of the function.

[0050] In this function, λ is the Lagrange multiplier associated with the equality constraint c(θ). For the geometric closure constraint c(θ) = h3(θ) - h1(θ) - h2(θ) = 0 for a three-antenna system, the corresponding multiplier λ reflects the "influence strength" of the constraint. When the observation data itself satisfies the closure relationship well, the value of λ is small; conversely, when there are inconsistencies in the observations, the optimization process will automatically adjust the value of λ to find a balance between satisfying the constraint and fitting the observations.

[0051] Similarly, the vector μ = [μ1, μ2, μ3, μ4, μ5] T It includes Lagrange multipliers corresponding to the five inequality constraints. These multipliers have a special property: when a certain constraint is not activated (e.g., g3(θ) = |h1(θ)| - 33.3 × 10⁻⁶), they are activated. -9 <0, meaning the time difference is much smaller than the physical limit), corresponding to μ3=0; only when the constraint approaches or reaches the boundary will the corresponding multiplier take a positive value, affecting the optimization process.

[0052] The established KKT (Karush-Kuhn-Tucker) optimality conditions provide the necessary and sufficient conditions for determining whether a solution is optimal for a constrained optimization problem. The gradient condition is ∇θL = ∇f + A. T λ+B T μ=0 requires that at the optimum, the gradient of the objective function and the gradient of the constraint function are in balance. Where ∇f=-2J T W(zh(θ)), where J is the Jacobian matrix of the observation function. In this embodiment, the first row of J is ∂h1 / ∂θ=(10 / 3×10 8 )×cos(θ) reflects the effect of angle change on the time difference between baselines 1 and 2.

[0053] The equality constraint Jacobian matrix A = ∂c / ∂θ describes how the closed constraint changes with angle. Calculated using the chain rule, A = ∂h³ / ∂θ - ∂h¹ / ∂θ - ∂h² / ∂θ = (10 / 3 × 10⁻¹⁰) 8 )×[cos(θ-120°)-cos(θ)-cos(θ-60°)]. Similarly, the inequality constraint Jacobian matrix B contains the gradient information of each inequality constraint.

[0054] Complementary relaxation condition μ T g(θ)=0 embodies the ingenuity of optimization theory: for each inequality constraint, either the constraint is inactive (g_i(θ)<0) and the corresponding multiplier is zero (μᵢ=0), or the constraint is on the boundary (g_i(θ)=0) and the multiplier may be positive. This ensures that only the constraints that actually take effect will influence the final solution.

[0055] S105 uses a constrained optimization algorithm to solve the constructed optimization problem, and iteratively calculates the angle estimate that satisfies all constraints and minimizes the objective function.

[0056] In this embodiment, it is assumed that the system is in the 5th iteration, the current angle estimate θ5 = 32.46°, the equality constraint multiplier λ5 = -0.015, and the inequality constraint multiplier vector μ5 = [0,0,0.002,0,0]ᵀ (indicating that the third inequality constraint is close to activation). Based on these current values, the system constructs a quadratic programming subproblem to determine the next search direction.

[0057] The construction process of the quadratic subproblem embodies the core idea of ​​sequential quadratic programming: approximating the original nonlinear objective function at the current point with a quadratic function, and approximating the nonlinear constraints with linear constraints. Specifically, the system calculates the gradient of the objective function at the current point, ∇f5=-2J5. T W(zh(θ5)), where J5 is the Jacobian matrix of the observation function at θ5. Assuming the principal component of the weighted residual vector W(zh(θ5)) comes from phase domain observations, the gradient value is calculated as ∇f5 = -3.2 × 10⁻⁶. -7 .

[0058] The Hessian matrix H5 approximates the second derivative information of the Lagrange function, determining the curvature characteristics of the quadratic subproblem. The system uses the BFGS method to maintain and update this matrix. From the 4th to the 5th iteration, the parameter changes s4 = θ5 - θ4 = 32.46° - 32.38° = 0.08°, and the gradient of the Lagrange function changes from ∇L4 = -2.8 × 10⁻⁶. -7 The value changes to ∇L5 = -3.5 × 10. -7 Therefore, y4 = -0.7 × 10 -7 Applying the BFGS update formula, the system calculates the rank-two correction term: the first term (y4y4) T ) / (y4 T s T Add curvature information along the y4 direction, the second term -(H4s4s4) T H4) / (s4 T H4s4) removes outdated curvature information, thus obtaining a more accurate approximation of the Hessian matrix H5.

[0059] Linearize the constraints at the current point. The equality constraint c(θ5) = h3(θ5) - h1(θ5) - h2(θ5) = 2.1 × 10 -10 (Minor violation), its Jacobian A5 = -0.85 × 10 -8 In the inequality constraints, since μ3 > 0, it indicates that g3(θ5) = |h1(θ5)| - 33.3 × 10 -9Approaching zero. The linearized constraints require the search direction d to satisfy A5d+c(θ5)=0 (to make the equality constraints tend to be satisfied) and B5d+g(θ5)≤0 (to keep the inequality constraints satisfied).

[0060] Solve the quadratic programming problem: min(1 / 2)dᵀH5d+∇f5 T The system uses the interior point method or the effective set method to obtain the search direction d5 = 0.18°. This direction comprehensively considers multiple requirements, including objective function descent, equality constraint satisfaction, and inequality constraint preservation.

[0061] However, traversing the full step length along the d5 direction may lead to constraint violations or an increase in the objective function. The system determines the appropriate step size α5 through line search. The value function merit(θ) = f(θ) + ρ||c(θ)||1 + ρΣmax(0, g i The quality of the solution is comprehensively evaluated using (θ) as a parameter, where the penalty parameter ρ = 150 ensures that constraint violations are appropriately penalized. The system checks the value function values ​​under different time lengths: when α = 1, merit(θ5 + d5) = 0.0125; when α = 0.5, merit(θ5 + 0.5d5) = 0.0098; and when α = 0.25, merit(θ5 + 0.25d5) = 0.0089. Since merit(θ5) = 0.0092, α5 = 0.5 is chosen to satisfy the Armijo criterion: merit(θ5 + 0.5d5) ≤ merit(θ5) + 10. -4 ×0.5×(-0.012).

[0062] The updated value is θ6 = θ5 + 0.5 × 0.18° = 32.55°, and the system continues to iterate. After several iterations, in the 12th iteration, the system checks the convergence condition: the Lagrange function gradient norm ||∇L 12 ||=7.3×10 -7 <10 -6 ; Parameter change ||θ 13 -θ 12 ||=5.6×10 -9 <10 -8 The objective function remains essentially unchanged |f(θ) 13 )-f(θ 12 )|=8.9×10 -11 <10 -10 Equality constraint ||c(θ) 13 )||=6.2×10 -9 <10 -8 ; All inequality constraints max(g(θ) 13 ))=-0.008<10 -8All values ​​are negative. When all convergence conditions are met simultaneously, the system is considered convergent and outputs the optimal angle estimate θ` = 32.5847°.

[0063] In one possible implementation of the present invention: After obtaining the optimal angle estimate θ' and the corresponding Jacobian matrix J and weight matrix W, the first step is to use F=J T WJ-∑λ i ∇ 2 c i -∑μ j ∇ 2 g j Calculate the Fisher information matrix F, where J is the partial derivative matrix of the observation equation with respect to θ, W is the observation weight matrix, and λ i μ j The equality constraints c are respectively i and inequality constraints g j Lagrange multipliers, ∇ 2 c i 、∇ 2 g j The second derivative matrix is ​​F; then, the inverse of F is taken to obtain the parametric covariance matrix C_θ=F. -1 The standard error σ_θ=√diag(C_θ) is extracted. Under the assumption of normal distribution, the confidence interval (1–α) is calculated as [θ`–z_{α / 2}·σ_θ,θ`+z_{α / 2}·σ_θ] using the critical value z_{α / 2} at the significance level α. If it is necessary to evaluate the uncertainty of the angle function f(θ), its standard error is calculated according to the error propagation formula σ_f=|∂f / ∂θ|·σ_θ.

[0064] In one possible implementation of the present invention: When performing statistical analysis on the optimal angle estimation result θ', the residual vector z–h(θ') is first standardized using the square root of the observation weight matrix W, resulting in the standardized residual r=W. 1 / 2 (z–h(θ`)) should theoretically approximately follow a standard normal distribution; subsequently, the chi-square statistic χ² is calculated. 2 =(z–h(θ`))ᵀW(z–h(θ`)) and perform a chi-square test with degrees of freedom df=m–p–q (where m, p, and q represent the number of observations, the number of parameters to be estimated, and the number of equality constraints, respectively), which can test the goodness of fit of the model to the overall observations. To identify the influence of individual observation points on the estimation results, this method further constructs a hat matrix H=J(J T WJ) -1 J T W, whose diagonal element h iiThe leverage effect of the i-th observation on the overall fit is measured and based on the standardized residual r. i and h ii Calculate Cook distance Di=(r i 2 / p)·\[h ii / (1–h ii ) 2 ], when D i When the threshold of 4 / m is exceeded, the observation is determined to be an abnormal influence point, so as to prompt the subsequent removal or recalibration of the point.

[0065] In one possible implementation of the present invention: In practical implementation, this method further robustly handles the residuals and matrix numerical stability: after each iteration, the standardized residual r=W 1 / 2 After calculating (z–h(θ`) , residuals exceeding the 3σ range are first marked as outliers, and then the Huber loss function is applied. ρ(r)={r 2 / 2,|r|≤k;k|r|-k 2 / 2,|r|>k} (Where the adjustment parameter k is set to 1.345) Instead of the traditional quadratic loss, the residuals are reweighted through the influence function ψ, thereby updating the weights during the iterative reweighted least squares (IRLS) process. w i (k+1) =ψ(r i (k) / σ) / (r i ( k) / σ) If, in any iteration, the condition number of the Jacobian matrix, cond(JᵀWJ), exceeds 10. 12 Then, by introducing the Tikhonov regularization term λI (λ being the regularization parameter) into the normal equation, we construct... (J T WJ+λI)θ=J T Wz To ensure numerical stability; if the preset convergence condition is not met after more than 100 iterations, the problem automatically switches to a trust region subproblem. min f(θ k ) + ∇f k T d + (1 / 2)d T H k d By dynamically adjusting the trust region radius ||d|| ≤ Δ k To ensure algorithm convergence, we can ultimately obtain an accurate angle estimate that is both robust to outlier interference and numerically stable.

[0066] A second embodiment of the present invention provides an angle estimation device for AIS signal-constrained least squares, comprising: The quality assessment unit is used to receive AIS signal time difference observation data from multiple observation domains and to perform quality assessment on each observation data to obtain the measurement accuracy of each observation data. The observation model construction unit is used to establish a unified observation equation based on the time difference of arrival observation data, and to construct an observation weight matrix based on the measurement accuracy of each observation data. The observation equation is used to represent the functional relationship between the observed value and the arrival angle of the signal to be estimated. The constraint establishment unit is used to establish constraints based on the physical characteristics and geometric configuration of the antenna array. The constraints include physical constraints that limit the physically feasible range of the estimated angle and geometric constraints that ensure the consistency of multi-antenna measurements. The optimization problem construction unit is used to construct a constrained least squares optimization problem by combining the observation equation, the observation weight matrix, and the constraints. The angle estimation calculation unit is used to solve the constructed optimization problem using a constrained optimization algorithm, and iteratively calculates the angle estimate that satisfies all constraints and minimizes the objective function.

[0067] The third embodiment of the present invention provides an AIS signal constrained least squares angle estimation device, including a memory and a processor. The memory stores a computer program, which can be executed by the processor to implement an AIS signal constrained least squares angle estimation method as described in any of the above embodiments.

[0068] The fourth embodiment of the present invention provides a computer-readable storage medium, characterized in that it stores a computer program, which can be executed by the processor of the device in which the computer-readable storage medium is located, to implement an AIS signal constrained least squares angle estimation method as described in any of the above claims.

[0069] Based on the constrained least squares angle estimation method, apparatus, device, and storage medium for AIS signals provided by this invention, three types of TDOA observation data—time domain, frequency domain, and phase domain—are collected, and the observation accuracy of each channel is evaluated based on signal-to-noise ratio, goodness of fit, and phase stability. Then, the multi-source observations are mapped to a unified nonlinear observation equation, and a weighted matrix is ​​constructed according to each observation accuracy, giving higher-precision channels greater weight. Next, the physical incident angle range constraint of dual / triple antenna arrays, the arrival time difference range constraint determined by baseline length and propagation speed, and the geometric constraint of three-antenna closed TDOA are combined to form a set of equality and inequality constraints. Finally, an optimization objective is constructed in the form of constrained least squares, and iteratively solved using a sequential quadratic programming (SQP) algorithm. Under the premise of satisfying all physical and geometric constraints, the optimal angle estimate is output, achieving high-precision estimation of the AIS signal arrival angle in complex electromagnetic environments.

[0070] Exemplary examples show that the computer program described in the third and fourth embodiments of the present invention can be divided into one or more modules, which are stored in the memory and executed by the processor to complete the present invention. The one or more modules can be a series of computer program instruction segments capable of performing specific functions, which describe the execution process of the computer program in implementing an AIS signal-constrained least squares angle estimation device. For example, the apparatus described in the second embodiment of the present invention.

[0071] The processor referred to can be a Central Processing Unit (CPU), or other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The general-purpose processor can be a microprocessor or any conventional processor. The processor is the control center of the AIS signal-constrained least squares angle estimation method, connecting the various parts implementing the AIS signal-constrained least squares angle estimation method through various interfaces and lines.

[0072] The memory can be used to store the computer program and / or modules. The processor implements various functions of an AIS signal-constrained least squares angle estimation method by running or executing the computer program and / or modules stored in the memory, and by calling the data stored in the memory. The memory may mainly include a program storage area and a data storage area. The program storage area may store the operating system, at least one application program required for a function (such as sound playback function, text conversion function, etc.), etc.; the data storage area may store data created according to the use of the mobile phone (such as audio data, text message data, etc.). In addition, the memory may include high-speed random access memory, and may also include non-volatile memory, such as hard disk, memory, plug-in hard disk, smart media card (SMC), secure digital card (SD card), flash card, at least one disk storage device, flash memory device, or other volatile solid-state storage device.

[0073] If the implemented module is implemented as a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, all or part of the processes in the above embodiments of the present invention can also be implemented by a computer program instructing related hardware. The computer program can be stored in a computer-readable storage medium, and when executed by a processor, it can implement the steps of the various method embodiments described above. The computer program includes computer program code, which can be in the form of source code, object code, executable files, or certain intermediate forms. The computer-readable medium can include: any entity or device capable of carrying the computer program code, recording media, USB flash drives, portable hard drives, magnetic disks, optical disks, computer memory, read-only memory (ROM), random access memory (RAM), electrical carrier signals, telecommunication signals, and software distribution media, etc. It should be noted that the content included in the computer-readable medium can be appropriately added or removed according to the requirements of legislation and patent practice in the jurisdiction. For example, in some jurisdictions, according to legislation and patent practice, computer-readable media do not include electrical carrier signals and telecommunication signals.

[0074] It should be noted that the device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs. Furthermore, in the accompanying drawings of the device embodiments provided by this invention, the connection relationships between modules indicate that they have communication connections, which can be specifically implemented as one or more communication buses or signal lines. Those skilled in the art can understand and implement this without any creative effort.

[0075] The above description is merely a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for angle estimation using AIS signal-constrained least squares, characterized in that, include: Receive AIS signal time difference of arrival observation data from multiple observation domains, and perform quality assessment on each observation data to obtain the measurement accuracy of each observation data; A unified observation equation is established based on the time difference of arrival observation data, and an observation weight matrix is ​​constructed according to the measurement accuracy of each observation data. Specifically, the quality of the observation data from the time domain, frequency domain, and phase domain is evaluated respectively: where the measurement accuracy σ of the time domain observation data is... i _time is determined based on the signal-to-noise ratio (SNR) and interpolation quality, and the measurement accuracy σ of the frequency domain observation data is also considered. i _freq is determined based on goodness of fit and coherence, and σ represents the measurement accuracy of the phase domain observation data. i The phase is determined based on phase stability and ambiguity confidence. The inverse square of the precision of each observation data is used to construct the initial weight matrix: W0=diag([1 / σ1 2 ,1 / σ2 2 ,...,1 / σ m 2 When there is correlation between observations, the weight matrix is ​​corrected to W=R. -1 Where R is the observation covariance matrix; during the iteration process, based on the residual... i The size of the weight is adaptively adjusted by a robust weight function. i =w i0 ×robust_weight(residual i ), wherein the observation equation is used to represent the functional relationship between the observed value and the angle of arrival of the signal to be estimated; Based on the physical characteristics and geometric configuration of the antenna array, constraints are established, including physical constraints that limit the physically feasible range of the estimated angle and geometric constraints that ensure the consistency of multi-antenna measurements. By combining the observation equation, the observation weight matrix, and the constraints, a constrained least squares optimization problem is constructed. The constrained optimization algorithm is used to solve the constructed optimization problem, and the angle estimate that satisfies all constraints and minimizes the objective function is obtained through iterative calculation.

2. The angle estimation method for AIS signal constrained least squares according to claim 1, characterized in that, The establishment of a unified observation equation based on the arrival time difference observation data is specifically as follows: The m observation values ​​are combined to form an observation vector z = [z1, z2, ..., z]. m ] T , where each observation z i The arrival time difference τ obtained by cross-correlation in the time domain i The arrival time difference τ obtained by fitting the time in the _time domain or the frequency domain is... i The arrival time difference τ obtained from _freq or carrier phase i _phase; Establish the theoretical model function h(θ) = [h1(θ), h2(θ), ..., h m (θ)] T , where θ is the arrival angle of the signal to be estimated; For a dual-antenna system, the theoretical model function is represented as h. i (θ)=(d / c)×sin(θ), where d is the antenna baseline length and c is the signal propagation speed; For a three-antenna system, the theoretical model function is represented as h1(θ) = (d 12 / c)×sin(θ),h2(θ)=(d 13 / c)×sin(θ-α 13 h3(θ)=(d 23 / c)×sin(θ-α 23 ), where d ij Let α be the baseline length between the i-th antenna and the j-th antenna. ij This is the angle between the corresponding baseline and the reference direction; Establish the observation equation z = h(θ) + ε, where ε is the observation error vector.

3. The angle estimation method for AIS signal constrained least squares according to claim 1, characterized in that, The constraints include: Angle range constraint: -90°≤θ≤90° for a dual-antenna system, and 0°≤θ<360° for a three-antenna system; Time difference of arrival range constraint: The absolute value of the time difference of arrival between any two antennas shall not exceed the ratio of the baseline length to the signal propagation speed, expressed as |τ|. ij |≤d ij / c; Three-antenna closure constraint: The time difference between the first antenna and the third antenna is equal to the time difference between the first antenna and the second antenna plus the time difference between the second antenna and the third antenna, expressed as τ. 13 =τ 12 +τ 23 ; The inequality constraint is expressed as g(θ)≤0, where g(θ)=[θ-θ max ,θ min -θ,|h(θ)|-τ max ] T The equality constraint is expressed as φ(θ)=0. For a three-antenna system, φ(θ)=h3(θ)-h1(θ)-h2(θ).

4. The angle estimation method for AIS signal constrained least squares according to claim 1, characterized in that, The constrained least squares optimization problem includes: Define the weighted least squares objective function as f(θ) = (zh(θ)). T W(zh(θ)) represents the weighted sum of squares of the differences between the observed values ​​and the theoretical model values; Construct the Lagrangian function L(θ,λ,μ)=f(θ)+λ T φ(θ)+μ T g(θ), where λ is the Lagrange multiplier vector of equality constraints and μ is the Lagrange multiplier vector of inequality constraints; Establish the optimality condition of KKT: gradient condition θL= f+A T λ+B T μ=0, equality constraint φ(θ)=0, inequality constraint g(θ)≤0, nonnegativity constraint μ≥0, complementary relaxation condition μ T g(θ) = 0, where The Jacobian matrix is ​​an equality constraint matrix. The Jacobian matrix is ​​used to constrain inequalities.

5. The angle estimation method for AIS signal constrained least squares according to claim 4, characterized in that, The optimization problem constructed by the constrained optimization algorithm is solved iteratively to obtain the angle estimate that satisfies all constraints and minimizes the objective function, specifically as follows: In the k-th iteration, θ is estimated based on the current angle. k and its multiplier λ k μ k Constructing a quadratic subproblem: min(1 / 2)d T H k d+ f k T d, the constraint is A k d+φ(θ k )=0 and B k d+g(θ k )≤0, where p=θ k+1 -θ k H is the search direction. k Let the Lagrangian function L(θ,λ,μ) be at θ k The approximate matrix of the Hessian matrix at that location, f k The gradient of the objective function; Update the Hessian matrix using the BFGS formula: H k+1 =H k +(y k y k T ) / (y k T s k )-(H k s k s k T H k ) / (s k T H k s k ), among them k =θ k+1 -θ k ,y k = L(θ k+1 ,l k+1 ,m k+1 )- L(θ k ,l k ,m k ) The step size α is determined using a line search based on the Armijo criterion. k This makes merit_function(θ) k +α k p k )≤merit_function(θ k )+c1α k merit function k T p k Where c1 is the line search parameter, and merit(·) is the comprehensive evaluation function that combines the objective function and the penalty term; The iteration terminates and the optimal estimate is output when the following convergence conditions are met simultaneously: Gradient norm || L||<ε0; Parameter increment ||θ k+1 -θ k ||<ε θ ; Change of objective function |f(θ) k+1 )-f(θ k )|<ε f ; The constraint satisfies ||φ(θ) k )||<ε c And max(g(θ) k ))<ε g ; Where, ε0, ε θ ε f ε c and ε g These are the preset convergence thresholds.

6. An angle estimation device based on AIS signal-constrained least squares, characterized in that, include: The quality assessment unit is used to receive AIS signal time difference observation data from multiple observation domains, and to perform quality assessment on each observation data to obtain the measurement accuracy of each observation data. The observation model construction unit is used to establish a unified observation equation based on the time difference of arrival observation data, and to construct an observation weight matrix based on the measurement accuracy of each observation data. Specifically, it performs quality assessments on the observation data from the time domain, frequency domain, and phase domain respectively: where the measurement accuracy σ of the time domain observation data is... i _time is determined based on the signal-to-noise ratio (SNR) and interpolation quality, and the measurement accuracy σ of the frequency domain observation data is also considered. i _freq is determined based on goodness of fit and coherence, and σ represents the measurement accuracy of the phase domain observation data. i The phase is determined based on phase stability and ambiguity confidence. The inverse square of the precision of each observation data is used to construct the initial weight matrix: W0=diag([1 / σ1 2 ,1 / σ2 2 ,...,1 / σ m 2 When there is correlation between observations, the weight matrix is ​​corrected to W=R. -1 Where R is the observation covariance matrix; during the iteration process, based on the residual... i The size of the weight is adaptively adjusted by a robust weight function. i =w i0 ×robust_weight(residual i ), wherein the observation equation is used to represent the functional relationship between the observed value and the angle of arrival of the signal to be estimated; The constraint establishment unit is used to establish constraints based on the physical characteristics and geometric configuration of the antenna array. The constraints include physical constraints that limit the physically feasible range of the estimated angle and geometric constraints that ensure the consistency of multi-antenna measurements. The optimization problem construction unit is used to construct a constrained least squares optimization problem by combining the observation equation, the observation weight matrix, and the constraints. The angle estimation calculation unit is used to solve the constructed optimization problem using a constrained optimization algorithm, and iteratively calculates the angle estimate that satisfies all constraints and minimizes the objective function.

7. An AIS signal-constrained least squares angle estimation device, characterized in that, It includes a memory and a processor, wherein the memory stores a computer program that can be executed by the processor to implement an AIS signal-constrained least squares angle estimation method as described in any one of claims 1 to 5.

8. A computer-readable storage medium, characterized in that, The device contains a computer program that can be executed by a processor of the device on which the computer-readable storage medium is located to implement an AIS signal-constrained least squares angle estimation method as described in any one of claims 1 to 5.

Images