Method for generating satellite communication beams

Abstract

The invention discloses a method for generating a satellite communication beam. The method comprises the following steps: acquiring boundary coordinates of a target geographic corridor and converting the boundary coordinates into satellite angular domain representation; extracting a skeleton curve of a corridor based on the angular domain boundary, and establishing a parameterized beam pattern model distributed along the skeleton curve; constructing a beam forming weight optimization problem which comprises a coverage gain constraint distributed along the skeleton curve and a continuity constraint between adjacent sampling positions on the skeleton curve; and solving the optimization problem to obtain a complex weight vector, and loading the complex weight vector to the beam forming network to drive the phased-array antenna. According to the invention, dimension reduction is carried out on a two-dimensional area coverage problem through skeleton curve modeling, and stable energy transition of a wave beam along a curved corridor is ensured by utilizing a continuous constraint explicit mode, so that the problems of difficult optimization convergence and discontinuous coverage in long and narrow non-convex area coverage are effectively solved.

Description

Technical Field

[0001] This invention belongs to the field of satellite communication, and in particular to a method for generating satellite communication beams. Background Technology

[0002] In satellite communication systems, multi-beam phased array antennas can flexibly change the beam direction and shape by controlling the amplitude and phase of the radiating elements, achieving precise coverage of specific geographical areas. For large-span, irregularly shaped strip-shaped geographical areas (such as winding rivers, coastlines, or transportation routes), generating shaped beams that highly match the geographical contours is of significant technical importance for improving the utilization rate of satellite transmission power, enhancing signal strength in edge areas, and improving the efficiency of spectrum resource utilization.

[0003] Currently, beamforming techniques for irregular regions commonly employ optimization methods based on two-dimensional grid sampling. These methods discretize the target coverage area into a dense set of two-dimensional latitude and longitude grid points, treating the gain requirement of each grid point as an independent constraint. They then utilize semi-definite relaxation (SDR) or heuristic algorithms (such as genetic algorithms) to solve for the antenna weight vector, satisfying the gain threshold for all grid points and suppressing sidelobes. This method is widely used when dealing with relatively regular, well-convex clustered regions (such as the outlines of provinces or countries).

[0004] However, for horseshoe-shaped or long, narrow, and curved strip regions, existing two-dimensional sampling methods suffer from the technical bottleneck of the curse of constraint dimensionality and the difficulty in simultaneously achieving coverage continuity. Specifically, to approximate the boundary of a long, narrow, and curved region, the two-dimensional mesh method requires an extremely high density of sampling points, leading to an exponential increase in the number of constraints in the optimization problem, excessive computational complexity, and a tendency to get trapped in local optima. In addition, since power gain constraints are inherently non-convex, and traditional methods lack utilization of the region's geometric and topological features, when dealing with highly non-convex curved sections, the constraints between discrete points lack intrinsic correlation. The algorithm struggles to automatically ensure a smooth transition of beam energy along the corridor direction, easily resulting in gain collapse or coverage breaks at the bends, making it difficult to meet the continuous communication requirements of high-speed mobile terminals. Summary of the Invention

[0005] The purpose of this invention is to provide a method for generating satellite communication beams to solve the aforementioned problems existing in the prior art.

[0006] Technical solution: A method for generating satellite communication beams, comprising:

[0007] Obtain the boundary coordinates of the target geographic corridor and convert them into an angular domain representation in a satellite coordinate system;

[0008] Based on the angular domain representation, the skeleton curve of the target geographic corridor is extracted, and a parameterized beam pattern model distributed along the skeleton curve is established.

[0009] Based on the parametric beam pattern model, a beamforming weight optimization problem is constructed, which includes coverage gain constraints distributed along the skeleton curve and continuity constraints between adjacent sampling positions located on the skeleton curve.

[0010] Solve the beamforming weight optimization problem to obtain the beamforming complex weight vector;

[0011] The beamforming complex weight vector is loaded into a pre-configured beamforming network to drive the phased array antenna to radiate a communication beam that matches the shape of the target geographic corridor.

[0012] Beneficial effects: This invention reduces the dimensionality of the two-dimensional region coverage problem by using skeleton curve modeling, and uses continuity constraints to explicitly ensure the smooth energy transition of the beam along the curved corridor, effectively solving the problems of difficult optimization convergence and coverage discontinuity in the coverage of narrow and non-convex regions. Attached Figure Description

[0013] Figure 1 This is a flowchart illustrating the steps of a method for generating a satellite communication beam, provided as an embodiment of this application.

[0014] Figure 2 A flowchart illustrating the steps for constructing a corridor width function as provided in this application embodiment.

[0015] Figure 3 A flowchart illustrating the steps for generating a set of overlay slice points provided in this application embodiment.

[0016] Figure 4 A flowchart illustrating the steps of dividing a target geographical corridor into predetermined sub-segments, provided for embodiments of this application. Detailed Implementation

[0017] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. 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 should fall within the scope of protection of the present invention.

[0018] It should be noted that the terms include and have, and any variations thereof, are intended to cover non-exclusive inclusion. For example, a process, method, system, product, or device that includes a series of steps or units is not necessarily limited to those steps or units that are explicitly listed, but may include other steps or units that are not explicitly listed or that are inherent to such process, method, product, or device.

[0019] According to one aspect of this application, a system for generating a satellite communication beam includes:

[0020] The beam control and management subsystem is configured to store the boundary coordinates of the target geographic corridor and executes a method for generating satellite communication beams to calculate the beamforming complex weight vector.

[0021] A digital beamforming network, connected to a beam control and management subsystem, is configured to receive beamforming complex weight vectors and generate multiple excitation signals accordingly.

[0022] A multi-beam phased array antenna, connected to a digital beamforming network, contains multiple radiating elements and is configured to radiate a communication beam that matches the shape of the target geographic corridor and provides continuous coverage along the route, driven by multiple excitation signals.

[0023] Specifically, the physical architecture of a system used to generate satellite communication beams mainly includes a beam control and management subsystem, a digital beamforming network (DBFN), and a multi-beam phased array antenna. The beam control and management subsystem is the core computing unit of the entire system. It can be deployed at a ground control center, communicating with the satellite via an uplink telemetry and control link; or it can be directly integrated into the satellite's onboard processing payload to achieve fully autonomous on-orbit beam reconfiguration. The beam control and management subsystem typically integrates high-performance computing modules (such as field-programmable gate arrays (FPGAs) or digital signal processors (DSPs)) to execute complex beam optimization algorithms. The digital beamforming network connects to the beam control and management subsystem and is responsible for converting the calculated complex weight vectors into specific amplitude and phase control signals, which are then distributed to each antenna element. The multi-beam phased array antenna consists of N radiating elements (e.g., a 32×32 rectangular planar array), each element possessing independent amplitude and phase control capabilities. Driven by the weights generated by the beam control and management subsystem, the antenna is able to synthesize a high-gain beam with a specific shape (such as a horseshoe) in space.

[0024] like Figure 1 As shown, a method for generating a satellite communication beam includes the following steps:

[0025] Obtain the boundary coordinates of the target geographic corridor and convert them into an angular domain representation in the satellite coordinate system.

[0026] In this embodiment, the target geographic corridor refers to a continuous area on the Earth's surface with significantly elongated and curved characteristics, such as an economic corridor spanning multiple countries, a winding coastal waterway, or a specific large-scale infrastructure belt. The process of obtaining boundary coordinates can be accomplished by reading pre-stored Geographic Information System (GIS) database files, which contain a set of latitude and longitude points describing the corridor's outline. i lati After obtaining latitude and longitude data, the system needs to combine the satellite's current orbital position (such as longitude, latitude, and altitude) and attitude information (such as roll, pitch, and yaw angles) with a coordinate transformation matrix to map these geographic coordinates into a spherical coordinate system under the satellite's body coordinate system, thus obtaining the angular domain representation. The angular domain representation typically consists of two dimensions: azimuth angle θ and elevation angle Φ. Since the antenna's beam pointing and shape are defined based on the angular domain from the satellite's perspective, rather than the latitude and longitude on the ground, this transformation is necessary. The transformed angular domain boundary data will serve as the input basis for all subsequent modeling and optimization.

[0027] Specifically, let the satellite's position in the geocentric inertial coordinate system be (X... s Y s Z s Given that the latitude and longitude of the target geographic point are (lon, lat), then the position of this point in the geocentric coordinate system (X... g Y g Z g This can be obtained through the conversion from standard geographic coordinates to geocentric coordinates. Furthermore, the unit vector pointing from the satellite to this point can be expressed as:

[0028] u=[(X g -X s ), (Y g -Y s ), (Z g -Z s )] / ||(X g -X s ), (Y g -Y s ), (Z g -Z s )||;

[0029] In the satellite's body coordinate system, the azimuth angle θ and elevation angle Φ of this vector can be calculated using standard spherical coordinate transformation formulas:

[0030] θ=arctan(u y / u x ); Φ=arcsin(u z );

[0031] Where u x u y u z The components of the unit vector u on the three axes of the satellite body coordinate system.

[0032] Based on angular domain representation, the skeleton curve of the target geographic corridor is extracted, and a parametric beam pattern model distributed along the skeleton curve is established.

[0033] In traditional beamforming methods, irregular regions are typically treated as two-dimensional surfaces for mesh sampling. However, when dealing with highly non-convex regions such as horseshoe or serpentine shapes, problems often arise such as an explosion in the number of constraints and difficulty in achieving convergence. This embodiment introduces the geometric concept of a skeleton curve, reducing the two-dimensional surface coverage problem to a one-dimensional line coverage problem. Specifically, the system calculates the centerline of the corridor, i.e., the skeleton curve, based on the left and right boundary points in the angular domain. Using the skeleton curve as a reference, a width function distributed along the curve is defined. The parameterized beam pattern model is no longer a static two-dimensional mask, but a set of geometric descriptions that vary along the skeleton curve parameter t, containing information in both the position (skeleton points) and extent (width) dimensions. This preserves the topological characteristics of the corridor, providing a natural mathematical framework for subsequently applying continuity constraints.

[0034] Based on the parametric beam pattern model, a beamforming weight optimization problem is constructed. The beamforming weight optimization problem includes coverage gain constraints distributed along the skeleton curve, as well as continuity constraints between adjacent sampling positions located on the skeleton curve.

[0035] Specifically, after establishing the parameterized model, the system needs to construct a mathematical optimization problem to solve for the antenna weights. The beamforming weight optimization problem is essentially the process of finding the optimal complex weight vector w, such that the generated antenna pattern meets specific performance indicators. The coverage gain constraint requires that within the target corridor area, i.e., the area defined by the parameterized model, the beam gain must be higher than a preset minimum service threshold G. min For example, 25 dBi is used to ensure the closure of the communication link. In the continuity constraint between adjacent sampling locations, since the skeleton curve discretizes the region into a series of ordered sampling points, the continuity constraint explicitly requires that the beam response (such as gain or field strength) between two adjacent sampling points along the skeleton direction must maintain a smooth transition, without abrupt jumps or breaks. This effectively prevents coverage holes at corridor bends, ensuring continuous communication service for users moving along the corridor. Compared to traditional methods that implicitly guarantee coverage by relying solely on dense sampling, the explicit continuity constraint reduces the number of sampling points required while maintaining effectiveness.

[0036] Solve the beamforming weight optimization problem to obtain the beamforming complex weight vector.

[0037] In this embodiment, the constructed optimization problem is typically a complex mathematical problem containing non-convex quadratic constraints. Numerical optimization algorithms can be used to find feasible or optimal solutions that satisfy all the above constraints. Depending on the specific application scenario and computing resources, various solution strategies can be adopted. As a basic implementation method, the semidefinite relaxation (SDR) technique can be used to transform the non-convex quadratic constraints into convex linear matrix inequalities, which are then solved using the interior-point method followed by rank-one restoration. As a preferred implementation method, the piecewise structure of the problem can also be utilized to perform efficient iterative solutions using the alternating direction multiplier method (ADMM). Regardless of the algorithm used, the final output is a complex weight vector w = [w1, w2, ..., w] of length N. N ] T Each element w n =a n e jΦn The amplitude and phase excitation correspond to the nth antenna element, where N is the total number of physical radiating elements in the array antenna. T Indicates transpose, a n Φ is the excitation amplitude of the nth antenna element. n For the phase excitation of the nth antenna element, a n e jΦn The complex excitation weight of the nth antenna element is represented by the amplitude a. n With phase Φ n Together they constitute.

[0038] The beamforming complex weight vector is loaded into a pre-configured beamforming network to drive the phased array antenna to radiate a communication beam that matches the shape of the target geographic corridor.

[0039] Specifically, after obtaining the numerical solution w, the beamforming complex weight vector is transmitted to the digital beamforming network (DBFN). The DBFN decomposes the complex weights into amplitude and phase commands, which control the variable gain amplifier (or attenuator) and phase shifter in each radiation channel, respectively. When all N units emit electromagnetic waves according to the commands, these electromagnetic waves interfere and superimpose with each other in space. Due to the weight design, constructive interference (signal enhancement) mainly occurs in the direction of the target geographic corridor, while destructive interference (signal cancellation) occurs in other directions. The energy distribution of the electromagnetic waves radiated by the antenna will exhibit a horseshoe-shaped or strip-shaped feature that highly matches the shape of the target geographic corridor, achieving efficient and precise coverage of a specific area.

[0040] like Figure 2 As shown, in one possible implementation, a parameterized beam pattern model distributed along the skeleton curve is established, including constructing a corridor width function, specifically:

[0041] Based on the corner domain representation, the left and right boundary point sets of the target geographic corridor are extracted, and the midpoint sequence of the corresponding boundary points is calculated;

[0042] By fitting the midpoint sequence, a continuous parameterized representation c(t) of the skeleton curve is obtained, where t is the normalized arc length parameter along the corridor direction.

[0043] Based on continuous parameterization, the angular spacing between the left and right boundary point sets from the satellite perspective is calculated, and a corridor width function distributed along the skeleton curve is constructed.

[0044] In this embodiment, the input data is the set of left boundary points P in the satellite angular coordinate system. L i and the right boundary point set P R i Where i=1,...,M represent the indices of discrete boundary points. The system calculates the geometric midpoint P for each pair of corresponding boundary points. C i =(P L i +P R i ) / 2. These midpoints constitute the initial discrete skeleton sequence. To facilitate subsequent differential calculations and piecewise processing, the discrete sequence needs to be transformed into a continuous function. Specifically, cubic spline interpolation or B-spline fitting methods can be used to fit the midpoint sequence into a smooth curve:

[0045] c(t) = (θ(t), Φ(t));

[0046] Where c(t) is the vector function of the skeleton curve; t is the normalized arc length parameter along the corridor direction, normalized to the interval [0, 1], i.e., t=0 corresponds to the starting point of the corridor, and t=1 corresponds to the ending point of the corridor; θ(t) is the azimuth component corresponding to parameter t; Φ(t) is the pitch component corresponding to parameter t. This not only ensures the continuity of the position but also the continuity of the tangent vector, i.e., the first derivative. Simultaneously, to describe the changes in the width of the corridor, the system calculates half the Euclidean distance w between each pair of boundary points (approximately the angular distance in the angular domain). i , i.e. w i =|P L i -P R iSimilarly, interpolation fitting is performed on the discrete width values ​​to obtain the continuous corridor width function w(t). For example, w(t) = d(t) / 2; where w(t) is the half-width of the corridor at parameter t, and d(t) is the angular distance between the left and right boundary points at parameter t from the satellite viewpoint. The entire horseshoe-shaped region is abstracted into a strip model centered at c(t) with a width of 2w(t), compressing the data volume while retaining the most crucial topological information.

[0047] like Figure 3 As shown, in a further embodiment, establishing a parameterized beam pattern model distributed along the skeleton curve also includes generating a set of coverage slice points, specifically:

[0048] At the sampling position of the skeleton curve, calculate the tangential vector and the normal vector, and establish a local coordinate system, where the tangential vector is along the extension direction of the skeleton curve, and the normal vector is perpendicular to the tangential vector.

[0049] Alternatively, multiple sampling positions t are determined on the skeleton curve based on a preset sampling interval or number of sampling points. k , where k is the index of the sampling position; at the sampling position t of the skeleton curve k Calculate the tangential and normal vectors and establish a local coordinate system.

[0050] Specifically, after obtaining the continuous curve c(t), the system selects K sampling positions t in the parameter space t∈[0,1] according to the accuracy requirements. k At each sampling position t k Calculate the tangential vector of the curve. According to the principles of differential geometry, the tangential vector e τ (t k The derivative vector of the curve at that point is defined as the normalized result.

[0051] e τ (t k )=c'(t k ) / ||c'(t k )||;

[0052] Where e τ (t k ) is the parameter t k The unit tangential vector at point c'(t) k ) represents the skeleton curve c(t) with respect to parameter t k The first derivative, i.e., c'(t) k )=(dθ / dt,dΦ / dt) | t=tk ;||c'(t k || is the magnitude of the first derivative vector. Define the normal vector e perpendicular to the tangent vector. n (tk In a two-dimensional angular domain plane, the following can be obtained by rotating the tangent vector counterclockwise by 90 degrees within the two-dimensional plane: if e τ =(u, v), then e n =(-v, u); where e n Let be the normal vector; v be the component of the tangential vector in the pitch direction; and u be the component of the tangential vector in the azimuth direction. The rotation operation ensures that the normal vector always points in the lateral direction of the corridor, forming a right-handed orthogonal coordinate system with the tangential vector. At each point on the skeleton curve, a locally orthogonal coordinate system {e} is established, rotating with the curve's curvature. τ e n The ability to adaptively follow the direction of the corridor is the geometric foundation for achieving structured sampling.

[0053] Based on the local coordinate system, sampling is performed in the direction of the normal vector and within the range determined by the corridor width function to generate a set of overlay slice points corresponding to the sampling positions.

[0054] In this embodiment, utilizing the established local coordinate system, the system no longer employs traditional global grid sampling, but instead uses slice-based sampling. Specifically, at the k-th skeleton point c(t) k At point ), the system only operates in the normal direction e. n (t k Linear sampling is performed on [-w(t). The sampling range is limited to [-w(t)]. k ), +w(t k The generated set of sampling points Ω is between )] in k It can be represented as:

[0055] Ω in k ={c(t k )+η·e n (t k )∣η∈[-w(t k ), w(t) k )]};

[0056] Where η is the offset along the normal direction.

[0057] In actual calculations, a finite number of uniformly distributed sampling points are taken as constraint points within a continuous interval. The number of sampling points can be determined based on the local width of the corridor and the required coverage accuracy. An appropriate sampling density can be selected according to actual needs.

[0058] The sampling points in the set constitute a cross-section or slice perpendicular to the corridor's orientation. Collecting all K slices covers the entire horseshoe-shaped region. This naturally adapts to the region's curved shape, and the ordered arrangement of the sampling points—that is, according to the skeleton sequence—makes defining adjacent slices and applying continuity constraints very intuitive and easy. In contrast, points generated by traditional global mesh sampling are disordered, making it difficult to identify which points are geometrically adjacent and to effectively apply continuity constraints.

[0059] According to one aspect of this application, it further includes: in the subsequent optimization process, it is necessary to calculate the response of the antenna array in any direction, which is determined by the array steering vector. The array steering vector a(θ, Φ) is a complex vector of length N, and its nth element is calculated using the following formula:

[0060] a n (θ,Φ)=g n (θ, Φ)·exp(j·2·π / λ*·(x n ·sin(θ)·cos(Φ)+y n ·sin(Φ)));

[0061] Where a n (θ, Φ) is the nth element of the steering vector, representing the complex response of the nth antenna element to electromagnetic waves in a specified direction; g n (θ, Φ) is the element pattern function of the nth antenna element, reflecting the radiation characteristics of that element; j is the imaginary unit; π is pi; λ* is the operating wavelength; x n Let y be the lateral position coordinate of the nth antenna element in the array coordinate system; n Let be the longitudinal position coordinates of the nth antenna element in the array coordinate system; θ be the azimuth angle in the target direction; Φ be the elevation angle in the target direction; exp be the natural exponential function; sin be the sine function; and cos be the cosine function. The steering vector contains all the information about the array geometry and the observation direction, and is the basis for calculating the antenna gain response.

[0062] Alternatively, the simplified formula for calculating the guiding vector a(θ, Φ) can also be:

[0063] a(θ, Φ) = [exp j·2π·d1·sin(θ)·cos(Φ) / λ* ,exp j·2π·d2·sin(θ)·cos(Φ) / λ* , ...,exp j ·2π·dN·sin(θ)·cos(Φ) / λ* ] T ;

[0064] Where a(θ, Φ) is the guiding vector corresponding to the direction (θ, Φ), j is the imaginary unit, π is pi, and d nθ represents the position offset of the nth array element relative to the reference array element, Φ represents the azimuth angle, Φ represents the elevation angle, λ* represents the operating wavelength, and N represents the total number of array elements.

[0065] It should be noted that the simplified representation of the steering vector applies to uniform linear arrays or when beam scanning is considered only in a single plane. For the general case of two-dimensional planar arrays, the steering vector is calculated using the aforementioned method involving x. n and y n The complete form of the coordinates of the two array faces.

[0066] like Figure 4 As shown, in a preferred implementation, constructing a beamforming weight optimization problem includes dividing the target geographic corridor into a predetermined number of sub-segments, specifically:

[0067] Calculate the curvature distribution of the skeleton curve along parameter t;

[0068] Based on curvature distribution, the target geographic corridor is divided into a predetermined number of sub-segments along the skeleton curve. The curvature integral value in each sub-segment is less than a preset curvature threshold, so that the local area corresponding to the sub-segment exhibits convex geometric characteristics.

[0069] In this embodiment, directly optimizing the entire curve when dealing with a highly curved horseshoe-shaped corridor is often very difficult because the entire region is non-convex. The principle of local convexity in differential geometry can be utilized, that is, any smooth curve can be approximated as a straight line at a sufficiently small scale. Specifically, the system calculates the curvature κ(t) of the skeleton curve c(t) at various points. The curvature reflects the degree of curvature of the curve. For example, the curvature κ(t) of the skeleton curve c(t) = (θ(t), Φ(t)) at parameter t is defined as:

[0070] κ(t) = |θ'(t)·Φ''(t) - Φ'(t)·θ''(t)| / [(θ'(t)) 2 + (Φ'(t)) 2 ] 3 / 2 ;

[0071] Where θ'(t) and Φ'(t) are the first derivatives of the curve coordinates with respect to the parameter t, and θ''(t) and Φ''(t) are the second derivatives. A larger curvature value indicates a more pronounced curvature at that point. An adaptive piecewise segmentation strategy is used to divide the parameter interval [0, 1] into S sub-segments [τ]. s-1 , τ s The division is based on controlling the cumulative bending amount within each sub-segment, which requires satisfying the integral inequality:

[0072] ∫ τs-1 τs |κ(t)|dt≤κ th;

[0073] Among them κ th This is a preset curvature integral threshold. It ensures that segments are denser in areas of sharp curvature and sparser in areas of flatness, and that each segment is geometrically approximately a straight line, thus being considered a convex region. This provides the geometric basis for subsequently converting non-convex gain constraints into convex constraints.

[0074] It should be noted that the curvature integral threshold κ th It is an adjustable engineering parameter used to control the cumulative degree of curvature within each segment. Its specific value can be determined based on the geometric characteristics of the corridor and the accuracy requirements of the local convex approximation.

[0075] According to one aspect of this application, after the corridor is segmented, it is necessary to calculate the geometric characteristic parameters of each sub-segment in order to subsequently establish a local coordinate system and constraints. Specifically, for the s-th sub-segment, its arc length L s It is obtained by integrating the magnitude of the derivative of the skeleton curve over this parameter interval. The specific calculation formula is as follows:

[0076] L s =∫ τs-1 τs |c'(t)|dt;

[0077] Where L s Let τ be the arc length of the s-th segment, with units consistent with angular coordinates; s-1 τ is the initial parameter value for the s-th segment; s Let be the termination parameter value of the s-th segment; c'(t) be the derivative vector of the skeleton curve at parameter t; || denotes the Euclidean norm of the vector; dt is the differential of the integration variable; ∫ is the sign of the definite integral. This integral can be calculated using numerical integration methods. Simultaneously, the average tangential direction of this sub-segment is calculated using the following formula:

[0078] e τ_bar_s =1 / (τ s -τ s-1 )·∫ τs-1 τs e τ (t)dt;

[0079] Where e τ_bar_s τ is the average tangential direction vector of the s-th segment; s τ is the termination parameter value for the s-th segment; s-1 e is the initial parameter value for segment s; τ(t) is the unit tangential vector at parameter t; dt is the differential of the integration variable; ∫ is the definite integral sign. The average tangential direction represents the overall orientation of the sub-segment and is used to establish the local reference coordinate system for that segment. In the local coordinate system, the coverage area of ​​the sub-segment can be approximately represented as a rectangular strip region. Let ξ be the coordinate along the local tangential direction and η be the coordinate along the local normal direction, then the coverage area of ​​the s-th segment is approximately represented as:

[0080] 0≤ξ≤L s , -w s ≤η≤w s ;

[0081] Where ξ is the local coordinate along the tangential direction; L s Let be the arc length of the s-th segment; η be the local coordinates along the normal direction; w s Let be the average half-width of the corridor in the s-th segment. The rectangular approximation makes the originally curved, non-convex region exhibit convex set characteristics locally, providing a geometric basis for subsequent convex optimization processing.

[0082] As an optional implementation, in each sub-segment, the coverage gain constraint is specifically configured as a field strength domain convexity constraint. The field strength domain convexity constraint requires that, after introducing the reference phase corresponding to the sub-segment, the real part of the product of the beamforming complex weight vector and the guide vector of the sampling point in the sub-segment after phase rotation is not lower than the field strength threshold corresponding to the preset gain threshold.

[0083] In other words, in each of the predetermined sub-segments, the coverage gain constraint is specifically configured as a field strength domain convexity constraint. The reference phase is determined based on the phase of the steering vector at the center of that sub-segment.

[0084] In this embodiment, since the original coverage constraint is based on power gain, i.e., |w H a| 2 ≥G min , is a non-convex quadratic inequality, where a is the guiding vector, and G is... min This represents the minimum power gain value corresponding to the gain threshold. The constraint can be transformed into a linear constraint in the field strength domain by utilizing the local convexity of the sub-segment. For the s-th sub-segment, the system introduces a reference phase ψ. s Preferably, the initial value of the reference phase corresponding to the sub-segment is determined by: calculating the parametric coordinates of the center position of the sub-segment on the skeleton curve; calculating the corresponding array steering vector based on the parametric coordinates; and extracting the phase of the array steering vector as the initial reference phase of the sub-segment. Specifically, the initial value of the reference phase can be calculated by calculating the center point c((τ) of the sub-segment. s-1 +τ s The guiding vector a at () / 2) center The initial value ψ of the reference phase is obtained from the phase. s 0= arg(w init H a center ), where arg() represents the argument of the complex number; w init The initial beamforming complex weight vector can be set as a uniform weight vector, meaning all elements have the same value, or as matched filter weights pointing towards the corridor center. A suitable initialization method can be selected based on the specific application scenario. Using the reference phase, the power constraint is transformed into a field strength real part constraint.

[0085] Re{e -jψs ·w H a(θ, Φ)}≥sqrt(G min );

[0086] Where a(θ, Φ) is the guiding vector of any sampling point within the sub-segment, sqrt(G min ) represents the corresponding field strength threshold, and w is the beamforming complex weight vector. H Let represent the conjugate transpose, and Re{} denote taking the real part of the complex number. The transformation process involves linear inequalities with respect to the weight vector w, defining a convex half-space. Through this transformation, the originally complex non-convex constraints are relaxed into a series of linear convex constraints, reducing the difficulty of the solution.

[0087] In a further embodiment, the continuity constraint between adjacent sampling positions is specifically configured as an inter-segment field strength continuity coupling constraint; the inter-segment field strength continuity coupling constraint requires that at the junction of two adjacent sub-segments, the difference between the complex field strength corresponding to the end of the previous sub-segment and the complex field strength corresponding to the beginning of the next sub-segment, or the difference between the complex field strength and the preset boundary reference field strength, does not exceed the preset coupling tolerance.

[0088] In this embodiment, since the entire corridor is divided into multiple independent segments, inter-segment coupling constraints need to be introduced to ensure that the final generated beam does not break or undergo phase abrupt changes at the junctions of these segments. Specifically, at the junction P of segment s and segment s+1... s =c(τ s At point ), the system introduces an auxiliary variable—the boundary reference field strength E. ref s The constraint condition requires that the actual complex field strength E(P) at this location be... s )=w H a(P s It must approximate the reference value: |w H a(P s )-E ref s | 2 , where a(P s Point P is the boundary point. sThe guiding vector at the point; or, under the augmented Lagrange framework, directly treated as a penalty term. The inter-segment field strength continuity coupling constraint acts as a glue, forcing adjacent sub-segments to maintain consistent amplitude and phase at the boundary, stitching all independent local convex optimization problems into a globally continuous solution.

[0089] Optionally, the inter-segment field strength continuity coupling constraint can be expressed as:

[0090] |E s，end -z s |≤ε and|E s+1，start -z s |≤ε;

[0091] Where E s，end Let E be the complex field strength at the end of the s-th sub-segment. s+1，start Let z be the complex field strength at the beginning of the (s+1)th subsegment. s Let be the boundary reference field strength at the s-th intersection, || denotes the complex modulus, and ε is the preset coupling tolerance.

[0092] In a further embodiment, solving the beamforming weight optimization problem includes:

[0093] An augmented Lagrangian function is constructed, which includes a pre-configured original optimization objective and an inter-segment field strength continuity coupling constraint as a penalty term. The augmented Lagrangian function is iteratively solved using the alternating direction multiplier method. Each round of the iteration includes alternating weight update steps, reference field strength update steps, and multiplier update steps. The weight update step, with fixed boundary reference field strength and Lagrangian multipliers, solves for the beamforming complex weight vector that minimizes the augmented Lagrangian function within the feasible region that satisfies the field strength domain convexity constraint.

[0094] In this embodiment, the system constructs an augmented Lagrangian function L of the following form. ρ :

[0095] L ρ (w,{E ref},{λ})=f(w)+∑(coupling constraint terms)+∑(penalty terms);

[0096] Where f(w) is the original optimization objective function, such as minimizing sidelobes; λ is the Lagrange multiplier; and ρ is the penalty parameter. The ADMM algorithm finds the optimal solution through the following three iterative steps:

[0097] Weight update steps: Fix boundary reference field strength E refUsing the Lagrange multiplier λ, solve for the beamforming complex weight vector w. Since the field strength constraint is convex and the coupling constraint is a quadratic term, the weight update is equivalent to solving a standard convex quadratic programming (QP) or second-order cone programming (SOCP) problem, which can be solved efficiently using optimizers (such as CVX, Gurobi).

[0098] Reference field strength update steps: Fix the beamforming complex weight vector w and the Lagrange multiplier λ, and update the boundary reference field strength E. ref It is an unconstrained quadratic minimization problem with an explicit closed-form solution:

[0099] E ref (s，k+1) =(1 / 2)·(w H a(P s )+E ref (s，k) -λ s k / ρ);

[0100] Where E ref (s，k+1) Let λ be the updated boundary reference field strength value corresponding to the s-th sub-segment in the (k+1)-th ADMM iteration; s k Let be the Lagrange multiplier corresponding to the s-th sub-segment in the k-th ADMM iteration. The closed-form solution is computed extremely quickly, which is key to the algorithm's efficiency. It should be noted that the above update uses an under-relaxed ADMM variant with a relaxation parameter of 0.5; the standard ADMM update formula can also be used.

[0101] Multiplier update step: Update the Lagrange multiplier λ to accumulate constraint violation errors and drive the algorithm to convergence, i.e.:

[0102] λ s k+1 =λ s k +ρ(w H a(P s )-E ref (s，k+1) ).

[0103] According to another aspect of this application, the specific mathematical form of the augmented Lagrange function is as follows:

[0104] L ρ (w,{E ref},{λ})=f(w)+Σ s=1 S-1 [λ s ·(E(p s ;w)-E ref_s )+ρ / 2·|E(ps ;w)-E ref_s | 2 ];

[0105] Where L ρ For the augmented Lagrangian function; ρ is the penalty parameter; w is the beamforming complex weight vector to be optimized; {E ref {} represents the set of all boundary reference field strengths; {λ} represents the set of all Lagrange multipliers; f(w) is the original optimization objective function, usually taken as the out-of-band maximum sidelobe level; Σ is the summation sign; s is the sub-segment index; S is the total number of sub-segments; λ s E(p) is the complex Lagrange multiplier corresponding to the s-th boundary coupling constraint; s w) represents the intersection point p between segment s and segment s+1. s The complex field strength E generated by the weight vector w at that point; ref_s Let be the reference field strength at the s-th boundary; || represents the modulus of the complex number.

[0106] It should be noted that the choice of the penalty parameter ρ directly affects the convergence speed and solution quality of the algorithm. A value that is too small will lead to slow convergence, while a value that is too large may cause numerical instability. In engineering practice, the typical range of ρ is 0.1 to 10, which can be adjusted according to the problem size and the number of constraints. It is usually advisable to start with ρ=1.

[0107] In a further embodiment, the iterative solution process also includes a phase adaptive update step, specifically:

[0108] Using the beamforming complex weight vector updated in the current round, the array response phase corresponding to the center position of the sub-segment is calculated; the array response phase is updated to the reference phase used by the sub-segment in the next round, and the definition of the field strength domain convexity constraint is corrected based on the reference phase to reduce the approximate error between the field strength domain convexity constraint and the preset original power gain constraint.

[0109] Specifically, the real part constraint of the electric field strength is an approximation of the original power constraint, i.e., it assumes that the phase is known. To eliminate the error introduced by the approximation, a closed-loop feedback mechanism is introduced. After each round or every few rounds of the ADMM iteration, the latest weight w calculated at the current time is used. k+1 The actual beam phase at the center point of each sub-segment is recalculated, and the new phase is used as the reference phase ψ for the next round. s That is, the new reference phase ψ s new for:

[0110] ψ s new =angle(w (k+1)H a center ).

[0111] As iterations proceed, the reference phase gradually approaches the true phase of the optimal solution, making the convex constraints increasingly accurate in describing the original non-convex power constraints. An adaptive update mechanism enables the algorithm to converge to a high-quality local optimum, overcoming the relaxation gap problem inherent in traditional convex relaxation methods.

[0112] According to another aspect of this application, the specific calculation formula for the multiplier update step is as follows:

[0113] λ s k+1 =λ s k +ρ·(E(p s ;w k+1 )-E ref_s k+1 );

[0114] Where λ s k+1 λ is the updated value of the s-th multiplier after the (k+1)-th iteration; s k The current value of the s-th multiplier in the k-th iteration is ρ; the penalty parameter is E(p). s ;w k+1 ) represents the weight vector w obtained using the (k+1)th iteration. k+1 The calculated boundary point p s Field strength at E ref_s k+1 This represents the reference field strength after the (k+1)th iteration update. It indicates that the multiplier is updated cumulatively according to the constraint violation amount, and the multiplier no longer changes when the constraint is precisely satisfied.

[0115] According to another aspect of this application, it also includes: a phase update period and a convergence determination. Specifically, in the actual iteration process, the update of the reference phase does not need to be performed in every iteration round, but can be performed once every few rounds to reduce computational overhead. Let the phase update period be N. phase That is, N per iteration phase A reference phase update is performed after each round. N phase Typical values ​​for are 5 to 10. Smaller values ​​result in more accurate phase tracking but increased computation, while larger values ​​improve computational efficiency but may slow convergence. The iteration process terminates using a dual convergence condition. The algorithm terminates and outputs the current weight vector as the final solution when both of the following conditions are met simultaneously. The first condition is the boundary coupling error convergence condition:

[0116] max s |E(p s ;w k+1 )-E ref_s k+1 |<εconv ;

[0117] Where max s This indicates the operation of taking the maximum value over all boundary points s; E(p) s ;w k+1 ) represents the boundary point p after the (k+1)th iteration. s Field strength at E ref_s k+1 The reference field strength after the (k+1)th iteration; || represents the modulus of the complex number; ε conv This is the field strength coupling convergence threshold, typically ranging from 1% to 5% of the target field strength.

[0118] The second condition is the convergence condition of the weight vector:

[0119] ||w k+1 -w k ||2<δ conv ;

[0120] Where w k+1 w is the weight vector obtained in the (k+1)th iteration. k Let be the weight vector for the k-th iteration; || ||2 represents the 2-norm of the vector; δ conv This is the convergence threshold for weight changes, typically set to 1×10. -4 Up to 1×10 -3 .

[0121] When both of the above conditions are met, it indicates that the algorithm has converged to a stable solution. At this point, the boundary coupling constraints have been effectively satisfied, and the continuity of the beam along the corridor direction is guaranteed.

[0122] The continuity constraint between adjacent sampling locations can also be configured as follows:

[0123] At any two adjacent sampling positions on the skeleton curve, calculate the absolute value of the difference between the beam gain corresponding to the previous sampling position and the beam gain corresponding to the next sampling position; the absolute value of the difference is required not to exceed the preset gain fluctuation tolerance along the path, so as to limit the drastic jump in signal strength along the skeleton curve direction.

[0124] In this embodiment, since a segmentation strategy is not employed, the system directly applies a global continuity constraint to all sampling points distributed along the skeleton curve. Specifically, for the k-th and (k+1)-th sampling positions on the skeleton curve, their corresponding beam gains are G and G, respectively. k and G k+1 The gain G is a quadratic function of the weight vector w, i.e., G = |w|. H a| 2 The system calculates the absolute value of the gain difference between these two points, |G|. k -G k+1| and requires that this value be less than or equal to the preset tolerance ΔG cont Tolerance ΔG cont The specific value is determined based on the communication system's requirements for gain smoothness along the path. Mathematically, this is represented by a series of quadratic inequality constraints: -ΔG cont ≤w H (A k -A k+1 w≤ΔG cont A k =a k a k H It is the guiding matrix at the k-th position, a k This is the steering vector at the k-th position. It forces the main lobe response of the beam to change smoothly along the corridor direction, avoiding sudden fluctuations and ensuring the stability of the received signal for the user terminal during movement.

[0125] Furthermore, solving the beamforming weight optimization problem includes:

[0126] By introducing positive semidefinite matrix variables to replace the beamforming complex weight vector, the non-convex quadratic constraint on the beamforming complex weight vector is transformed into a linear constraint on the positive semidefinite matrix variables.

[0127] Ignoring the rank-one constraint of the positive semidefinite matrix variables, the beamforming weight optimization problem is relaxed into a semidefinite programming problem and solved to obtain the optimal positive semidefinite matrix.

[0128] If the rank of the optimal positive semi-definite matrix is ​​greater than one, then the beamforming complex weight vector that satisfies the coverage gain constraint and the continuity constraint can be recovered from the optimal positive semi-definite matrix using the Gaussian randomization method.

[0129] This embodiment describes the core process of the semi-definite relaxation (SDR) algorithm. The original beamforming problem includes a form w H Aw≥γ is a non-convex quadratic constraint, where A is the steering matrix and γ is the threshold of the coverage gain constraint. To solve this problem, a new matrix variable W=ww is introduced. H By definition, W is a positive semi-definite (PSD) matrix of beamforming weights, and its rank is 1. Using the property of the matrix trace, i.e., w... H Aw=Tr(Aww H =Tr(AW), where Tr() represents the trace of the matrix; all quadratic constraints can be transformed into linear constraints on W. For example, the cover gain constraint becomes Tr(AW)≥γ, and the continuity constraint becomes |Tr((AW)|A|W ... k -A k+1 W)|≤ΔG contAt this point, the optimization problem is transformed into finding a positive semi-definite matrix W that satisfies the linear constraint and has a rank of 1. The rank-1 constraint rank(W) = 1 is still non-convex. The SDR method directly ignores (relaxes) the rank constraint, retaining only the convex constraint that W is a positive semi-definite matrix. After relaxation, the original problem becomes a standard semi-definite programming (SDP) problem, which can be solved in polynomial time using convex optimization toolboxes (such as CVX, SeDuMi, or MOSEK). If the rank of the obtained W* is exactly 1, the optimal weight vector w* can be directly obtained through eigenvalue decomposition. However, in practical applications, especially when there are many constraints, the rank of W* is often greater than 1. In this case, Gaussian randomization is needed to recover the solution. The specific process is as follows: generate a large number of random vectors z based on the distribution z ~ CN(0, W*). m For example, 1000, each z m Mapping to the feasible region, for example by scaling it to satisfy power constraints, the vector that minimizes the objective function such as the sidelobe level is selected as the final approximate solution w.

[0130] According to one aspect of this application, the concept of a direction matrix needs to be introduced during the positive semidefinite relaxation process. The direction matrix A(θ, Φ) is defined as the outer product of the guiding vector and its conjugate transpose, and the specific calculation formula is as follows:

[0131] A(θ,Φ)=a(θ,Φ)·a(θ,Φ) H ;

[0132] Where A(θ, Φ) is an N×N Hermitian matrix; N is the number of antenna elements; and a(θ, Φ) is the steering vector in the direction (θ, Φ). H This represents the conjugate transpose operation, i.e., taking the conjugate first and then transposing. Using the direction matrix, the original power gain expression can be rewritten in the form of matrix trace. Let the weight matrix W = w·w H Then the formula for calculating the power gain in the direction (θ, Φ) is:

[0133] G(θ, Φ)=Tr(A(θ, Φ)·W);

[0134] Where G(θ, Φ) is the power gain in the direction (θ, Φ); Tr() represents the trace operation of the matrix, i.e., the sum of the diagonal elements; A(θ, Φ) is the direction matrix; and W is the weight matrix. Transforming the quadratic function in w into a linear function in W is the key step in achieving convexity transformation.

[0135] In a further embodiment, after recovering the beamforming complex weight vector using the Gaussian randomization method, the method further includes:

[0136] Using the recovered beamforming complex weight vector as the initial point, a local search is performed on the beamforming weight optimization problem using sequential convex programming or gradient projection methods to further reduce the maximum beam gain of the out-of-band sampling point set.

[0137] In this embodiment, the solution w obtained by Gaussian randomization is only an approximate solution to the original non-convex problem and may not be a local optimum. To further improve performance, a local refinement process can be introduced. The weight vector recovered by randomization is used as the initial point w0, and iterative optimization is performed using the Sequential Convex Programming (SCP) method. Specifically, in each iteration, the non-convex constraints are linearized and approximated near the current point, for example, by using a first-order Taylor expansion, to construct a local convex optimization subproblem. Solving this subproblem yields the descent direction, and the weight vector is updated. This process is repeated until convergence. This effectively corrects the error caused by SDR relaxation, further reduces the sidelobe level, and can bring a performance gain of 1dB to 3dB, ensuring the rigor and high performance of the final solution in engineering applications.

[0138] This embodiment can effectively handle non-convex quadratic constraint problems and is suitable for scenarios where the corridor curvature is relatively small or where computing resources allow for high-dimensional matrix operations.

[0139] According to one aspect of this application, constructing the beamforming weight optimization problem further includes setting out-of-band sidelobe suppression constraints; setting out-of-band sidelobe suppression constraints includes:

[0140] Based on the local coordinate system, sampling is performed in the direction of the normal vector, and in the area beyond the coverage area corresponding to the coverage slice and separated by a preset protection distance, to generate an out-of-band sampling point set; it is required that the beam gain of each sampling point in the out-of-band sampling point set does not exceed the preset sidelobe level upper limit, or that the maximum beam gain in the out-of-band sampling point set is minimized.

[0141] In other words, the beam gain of each sampling point in the out-of-band sampling point set is required not to exceed the preset upper limit of the sidelobe level, or the maximum beam gain in the out-of-band sampling point set is minimized as the optimization objective.

[0142] In this embodiment, to prevent satellite signals from interfering with non-target areas, power limiting is required for the region outside the main lobe. Preferably, a local coordinate system is used to accurately define the out-of-band region. Specifically, for each position t on the skeleton curve... k The out-of-band region is defined as the region whose absolute value of the normal distance η is greater than the corridor half-width w(t). k ) plus protection interval Δ guard The area; where the preset protection distance Δ guard This is used to provide physical roll-off space between the main lobe and side lobes, typically ranging from 0.1 to 0.5 degrees. Sample points collected in these regions constitute the out-of-band sampling set Ω.out For these sampling points, two forms of constraints can be imposed. The first is a hard constraint, which requires that the gain G(θ, Φ) ≤ G at all out-of-band points. SL G SL The first is a preset upper limit for the sidelobe, such as -20dB. The second is the Min-Max criterion, which uses minimizing the maximum out-of-band gain as the objective function of the optimization problem: minmaxG(θ, Φ), (θ, Φ)∈Ω out In practical engineering, they are usually used in combination, that is, minimizing the maximum sidelobe in the objective function while setting a bottom-line threshold that must be met to ensure compliance with the radio regulations of the International Telecommunication Union (ITU).

[0143] In one possible embodiment, the beamforming weight optimization problem further includes a power constraint on the beamforming complex weight vector, configured as follows:

[0144] Total power constraint: The square of the L2 norm of the beamforming complex weight vector must not exceed the preset maximum transmit power;

[0145] Constant mode constraint: requires that the modulus of each element in the beamforming complex weight vector be constant, and only the phase is adjusted.

[0146] In this embodiment, two power constraint configurations are provided for different satellite payload hardware characteristics. The first is a total power constraint, suitable for systems with limited total power but dynamically allocable unit power. Its mathematical expression is ||w|| 2 ≤P max , where ||w|| 2 Let P be the 2-norm of the beamforming complex weight vector w. max This is the preset maximum transmit power. This is a convex constraint (inside the sphere) that can be directly incorporated into the SDR or ADMM solution process.

[0147] The second type is constant mode constraint, which is particularly critical for phased array antennas employing solid-state power amplifiers (SSPA) or traveling wave tube amplifiers (TWTA) operating in saturation mode. Constant mode constraint requires |w n |=1 (assuming normalized magnitude), meaning the weight vector consists only of the phase term e jΦn Composition, where w n Let w be the nth element of the beamforming complex weight vector w. This is a non-convex equality constraint. Handling constant mode constraints is very intuitive within the ADMM framework. Specifically, the intermediate solution w is obtained during the weight update step. opt Then, it is mapped onto the unit circle through a projection operation, that is, the updated beamforming complex weight vector w is obtained. new For: w new =exp(j*angle(w optThis avoids complex non-convex optimization solutions, allowing the algorithm to maintain efficient convergence while meeting hardware requirements. By supporting constant modulus constraints, it can be directly applied to low-cost, high-energy-efficiency spaceborne phased array systems.

[0148] According to one aspect of this application, the fitting method for the skeleton curve is preferably implemented using cubic spline interpolation. Cubic spline interpolation can ensure the second-order continuous differentiability of the curve, so that a smooth result can be obtained when calculating the curvature later. For cases with a large number of control points or requiring local adjustments, B-spline fitting can also be used. The above methods can be implemented using existing technologies.

[0149] In one embodiment of this application, a method for generating satellite communication beams may further include: obtaining the boundary coordinates of a target geographic corridor and converting them into azimuth and elevation angle ranges in a satellite coordinate system; defining a desired horseshoe-shaped beam pattern based on the converted angular domain range, wherein the gain requirement is high and uniform within the target angular domain, and the gain requirement is significantly suppressed outside the target angular domain; constructing an optimization problem with a beamforming weight vector as the variable, the objective function of which is to maximize the average gain within the target coverage area while minimizing the maximum sidelobe level in non-target areas (especially sensitive interference areas); solving the optimization problem to obtain the optimal complex weight coefficients; and loading the complex weight coefficients onto a digital beamforming network to drive a phased array antenna to radiate a horseshoe-shaped beam that highly matches the shape of the geographic corridor. Preferably, the horseshoe-shaped beam can be a closed approximate U-shape or C-shape, or an open, curved strip, the shape of which is entirely determined by the input geographic boundary information. The system supports the simultaneous generation of multiple horseshoe-shaped beams to cover multiple independent geographic corridors, and the isolation between beams is ensured through optimization algorithms. The beam control and management system can dynamically update target geographic information according to schedules or instructions, thereby enabling dynamic reconstruction of coverage areas and allowing satellite communication resources to be scheduled on demand.

[0150] In another embodiment of this application, a system for generating satellite communication beams is also provided, comprising:

[0151] Multi-beam phased array antenna: Deployed on communication satellites, it contains multiple radiating elements whose phase and amplitude can be independently controlled;

[0152] Digital beamforming network: electrically connected to a multi-beam phased array antenna, used to calculate and generate excitation signals for each radiating element based on beam weighting coefficients;

[0153] The beam control and management system includes a geographic information database, a beam shape calculation module, and a weight generation module. The geographic information database stores the boundary coordinates of the target coverage area. The beam shape calculation module maps the boundary coordinates to the angular domain coverage requirements from a satellite perspective and calculates the desired horseshoe-shaped beam spatial energy distribution map. Based on the spatial energy distribution map, the weight generation module uses a convex optimization algorithm or an iterative algorithm to calculate the complex weight coefficients required for the digital beamforming network. The complex weight coefficients make the main lobe shape of the synthesized beam match the target geographic area, while forming deep nulls or extremely low sidelobes in non-target areas.

[0154] In a detailed embodiment, a high-throughput satellite (HTS) equipped with an active phased array antenna (AESA) is used as a platform to provide communication services for a key waterway along the Maritime Silk Road. The ground operations center delineates the curved waterway areas requiring coverage based on the waterway map and generates boundary coordinate files. These boundary coordinate files are uploaded to the satellite's beam control and management system via a gateway station. The beamform calculation module in the system, combined with the satellite's real-time orbit and attitude data, converts the ground coordinates into a two-dimensional coverage matrix of azimuth and elevation angles from the perspective of the onboard antenna. The weight generation module initiates an optimization algorithm. The algorithm uses the excitation weights of each phased array element as variables, with the joint optimization objective of maximizing gain flatness within the corresponding angular domain of the waterway and minimizing peak sidelobe levels in the area outside the waterway, and performs iterative calculations. After convergence, a set of optimal weight coefficients is obtained and sent to the digital beamforming network (DBFN). The DBFN generates corresponding amplitude and phase control signals based on the weight coefficients, driving each radiating element of the multi-beam phased array antenna. The antenna ultimately radiates a horseshoe-shaped beam, whose main lobe perfectly envelops the target flight path, while signal energy in other directions is drastically suppressed. Results show that, with the same total radiated power, the horseshoe-shaped beam in this embodiment provides an average signal strength approximately 5 dB higher than a traditional circular beam on the target flight path. In neighboring countries (interference directions), the signal strength is reduced by more than 15 dB, demonstrating a significant advantage. This invention is applicable not only to fixed-orbit satellites but also to low Earth orbit (LEO) constellations. Within a constellation, each satellite can dynamically generate a horseshoe-shaped beam matching its instantaneous coverage area, thereby achieving continuous relay coverage of global strip-shaped regions.

[0155] In an exemplary embodiment, the satellite orbital altitude is set to 35,786 km (geosynchronous orbit) in the simulation scenario. The operating frequency is 20 GHz (Ka band), corresponding to a wavelength λ* of 15 mm. The onboard multi-beam phased array antenna uses a 32×32 uniform rectangular planar array with an element spacing of 0.5λ*. The target geographic corridor is set as an S-shaped curved area, with an azimuth span of approximately 10 degrees and an elevation span of approximately 5 degrees within the angular domain of the satellite's view. The average angular width of the corridor is set to 0.5 degrees. The optimization parameters set in the simulation are as follows: in-band minimum gain threshold G min Set to 25 dBi; set ΔG cont Corresponding to a gain ripple of 1.5 dB; sidelobe suppression target G SL Set to 0 dBi, which is approximately 25 dB lower than the main lobe peak value. Protection interval Δ guard The curvature was set to 0.2 degrees. A piecewise convexity transformation (PPT) + ADMM algorithm was used for solving the problem. The skeleton curve was extracted, and the corridor was divided into 5 sub-segments based on the curvature threshold. During the ADMM iteration process, the penalty parameter ρ was set to 1.0, and the number of iterations was set to 200.

[0156] Simulation results show that after the optimized complex weighted vector drives the antenna, the resulting main lobe perfectly matches the S-shaped corridor. In terms of coverage performance, the average gain within the target area reaches 26.5 dBi, and the minimum gain is 25.1 dBi, meeting the coverage requirements. Regarding continuity, the maximum gain jump between adjacent sampling points along the skeleton curve is only 0.8 dB, far below the 1.5 dB tolerance, proving the effectiveness of the continuity constraint. In contrast, if this continuity constraint is removed, i.e., only traditional point sampling optimization is used, a 3.5 dB gain pit appears at the point of most severe corridor curvature, forming a significant coverage discontinuity. Regarding sidelobe suppression, this embodiment reduces the maximum out-of-band sidelobe level to -2.3 dBi, achieving excellent isolation. Compared to the SDR method, the ADMM method reduces computation time while achieving the same sidelobe level because ADMM decomposes large-scale matrix operations into multiple small-scale sub-problems. Meanwhile, the ADMM method can directly handle constant modulus constraints, while the SDR method requires an extremely complex rank-one recovery process when handling constant modulus constraints, often resulting in significant performance loss. In summary, this embodiment balances computational efficiency and hardware adaptability while ensuring the complexity and continuity of the covered shape.

[0157] This invention employs skeleton curve parametric modeling technology. By extracting the corridor centerline, the originally complex two-dimensional irregular region coverage problem is reduced in dimensionality to a one-dimensional ordered parametric problem distributed along a curve. This reduces the number of optimization variables and constraints, lowering the computational load of the algorithm. An adaptive piecewise segmentation and field strength domain convexity transformation strategy based on curvature is introduced. Geometric segmentation decomposes the overall non-convex region into multiple locally approximately convex sub-segments, and the non-convex power gain constraint is transformed into a linear field strength real part constraint. This is combined with the ADMM algorithm for collaborative solving, eliminating the convergence barrier caused by non-convexity and enabling the algorithm to quickly find a high-quality global optimal solution. Explicit piecewise continuity constraints are applied, including gain difference limits between adjacent points or field strength coupling between segments, forcing a smooth transition of beam energy along the corridor direction. This effectively prevents coverage discontinuities and ensures the continuity of communication services.

[0158] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A method for generating satellite communication beams, characterized in that, include: Obtain the boundary coordinates of the target geographic corridor and convert them into an angular domain representation in a satellite coordinate system; Based on the angular domain representation, the skeleton curve of the target geographic corridor is extracted, and a parameterized beam pattern model distributed along the skeleton curve is established. Based on the parametric beam pattern model, a beamforming weight optimization problem is constructed, which includes coverage gain constraints distributed along the skeleton curve and continuity constraints between adjacent sampling positions located on the skeleton curve. Solve the beamforming weight optimization problem to obtain the beamforming complex weight vector; The beamforming complex weight vector is loaded into a pre-configured beamforming network to drive the phased array antenna to radiate a communication beam that matches the shape of the target geographic corridor.

2. The method according to claim 1, characterized in that, Establish a parameterized beam pattern model distributed along the skeleton curve, including constructing the corridor width function, specifically: Based on the corner domain representation, the left and right boundary point sets of the target geographic corridor are extracted, and the midpoint sequence of the corresponding boundary points is calculated; By fitting the midpoint sequence, a continuous parameterized representation c(t) of the skeleton curve is obtained, where t is the normalized arc length parameter along the corridor direction. Based on continuous parameterization, the angular spacing between the left and right boundary point sets from the satellite perspective is calculated, and a corridor width function distributed along the skeleton curve is constructed.

3. The method according to claim 2, characterized in that, Establishing a parameterized beam pattern model distributed along the skeleton curve also includes generating a set of coverage slice points, specifically: At the sampling position of the skeleton curve, calculate the tangential vector and the normal vector, and establish a local coordinate system, where the tangential vector is along the extension direction of the skeleton curve, and the normal vector is perpendicular to the tangential vector. Based on the local coordinate system, sampling is performed in the direction of the normal vector and within the range determined by the corridor width function to generate a set of overlay slice points corresponding to the sampling positions.

4. The method according to claim 2, characterized in that, The beamforming weight optimization problem is constructed, which includes dividing the target geographic corridor into a predetermined number of sub-segments, specifically: Calculate the curvature distribution of the skeleton curve along parameter t; Based on curvature distribution, the target geographic corridor is divided into a predetermined number of sub-segments along the skeleton curve. The curvature integral value in each sub-segment is less than a preset curvature threshold, so that the local area corresponding to the sub-segment exhibits convex geometric characteristics.

5. The method according to claim 4, characterized in that, In each sub-segment, the coverage gain constraint is specifically configured as a field strength domain convexity constraint; Field strength domain convexity constraint requirement: After introducing the reference phase corresponding to the sub-segment, the real part of the product of the beamforming complex weight vector and the guide vector of the sampling point in the sub-segment after phase rotation shall not be lower than the preset field strength threshold corresponding to the gain threshold.

6. The method according to claim 5, characterized in that, The continuity constraint is specifically configured as an inter-segment field strength continuity coupling constraint; Inter-segment field strength continuity coupling constraint requirement: at the boundary between two adjacent sub-segments, the difference between the complex field strength corresponding to the end of the previous sub-segment and the complex field strength corresponding to the beginning of the next sub-segment, or the difference between the complex field strength and the preset boundary reference field strength, shall not exceed the preset coupling tolerance.

7. The method according to claim 6, characterized in that, Solving the beamforming weight optimization problem includes: Construct an augmented Lagrangian function, which includes a pre-configured original optimization objective and an inter-segment field strength continuity coupling constraint as a penalty term; The augmented Lagrangian function is solved iteratively using the alternating direction multiplier method. Each round of the iterative solution includes alternating weight update steps, reference field strength update steps, and multiplier update steps. The weight update step involves solving for the beamforming complex weight vector that minimizes the augmented Lagrange function within the feasible region that satisfies the field strength domain convexity constraint, while keeping the boundary reference field strength and Lagrange multipliers fixed.

8. The method according to claim 7, characterized in that, The iterative solution process also includes a phase adaptive update step, specifically: Using the beamforming complex weight vector updated in the current round, calculate the array response phase corresponding to the center position of the sub-segment; The array response phase is updated to the reference phase used by the segment in the next round, and the definition of the field strength domain convexity constraint is corrected accordingly to reduce the approximate error between the field strength domain convexity constraint and the preset original power gain constraint.

9. The method according to claim 1, characterized in that, Continuity constraints can also be configured as follows: At any two adjacent sampling positions of the skeleton curve, calculate the absolute value of the difference between the beam gain corresponding to the previous sampling position and the beam gain corresponding to the next sampling position. The absolute value of the difference must not exceed the preset gain fluctuation tolerance along the path.

10. The method according to claim 9, characterized in that, Solving the beamforming weight optimization problem includes: By introducing positive semidefinite matrix variables to replace the beamforming complex weight vector, the non-convex quadratic constraint on the beamforming complex weight vector is transformed into a linear constraint on the positive semidefinite matrix variables. Ignoring the rank-one constraint of the positive semidefinite matrix variables, the beamforming weight optimization problem is relaxed into a semidefinite programming problem and solved to obtain the optimal positive semidefinite matrix. If the rank of the optimal positive semi-definite matrix is ​​greater than one, then the beamforming complex weight vector that satisfies the coverage gain constraint and the continuity constraint can be recovered from the optimal positive semi-definite matrix using the Gaussian randomization method.

Images