Adaptive point cloud reconstruction and ISAR simulation imaging method for space target

Abstract

The application relates to the technical field of radar signal processing, and provides a kind of adaptive point cloud reconstruction and ISAR simulation imaging method based on space target, comprising: obtaining three-dimensional model data of space target;Adopting adaptive point cloud rule sampling strategy based on geometric features, extracting key structural points with physical scattering significance from the surface of three-dimensional model;Performing structured cleaning and statistical filtering operation on the key structural points obtained by sampling;The point cloud after the structured cleaning and statistical filtering operation is used as an ideal scattering center set, and a point target model is used to simulate received echo signals;Distance alignment and phase compensation are carried out;Using an improved minimum variance phase correction method, motion compensation and phase correction are carried out;Two-dimensional Fourier transform is carried out on the data after motion compensation and phase correction, and the final two-dimensional ISAR image is output.The application can provide high-precision and high-efficiency simulation and analysis means for space target identification and tracking.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of radar signal processing, and in particular to a method for adaptive point cloud reconstruction and ISAR simulation imaging based on space targets. BACKGROUND

[0002] With the acceleration of low-orbit constellation deployment, the surge in the number of on-orbit spacecraft, and the increasingly serious problem of space debris, automatic identification, classification and state perception of space targets have become one of the core tasks in the field of current space situational awareness (SSA). Space-based ISAR (Inverse Synthetic Aperture Radar) imaging is widely considered as a key means to realize the structural perception and identity recognition of space targets due to its advantages of long distance, all-weather, high resolution, etc.

[0003] However, the development of automatic recognition algorithms for complex space targets (such as spacecraft bodies, component debris, combined structures, etc.) highly depends on large-scale, high-quality, and clearly-structured radar image data. These data are not only used for algorithm training, but also support performance verification and robustness evaluation. Currently, real ISAR images are extremely limited in public access due to orbital conditions, target controllability, payload constraints, and security factors, which makes it difficult to meet the needs of data-driven methods such as deep learning.

[0004] Therefore, it is urgent to build an ISAR simulation imaging method for real structure modeling, which can generate large-scale high-fidelity radar image data under controllable conditions, and provide data support for the development and verification of target recognition, classification, structure analysis, etc.

[0005] However, the current ISAR simulation imaging technology still faces key challenges: (1) Most methods only use regular geometric bodies or simplified satellite models for scattering modeling, which is difficult to accurately reflect the complex structure, fine components and multi-scale features of real space targets, resulting in distortion of the simulated image and inconsistency of the structure distribution with the actual one. (2) In the face of irregular target topography and complex scattering structure, traditional sidelobe suppression methods (such as frequency domain windowing, CLEAN algorithm, etc.) cannot effectively suppress artifacts and cross-coupling terms, and the image has serious speckle and structure aliasing, which limits the effectiveness of the image for algorithm training. (3) It relies on special modeling tools and specific simulation platforms, with single format support, lack of automatic modeling and batch simulation capabilities, and is difficult to extend to multi-target, multi-scenario, and multi-frequency tasks.

[0006] Therefore, it is necessary to provide an ISAR imaging method capable of realizing accurate expression of a target based on a real target structure and combining point cloud reconstruction technology, and on this basis, capable of effectively suppressing sidelobes, improving imaging quality, and solving problems such as image blurring and artifact enhancement, which has important strategic significance for enhancing space situation awareness capability.

[0007] A patent with the application publication number CN109188384A discloses an electromagnetic simulation method for dynamic observation of space targets, which relates to establishing multiple coordinate systems using real orbits and observation geometry, obtaining electromagnetic scattering of space targets through large facet physical optics combined with FEKO electromagnetic simulation at the best observation time and viewing angle, and realizing two-dimensional image reconstruction based on range-Doppler algorithm. The method provides a relatively complete implementation path in terms of orbit modeling, attitude control, and basic electromagnetic echo simulation process, and has high engineering feasibility and simulation consistency. However, when the method is applied to real target modeling and imaging simulation at high frequency bands (such as terahertz bands), the following significant deficiencies still exist: (1) FEKO only supports a small number of model formats (such as.sat and.stl), among which the.stl format can only be imported in grid form and does not support post-editing or repair; it is not compatible with mainstream three-dimensional model formats (such as.glb,.ply,.obj, etc.), resulting in a tedious actual engineering model import process; (2) In terahertz or high frequency band simulation, the model facet size is required to be much smaller than the wavelength (typically λ / 3-λ / 10), and tens of thousands or even millions of facets need to be divided for complex spacecraft models. Too many facets will result in high simulation memory occupation, long calculation time, convergence difficulty or failure, which is particularly obvious in high-precision imaging simulation.

[0008] Therefore, there is an urgent need for a flexible, efficient, and scalable space target ISAR imaging simulation method that supports multi-format three-dimensional model input, combines adaptive point cloud construction and simplified scattering modeling mechanism, and realizes high-precision motion compensation to adapt to the needs of complex space target identification and algorithm training under future space-based platforms, terahertz bands, and high dynamic environments, while maintaining physical interpretability. SUMMARY

[0009] The present application mainly solves the technical problems of the prior art, such as excessive calculation amount of high-frequency imaging, difficulty in constructing a scattering point model, insufficient precision, and significant difference between simulated imaging and real imaging, and proposes an adaptive point cloud reconstruction and ISAR simulation imaging method for space targets, extracts key structure point clouds by compatible multi-format models (such as.glb,.ply, etc.), and performs scattering modeling in combination with an ideal scattering point model, thereby significantly reducing modeling complexity and improving image realism; meanwhile, an improved motion compensation and phase correction mechanism is introduced, an sidelobe suppression strategy is optimized, and the structure fidelity and imaging consistency of the simulated image are improved. The method has high flexibility, high simulation efficiency, and good engineering adaptability, and is suitable for space target identification and algorithm training tasks under space-based platforms, terahertz bands, and high dynamic environments.

[0010] The present application provides an adaptive point cloud reconstruction and ISAR simulation imaging method for space targets, comprising:

[0011] Step S1, acquiring three-dimensional model data of a space target;

[0012] Step S2, adopting an adaptive point cloud rule sampling strategy based on geometric features to extract key structure points with physical scattering significance from a three-dimensional model surface;

[0013] Step S3, performing a structured cleaning and statistical filtering operation on the key structure points obtained in step S2;

[0014] Step S4, taking the point cloud after the structured cleaning and statistical filtering operation in step S3 as an ideal scattering center set, and adopting a point target model to simulate received echo signals;

[0015] Step S5, adopting an improved cross-correlation-based distance alignment and phase compensation method to perform distance alignment and phase compensation;

[0016] Step S6, considering signal-to-noise ratio, phase stability, and multi-point weighting mechanism, adopting an improved minimum variance phase correction method to perform motion compensation and phase correction;

[0017] Step S7, performing two-dimensional Fourier transform on the data after motion compensation and phase correction, and outputting a final two-dimensional ISAR image.

[0018] Further, step S1 comprises steps S11-S13:

[0019] Step S11, acquiring a three-dimensional model of a spacecraft through a public data source;

[0020] Step S12, acquiring a three-dimensional model of a space debris through a ground super-high-speed impact test;

[0021] Step S13, uniform scale control is performed on the space target of different scales.

[0022] Further, step S13 includes steps S131-S132:

[0023] Step S131, for the space target with size greater than the size threshold range, downsizing processing is performed under the premise of keeping the structural proportion consistent;

[0024] Step S132, for the space target with size less than the size threshold range, equal proportion enlargement is performed.

[0025] Further, step S2 includes steps S21-S23:

[0026] Step S21, edge feature extraction and edge point sampling: based on the normal gradient or curvature change of the surface of the three-dimensional model, edge detection is performed, and structural edge regions are identified, and high-density sampling is performed in the structural edge regions;

[0027] Step S22, geometric feature point extraction and high-curvature point sampling: geometric feature points in the target three-dimensional model are detected to form a geometric feature point set; the geometric feature points include corner points, sharp corners, and protrusion intersection points;

[0028] Step S23, regular surface sampling of planar regions: for the flat regions in the surface of the space target, an equidistant uniform sampling method is used to generate a regular surface sampling point set in the flat region.

[0029] Further, step S21 includes:

[0030] Let the three-dimensional mesh model of the space target be a triangular mesh M=(V, F), where V={v1, v2,..., v N} is the vertex set; F={f1, f2,..., f M} is the triangular facet set; each pair of adjacent facets f i ~ f j share an edge e ij =(v p , v q ), and the corresponding unit normal vectors are n i and n j , respectively;

[0031] Based on the normal vector angle between adjacent facets in the triangular mesh model, it is judged whether there is a significant structural mutation, and edge detection is performed; when the normal angle of two adjacent facets is greater than a preset normal angle threshold and the edge length is greater than a minimum effective length, the corresponding boundary is determined as a structural edge;

[0032] The normal angle between adjacent facets is defined as:

[0033] θ ij = cos -1 (n i · n j )

[0034] The length of the edge is defined as:

[0035] L ij = ||v p -v q ||

[0036] If the following conditions are met, the edge is determined to be a characteristic edge:

[0037] θ ij > θ th , L ij > L min ,

[0038] where θ th is the normal angle threshold, L min is the minimum effective edge length;

[0039] For the extracted structure edges, an adaptive point assignment strategy based on edge length weight is used to assign the number of sampled edge points; the number of sampled points obtained by each edge is proportional to its edge length; the coordinates of the sampled points on each edge are obtained by linear interpolation;

[0040] Let the total number of edge points to be sampled be N total , the set of all characteristic edges that meet the conditions be E, and the length of each edge be L k , then the number of sampled points assigned to the kth edge is:

[0041]

[0042] Further, step S22 includes:

[0043] Let the three-dimensional model mesh be M = (V, F), where V represents the vertex set and F represents the face set. Let the set of adjacent face indices associated with each vertex v i ∈ V be F , and the corresponding normal vector set be {n j |j∈F i};

[0044] Define the curvature response of each vertex as the standard deviation of its normal distribution, denoted as:

[0045]

[0046] where n j is the dth component of the normal vector n i , and std(·) represents the standard deviation;

[0047] Selecting curvature response C from all vertices i The maximum N corner vertices as corner point sampling results, to form a corner point set

[0048] Further, step S23 includes:

[0049] For each triangular facet f i =(v1, v2, v3)∈F, whose area is:

[0050]

[0051] Let the total number of required sampling facet points be N face , then the number of sampling points allocated to the i-th facet is:

[0052]

[0053] In each facet, a regular triangle resampling algorithm is executed, a regular grid is constructed in the facet, and for any sampling point, its barycentric coordinates are:

[0054] p=αv1+βv2+γv3

[0055] Wherein, α+β+γ=1, α, β, γ≥0

[0056] The sampling points are generated by linear interpolation of the grid, and the following sampling strategy is adopted:

[0057]

[0058] Wherein, is the grid sampling resolution in the facet.

[0059] Further, step S3 specifically includes steps S31-S33:

[0060] Step S31, repeated point removal and point cloud simplification: repeated point detection and removal are performed on the sampling point cloud to eliminate coordinate redundancy caused by sampling overlap or model error;

[0061] Step S32, voxel grid division and small cluster removal: the entire point cloud space is divided into a three-dimensional grid according to a fixed voxel resolution, and the point group in each voxel is analyzed by clustering; isolated small clusters or low-density point sets in the voxel are identified and removed;

[0062] Step S33, outlier removal and statistical filtering: the point density, distance distribution, and statistical mean in the neighborhood of each point are analyzed, and outliers are removed; and combined with a statistical filtering method based on local mean and variance, further remove abnormal points that deviate significantly from the neighborhood statistical characteristics.

[0063] Further, step S5 specifically includes steps S51-S53:

[0064] Step S51, reference profile selection: perform energy statistics on the distance spectrum corresponding to all pulses, select the pulse with the strongest energy as the reference pulse; take the amplitude spectrum of the pulse as the reference distance profile;

[0065] Let the original echo matrix be:

[0066] X∈C M×N

[0067] Where M is the number of pulses, and N is the number of distance sampling points corresponding to each pulse; perform distance-wise fast Fourier transform (FFT) on each row of signals to obtain a distance spectrum matrix:

[0068] G=FFT(X,axis=2)

[0069] Calculate the energy of each pulse distance spectrum:

[0070]

[0071] Select the rth row with the maximum energy as the reference profile G ref :

[0072] G ref =|G r，： |

[0073] Step S52, cross-correlation delay estimation: perform normalized cross-correlation calculation on each pulse and the reference profile, search for the maximum cross-correlation position to estimate the relative distance drift;

[0074] For each pulse i∈[1,M], extract its distance spectrum amplitude |G i,： |, and perform normalized cross-correlation calculation with the reference profile G ref :

[0075]

[0076] Step S53, obtain the aligned distance compressed data: apply a compensation factor to the original signal matrix element by element to eliminate the distance phase drift caused by motion for each pulse, and obtain the aligned distance compressed data;

[0077] Unfold the estimated delay sequence and perform polynomial fitting:

[0078]

[0079] Define the distance-wise normalized frequency index as:

[0080]

[0081] Constructing a two-dimensional phase compensation matrix:

[0082]

[0083] Element-wise compensating the original echo signal matrix:

[0084] X′ i，n = X i，n ·φ i，n

[0085] On this basis, the distance compression data after distance compensation and phase compensation is obtained:

[0086] rngPro i，： = FFT(X′ i，： ).

[0087] Further, step S6 specifically includes steps S61-S63:

[0088] Step S61, reference distance unit selection: a quality index Q n is proposed for reference distance unit screening, and a plurality of high-quality reference distance units are selected according to the quality index as the reference for phase estimation;

[0089] Let the distance compressed data matrix be:

[0090] Y∈C M×N

[0091] where M is the number of pulses, and N is the number of distance units;

[0092] Calculate the signal-to-noise ratio SNR n , phase stability S n and average amplitude of each distance unit, and construct the quality index as follows:

[0093]

[0094] where ∠Y ：，n represents the phase sequence of the nth column, and ε is a constant. Select the first R distance units with the highest quality index {n1, n2,..., n R} as the reference distance unit set R;

[0095] Step S62, weighted phase estimation: the phase difference between the reference distance unit and the reference distance unit is unwrapped and weighted average (not equal weight average) is generated to generate the phase correction function;

[0096] Taking the first reference distance unit as the benchmark, the phase difference of other reference distance units is unwrapped, that is:

[0097] Suppose the main reference distance unit is n1, and its phase sequence is:

[0098]

[0099] For other reference distance units n r ∈R, r≥2, calculate the phase difference relative to the main reference distance unit:

[0100]

[0101] Assign weights ω r , satisfying ∑ r w r =1, and perform weighted averaging:

[0102]

[0103] Step S63, phase compensation: generate a two-dimensional phase compensation matrix according to the estimated phase error vector, and perform phase adjustment on the current pulse compression data matrix; apply a two-dimensional window function to the adjusted data to suppress sidelobe leakage;

[0104] Generate a one-dimensional phase vector compensation factor:

[0105] Φ i =exp(-j·ψ(i)), i=1,...,M

[0106] Construct a two-dimensional phase compensation matrix:

[0107] P=Φ·1 1×N

[0108] Correct the original data:

[0109] Z=Y·P

[0110] Where Z represents the corrected data.

[0111] To further suppress sidelobe leakage, perform a windowing operation on the corrected data Z:

[0112] W(m,n)=w az (m)·w rg (n)

[0113] The final corrected output data is:

[0114] Z'(m,n)=Z(m,n)·W(m,n).

[0115] The application provides a space target-oriented adaptive point cloud reconstruction and ISAR simulation imaging method, which has the following beneficial effects.

[0116] 1. The space target-oriented adaptive point cloud reconstruction and ISAR simulation imaging method of the application solves the problems of poor three-dimensional model format compatibility, low grid modeling efficiency, long imaging calculation time and the like in the prior art, supports various general model formats such as.glb,.ply and.obj, adopts adaptive point cloud reconstruction and replaces the traditional electromagnetic grid calculation with an ideal scattering point model, and thus high-efficiency, fast and accurate ISAR simulation imaging is realized.

[0117] 2. The space target-oriented adaptive point cloud reconstruction and ISAR simulation imaging method provided by the application solves the defect that the simulation image is significantly different from the real observation image due to the replacement of a complex target with a simple regular geometric body in the prior art, and proposes an adaptive sampling strategy based on geometric features. The strategy can automatically extract key scattering points with physical scattering significance from the original three-dimensional model, and mainly retains the edge, corner and surface structure features of the target, so that the structure consistency and imaging accuracy of the ISAR simulation image are effectively improved while the scattering physical authenticity is ensured.

[0118] 3. The space target-oriented adaptive point cloud reconstruction and ISAR simulation imaging method provided by the application solves the problems of limited applicability, image blurring, enhanced artifacts and easy covering of weak targets in the terahertz frequency band in the prior art mainly applicable to low-frequency radar systems, proposes an improved cross-correlation distance alignment and phase compensation method, and further introduces a multi-point weighted minimum variance phase correction algorithm. The method can effectively suppress sidelobe artifacts in the imaging process, improve the visibility of weak targets, and significantly enhance the focusing performance and image clarity of the ISAR simulation image. BRIEF DESCRIPTION OF DRAWINGS

[0119] Figure 1 It is an implementation flowchart of the space target-oriented adaptive point cloud reconstruction and ISAR simulation imaging method provided by the application.

[0120] Figure 2a It is a schematic view of a 3D model of a Jason-1 satellite jointly developed by the United States NASA and France CNES in the embodiment of the application.

[0121] Figure 2b It is a schematic view of a 3D model of a TRMM satellite jointly developed by the United States NASA and Japan JAXA in the embodiment of the application.

[0122] Figure 3a It is a schematic view of a 3D scanning model of No. 1 fragment generated by a super-high-speed impact test in the embodiment of the application.

[0123] Figure 3b is a schematic diagram of the 3D scanning model of the No. 6 fragment generated by the ultra-high speed impact test in the embodiment of the present application;

[0124] Figure 4a is a schematic diagram of the 3D point cloud of the Jason-1 satellite collected by the adaptive point cloud method in the embodiment of the present application;

[0125] Figure 4b is a schematic diagram of the 3D point cloud of the TRMM satellite collected by the adaptive point cloud method in the embodiment of the present application;

[0126] Figure 4c is a schematic diagram of the 3D point cloud of the No. 1 fragment collected by the adaptive point cloud method in the embodiment of the present application;

[0127] Figure 4d is a schematic diagram of the 3D point cloud of the No. 6 fragment collected by the adaptive point cloud method in the embodiment of the present application;

[0128] Figure 5a is a schematic diagram of the 3D point cloud of the No. 26 fragment before being optimized and processed by the adaptive point cloud method in the embodiment of the present application;

[0129] Figure 5b is a schematic diagram of the 3D point cloud of the No. 26 fragment after being optimized and processed by the adaptive point cloud method in the embodiment of the present application;

[0130] Figure 6 is a schematic diagram of the ISAR imaging of the Jason-1 satellite at different angles in the embodiment of the present application;

[0131] Figure 7 is a schematic diagram of the ISAR imaging of the TRMM satellite at different angles in the embodiment of the present application;

[0132] Figure 8 is a schematic diagram of the ISAR imaging of the No. 1 fragment at different angles in the embodiment of the present application;

[0133] Figure 9 is a schematic diagram of the ISAR imaging of the No. 6 fragment at different angles in the embodiment of the present application. DETAILED DESCRIPTION

[0134] In order to make the technical problems solved by the present application, the technical solutions adopted and the technical effects achieved more clear, the present application will be further described in detail below with reference to the drawings and embodiments. It can be understood that the specific embodiments described herein are only used to explain the present application, but not to limit the present application. In addition, it should be noted that, in order to facilitate the description, only the parts related to the present application are shown in the drawings, but not all the contents.

[0135] As Figure 1 shown, the space target-oriented adaptive point cloud reconstruction and ISAR simulation imaging method provided by the embodiment of the application includes the following processes:

[0136] Step S1, obtaining three-dimensional model data of a space target.

[0137] The space target includes spacecraft (on-orbit, retired or failed satellites, space stations, spaceships, deep space probes, etc.) and space debris (residues generated by launch missions, such as rocket upper stages, fairings, etc.). The three-dimensional model data format can include common high-precision three-dimensional formats such as.glb,.ply,.obj; and the three-dimensional model is subjected to unified scale control for subsequent point cloud sampling and imaging simulation.

[0138] Step S1 specifically includes the following steps S11-S13:

[0139] Step S11, obtaining a three-dimensional model of a spacecraft through a public data source. Some examples can be seen in Figure 2a -b schematic diagram.

[0140] In this embodiment, the three-dimensional model of the spacecraft is sourced from the NASA 3D library. The library provides spacecraft three-dimensional model data in multiple formats, including STL (Stereolithography), GLB (GL Transmission Format Binary File), and other common three-dimensional model data formats, facilitating flexible application in various scenarios such as fine modeling, simulation analysis, and imaging verification of space targets.

[0141] Step S12, obtaining a three-dimensional model of space debris through a ground hypervelocity impact test. Some examples can be seen in Figure 3a -b schematic diagram.

[0142] In this embodiment, the three-dimensional model of the space debris is derived from the ground hypervelocity impact test; the impact target is a typical spacecraft load or component, including but not limited to cameras, batteries, and other equipment. Through the hypervelocity impact test, the representative space debris generated after the impact is three-dimensionally reconstructed, and a grid model is obtained by 3D scanning to achieve real modeling of irregular targets.

[0143] Step S13, performing unified scale control on space targets of different scales. Step S13 specifically includes the following steps S131-S132:

[0144] Step S131, for a space target (such as a decommissioned satellite) with a size greater than the size threshold range (a relatively large size), downsizing processing is performed while keeping the structural proportions consistent.

[0145] The moderate downsizing processing can reduce the simulation burden caused by the over-dense grid.

[0146] Step S132, for a space target (such as an impact fragment) with a size smaller than the size threshold range (a relatively small size), isometric magnification is performed.

[0147] The isometric magnification can avoid the problem that the target features in the image are too small to be distinguished under high-resolution simulation conditions.

[0148] Step S2, an adaptive point cloud regular sampling strategy based on geometric features is adopted to extract key structural points with physical scattering significance from the surface of the three-dimensional model.

[0149] The key structural points include edge points, geometric feature points, and regular face sampling points. This step performs adaptive point cloud regular sampling, and some examples can be seen in Figure 4a - a schematic diagram. Step S2 specifically includes the following steps S21-S23:

[0150] Step S21, edge feature extraction and edge point sampling: based on the normal gradient or curvature change of the surface of the three-dimensional model, edge detection is performed to identify structural edge regions, and high-density sampling is performed in these structural edge regions.

[0151] Specifically, let the three-dimensional grid model of the space target be a triangular mesh M=(V, F), where V={v1, v2,..., v N} is the vertex set; F={f1, f2,..., f M} is the triangular face set; each pair of adjacent faces f i ~ f j share an edge e ij =(v p , v q ), and the corresponding unit normal vectors are n i and n j .

[0152] Step S211, based on the normal vector angle between adjacent faces in the triangular mesh model, it is judged whether there is a significant structural mutation, and edge detection is performed. Specifically, when the normal angle of two adjacent faces is greater than a preset normal angle threshold and the edge length is greater than a minimum effective length, the corresponding boundary is determined as a structural edge.

[0153] The normal angle between adjacent faces is defined as:

[0154] θ ij =cos-1 (n i , n j )

[0155] The length of the edge is defined as:

[0156] L ij =||v p -v q ||

[0157] If the following conditions are met, the edge is determined to be a characteristic edge:

[0158] θ ij >θ th , L ij >L min ,

[0159] wherein θ th is a normal angle threshold (set to 45° in this embodiment), and L min is a minimum effective edge length.

[0160] In step S212, for the extracted structure edges, an adaptive point assignment strategy based on edge length weight is used to assign the number of sampled edge points. Specifically, the number of sampled points obtained by each edge is proportional to its edge length, preventing the imbalance of boundary point distribution density. The coordinates of the sampled points on each edge are obtained by linear interpolation.

[0161] Let the total number of edge points to be sampled be N total , the set of all characteristic edges satisfying the conditions be E, and the length of each edge be L k , then the number of sampled points assigned to the kth edge is:

[0162]

[0163] In step S22, geometric feature point extraction and high-curvature point sampling: geometric feature points in the target three-dimensional model are detected to form a geometric feature point set. The geometric feature points include corner points, sharp corners, protrusion intersection points, and other geometric feature points with high curvature responses.

[0164] As a specific example, in the geometric feature point sampling process, a curvature estimation method based on the standard deviation of the face normal is used to identify high-curvature regions: for each vertex in the three-dimensional model, a set of associated adjacent faces is constructed; the standard deviation of the vector norm of all adjacent face normal vectors corresponding to the vertex is calculated as the local normal variation degree of the vertex; the above response value is used as a "corner scoring indicator" and is sorted from high to low; and the top several vertices with the highest curvature scores are selected as the corner sampling results.

[0165] Specifically, let the three-dimensional model mesh be M=(V, F), where V represents a vertex set, and F represents a face set. Let a set of adjacent face indexes associated with each vertex v i ∈V be a corresponding set of normal vectors be {n j |j∈F i}.

[0166] Define the curvature response of each vertex as the standard deviation of its normal distribution, denoted as:

[0167]

[0168] wherein is the dth component of the normal vector n j , and std(·) represents the standard deviation.

[0169] Select the first N i vertices with the maximum curvature response C corner from all vertices as the corner point sampling results to form a corner point set

[0170] Step S23, regular face sampling of the planar region: for the planar region (such as a satellite cabin wall, a solar panel, etc.) in the space target surface, a regular face sampling point set is generated in the planar region by using an equal-interval uniform sampling method. The regular face sampling points are uniformly distributed in the plane according to a set sampling interval, and are used to supplement the shortage of geometric feature points and high-curvature points in the planar region, so as to ensure uniform coverage and spatial resolution of the point cloud on the overall structure.

[0171] The planar region is sampled to ensure the overall surface coverage density and avoid scattering blind areas.

[0172] The application preferably adopts the following planar region regular sampling strategy: in a large smooth region of the three-dimensional model surface, in order to preserve the overall geometric distribution information, the point cloud distribution is sampled according to the area weight of each triangular face. The number of sampling points required for each face is proportional to the area proportion, so as to ensure that the point cloud density is consistent with the target structure proportion. For each triangular face, the grid filling method is used to distribute the inside into approximately uniform sampling points, so that the sampling points are uniformly distributed on the entire target outer surface.

[0173] Specifically, for each triangular face f i =(v1, v2, v3)∈F, the area is:

[0174]

[0175] Let the total number of required sampling face points be N face , and the number of sampling points allocated to the ith face be:

[0176]

[0177] A regular triangle resampling algorithm is executed within each patch, and a regular mesh is constructed within the patch. For any sampling point, the barycentric coordinates are:

[0178] p = a v1 + b v2 + g v3

[0179] where a + b + g = 1, a, b, g > 0

[0180] The sampling points are generated by linear interpolation of the mesh, and the following sampling strategy is adopted:

[0181]

[0182] wherein is the mesh sampling resolution within the patch.

[0183] Step S3, performing a structured cleaning and statistical filtering operation on the key structure points obtained by sampling in step S2. The key structure points (including edge points, geometric feature points and regular face sampling points) obtained by sampling in step S2 form point cloud data. The structured cleaning includes repeated point removal, voxel mesh clustering and small cluster removal, outlier rejection, etc. Some examples can be seen in Figure 5a - a schematic diagram. Step S3 specifically includes the following steps S31-S33:

[0184] Step S31, repeated point removal and point cloud simplification: repeated point detection and rejection is performed on the sampling point cloud to eliminate coordinate redundancy caused by sampling overlap or model error. It can ensure that each point cloud element is unique in space and avoid repeated scattering contribution.

[0185] Step S32, voxel mesh division and small cluster rejection: the entire point cloud space is divided into a three-dimensional mesh according to a fixed voxel resolution, and the point group in each voxel is analyzed for clustering; isolated small clusters or low-density point sets in the voxel are identified and removed. It can remove structure noise and reduce non-structural noise.

[0186] Step S33, outlier rejection and statistical filtering: the point density, distance distribution and statistical mean in the neighborhood of each point are analyzed, and outliers are removed; and combined with a statistical filtering method based on local mean and variance, further remove abnormal points that deviate significantly from the neighborhood statistical characteristics. It can improve the spatial continuity and smoothness of the point cloud.

[0187] The outlier is a point with an abnormally large distance or extremely low density. As a specific example, the point cloud data is subjected to outlier detection and rejection processing. First, according to the distance distribution between each point in the point cloud and other points in its neighborhood, the average neighborhood distance is calculated And set the radius threshold R in the range of 0.005m to 0.02m, the value is When a certain point in the radius R range of the number of adjacent points is less than 3 or more than 10 outside the range, or the point and the average distance of the neighborhood is more than (σ is the neighborhood distance standard deviation) is determined as an outlier and is removed.

[0188] On this basis, using statistical filtering method, with the average distance of the neighborhood as the benchmark, when the point distance deviation More than 2σ, further determine the abnormal point and delete.

[0189] Through the above method, the points with too large distance or too low neighborhood density can be effectively removed, and the spatial continuity and smoothness of the point cloud are improved.

[0190] Step S4, the point cloud after the structured cleaning and statistical filtering operation in step S3 is used as the ideal scattering center set, and the point target model is used to simulate the received echo signal.

[0191] As a specific example, for a high-resolution terahertz band radar system, assuming that the target motion has been compensated accurately, and the residual effects of velocity v and acceleration a are ignored, the echo signal can be modeled as:

[0192]

[0193] Where j represents the imaginary unit; R(t) is the instantaneous distance from the radar platform to the target geometric center; c represents the speed of light; f represents the frequency; ρ(x, y) is the scattering intensity function of the target on the imaging plane; f x (t), f y (t) are time-varying spatial frequency components, which are defined as:

[0194] f x (t) = 2fcosθ(t) / c, f y (t) = 2fsinθ(t) / c

[0195] The echo signal modeling is a conventional technical means, which is not described here.

[0196] Step S5, using the improved cross-correlation-based distance alignment and phase compensation method, distance alignment and phase compensation are performed. The distance drift between radar echo pulses can be eliminated, and the ISAR image focusing quality is improved.

[0197] In the traditional ISAR imaging method, the rough motion compensation of the target is usually realized by range alignment, specifically, the energy peak value matching or cross-correlation calculation is performed on the range profile line corresponding to each echo pulse, so that the distance drift between different pulses is estimated and aligned. However, the alignment method only processes the amplitude information, and after the execution, there is still a significant phase nonlinear error, especially near the frequency point in the range direction, which is manifested as inconsistent phase jump, thereby causing the final ISAR image to appear fuzzy, sidelobe enhancement or artifact stretching and the like.

[0198] To overcome the above defects, the improved distance alignment and phase compensation method based on envelope cross-correlation and polynomial fitting is further introduced on the basis of the traditional distance alignment. Step S5 specifically includes the following steps S51-S53:

[0199] Step S51, reference profile selection: the energy of the distance spectrum corresponding to all pulses is counted, and the pulse with the strongest energy is selected as the reference pulse; the amplitude spectrum of the pulse is taken as the reference distance profile.

[0200] Let the original echo matrix be:

[0201] X∈C M×N

[0202] Wherein, M is the pulse number, and N is the distance sampling point number corresponding to each pulse. The distance spectrum matrix is obtained by performing distance direction fast Fourier transform (FFT) on each row signal:

[0203] G=FFT(X,axis=2)

[0204] The energy of each pulse distance spectrum is calculated:

[0205]

[0206] The rth row with the maximum energy is selected as the reference profile G ref :

[0207] G ref =|G r，： |

[0208] Step S52, cross-correlation delay estimation: the normalized cross-correlation calculation is performed on each pulse and the reference profile, and the maximum cross-correlation position is searched to estimate the relative distance drift.

[0209] For each pulse i∈[1,M], the distance spectrum amplitude |G i，： | of the pulse is extracted, and the normalized cross-correlation calculation is performed on the reference profile G ref :

[0210]

[0211] Step S53, obtaining the aligned range-compressed data: applying the compensation factor to the original signal matrix element by element to eliminate the distance phase drift caused by motion for each pulse to obtain the aligned range-compressed data.

[0212] The estimated delay sequence Unfold and perform polynomial fitting:

[0213]

[0214] Define the range-to-normalized frequency index as:

[0215]

[0216] Construct a two-dimensional phase compensation matrix:

[0217]

[0218] Element-wise compensation is performed on the original echo signal matrix:

[0219] X′ i，n =X i，n ·φ i，n

[0220] On this basis, the range-compressed data after range compensation and phase compensation is obtained:

[0221] rngPro i，： =FFT(X′ i，： )

[0222] The method not only realizes the synchronization alignment of the radar echo in the amplitude domain, but also introduces a phase domain compensation mechanism, which fundamentally solves the imaging degradation problem of the traditional method in the nonlinear phase jump area.

[0223] Step S6, considering the signal-to-noise ratio, phase stability and multi-point weighting mechanism, an improved minimum variance phase correction method is used for motion compensation and phase correction. This step can realize accurate compensation of residual phase error of high dynamic target, and further improve the clarity and focusing degree of ISAR image. Step S6 specifically includes the following steps S61-S63:

[0224] Step S61, reference distance unit selection: a quality index Q n is proposed for reference distance unit screening, and a plurality of high-quality reference distance units are selected as the reference for phase estimation according to the quality index.

[0225] Let the range-compressed data matrix be:

[0226] Y∈C M×N

[0227] where M is the number of pulses, and N is the number of distance cells;

[0228] Calculate the signal-to-noise ratio SNR of each distance cell n , phase stability S n , and average amplitude The quality index is constructed as follows:

[0229]

[0230] where ∠Y ：，n represents the phase sequence of the nth column, and ε is a constant. From among them, the first R distance cells {n1, n2,..., n R} with the highest quality index are selected as the reference distance cell set R.

[0231] Step S62, weighted phase estimation: the phase difference between the reference distance cells and the reference distance cells is unwrapped, and a phase correction function is generated by weighted average (not equal weight average). Specifically:

[0232] Take the first reference distance cell as the reference, and unwrap the phase difference of other reference distance cells, that is,

[0233] Let the main reference distance cell be n1, and its phase sequence be:

[0234]

[0235] For other reference distance cells n r ∈R, r≥2, calculate the phase difference relative to the main reference distance cell:

[0236]

[0237] Assign a weight ω r , satisfying∑ r w r =1, and perform weighted average:

[0238]

[0239] Step S63, phase compensation: generate a two-dimensional phase compensation matrix according to the estimated phase error vector, and perform phase adjustment on the current pulse compression data matrix; apply a two-dimensional window function to the adjusted data to suppress sidelobe leakage. The specific calculation is as follows:

[0240] Generate a one-dimensional phase vector compensation factor:

[0241] Φ i =exp(-j·ψ(i)),i=1,...,M

[0242] Constructing a two-dimensional phase compensation matrix:

[0243] P = Φ · 1 1×N

[0244] Correcting the original data:

[0245] Z = Y · P

[0246] wherein Z represents the corrected data.

[0247] To further suppress sidelobe leakage, a windowing operation (Hamming window) is performed on the corrected data Z:

[0248] W(m, n) = w az (m) · w rg (n)

[0249] The final corrected output data is:

[0250] Z'(m, n) = Z(m, n) · W(m, n)

[0251] Step S7, performing a two-dimensional Fourier transform on the data after motion compensation and phase correction, and outputting the final two-dimensional ISAR image.

[0252] In one specific embodiment, to verify the effectiveness of the method described in the present application, a simulation experiment was carried out. As shown in FIG. 4, the space targets selected in the experiment were two retired satellites Jason-1 and TRMM, and two fragments No. 1 and No. 6. After processing the three-dimensional models, 39770|88295|23325|23887 point cloud data with scattering characteristics were generated, respectively, as an ideal scattering center set for simulation imaging.

[0253] Table 1 shows the specific parameters of the ISAR imaging process in the present example, Figure 6 shows the ISAR imaging results of Jason-1 under multi-angle observation obtained by the ISAR simulation imaging method described in the present application. Figure 7 shows the ISAR imaging results of TRMM under multi-angle observation obtained by the ISAR simulation imaging method described in the present application. Figure 8 shows the ISAR imaging results of fragment No. 1 under multi-angle observation obtained by the ISAR simulation imaging method described in the present application. Figure 9 shows the ISAR imaging results of fragment No. 6 under multi-angle observation obtained by the ISAR simulation imaging method described in the present application.

[0254] Table 1 ISAR imaging parameter settings

[0255]

[0256]

[0257] It can be seen that the method effectively captures the key scattering characteristics of the spatial target through the adaptive sampling strategy based on geometric characteristics, and combines the improved distance alignment and phase compensation algorithm and the multi-point weighted minimum variance phase correction technology, significantly improves the focusing degree and structural consistency of the ISAR simulation image, effectively suppresses the image blur and sidelobe artifact problems, improves the imaging quality and weak target detection ability at the terahertz band, and verifies the effectiveness and practicability of the method.

[0258] Finally, it should be noted that: the above embodiments are only used to illustrate the technical solutions of the present application, and not to limit them; although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that: the modification of the technical solutions recorded in the foregoing embodiments, or the equivalent replacement of part or all of the technical features, does not make the essence of the corresponding technical solution deviate from the scope of the technical solutions of the embodiments of the present application.

Claims

1. An adaptive point cloud reconstruction and ISAR simulation imaging method for space targets, characterized in that, include: Step S1: Obtain the three-dimensional model data of the space target; Step S2: Adopt an adaptive point cloud rule sampling strategy based on geometric features to extract key structural points with physical scattering significance from the surface of the 3D model; Step S3: Perform structured cleaning and statistical filtering operations on the key structural points obtained in step S2; Step S4: Using the point cloud after the structured cleaning and statistical filtering operations in step S3 as an ideal scattering center set, a point target model is used to simulate the received echo signal. Step S5: Perform distance alignment and phase compensation using a distance alignment and phase compensation method based on improved cross-correlation. Step S6: Considering the signal-to-noise ratio, phase stability, and multi-point weighting mechanism, an improved minimum variance phase correction method is adopted to perform motion compensation and phase correction. Step S6 specifically includes the following steps S61-S63: Step S61, Reference distance cell selection: A quality index for selecting reference distance cells is proposed. Multiple high-quality reference distance cells are selected based on quality indicators as the benchmark for phase estimation. Let the distance-compressed data matrix be: ； in, The number of pulses. This represents the number of distance units. Calculate the signal-to-noise ratio of each distance cell. Phase stability ( and average amplitude The quality indicators are constructed as follows: ； in, Indicates the first The phase sequence of the column, The value is a constant; select the top values ​​with the highest quality indicators. distance units As a set of reference distance cells ; Step S62, Weighted Phase Estimation: Unwrap the phase difference between the reference range cell and the baseline range cell, and generate a phase correction function by weighted averaging. Using the first reference range cell as a reference, the phase difference of other reference range cells is unwrapped, i.e.: Let the master reference distance cell be Its phase sequence is: ； For other reference distance cells , Calculate its phase difference relative to the master reference range cell: ； Assign weights ,satisfy Perform a weighted average: ； Step S63, Phase Compensation: Generate a two-dimensional phase compensation matrix based on the estimated phase error vector, and perform phase adjustment on the current pulse compressed data matrix; apply a two-dimensional window function to the adjusted data to suppress sidelobe leakage; Generate a one-dimensional phase vector compensation factor: ； Construct a two-dimensional phase compensation matrix: ； Correct the original data: ； Where Z represents the corrected data; To further suppress sidelobe leakage, a windowing operation is applied to the corrected data Z: ； The final corrected output data is as follows: ； Step S7: Perform a two-dimensional Fourier transform on the motion-compensated and phase-corrected data to output the final two-dimensional ISAR image.

2. The adaptive point cloud reconstruction and ISAR simulation imaging method for space targets according to claim 1, characterized in that, Step S1 includes the following steps S11-S13: Step S11: Obtain a 3D model of the spacecraft from publicly available data sources; Step S12: Obtain a three-dimensional model of the space debris through a ground-based ultra-high-speed impact test; Step S13: Perform unified scale control on spatial targets of different scales.

3. The adaptive point cloud reconstruction and ISAR simulation imaging method for space targets according to claim 2, characterized in that, Step S13 Includes the following steps S131-S132: Step S131: For spatial targets whose size is larger than the size threshold range, reduce their size while maintaining the structural proportions. Step S132: For spatial targets whose size is smaller than the size threshold range, scale them up proportionally.

4. The adaptive point cloud reconstruction and ISAR simulation imaging method for space targets according to claim 1, characterized in that, Step S2 includes the following steps S21-S23: Step S21, Edge Feature Extraction and Edge Point Sampling: Based on the normal gradient or curvature change of the 3D model surface, edge detection is performed to identify structural edge regions, and high-density sampling is performed in these structural edge regions; Step S22, Geometric feature point extraction and high curvature point sampling: Detect geometric feature points in the target 3D model to form a set of geometric feature points; the geometric feature points include corner points, sharp corners, and protruding intersections; Step S23: Regular surface sampling of planar regions: For flat regions on the surface of spatial targets, a set of regular surface sampling points is generated within the flat region using a uniform sampling method with equal spacing.

5. The adaptive point cloud reconstruction and ISAR simulation imaging method for space targets according to claim 4, characterized in that, Step S21 includes: Assume the three-dimensional mesh model of the spatial target is a triangular mesh. ,in, For vertex set; It is a set of triangular facets; each pair of adjacent facets Shared edge The corresponding unit normal vectors are respectively and ; The presence of obvious structural abrupt changes is determined based on the angle between the normal vectors of adjacent faces in the triangular mesh model, and edge detection is performed. When the angle between the normal vectors of two adjacent faces is greater than a preset threshold and the side length of the common edge of the two adjacent faces is greater than the minimum effective length, the corresponding boundary is determined as a structural edge. Define the angle between the normals of adjacent facets as: ； Define the length of this side as: ； An edge is considered a feature edge if the following conditions are met: ， in, It is the threshold of the included normal angle. It is the minimum effective side length; For the extracted structural edges, an adaptive point allocation strategy based on edge length weights is used to allocate the number of sampling edge points; the number of sampling points obtained for each edge is proportional to its edge length; the coordinates of the sampling points on each edge are obtained through linear interpolation. Let the total number of edge points to be sampled be... The set of all feature edges that satisfy the conditions is The length of each side is Then the first The number of sampling points assigned to each edge is: 。 6. The adaptive point cloud reconstruction and ISAR simulation imaging method for space targets according to claim 4, characterized in that, Step S22 includes: Let the 3D model mesh be ,in, Represents the set of vertices. Let each vertex represent a set of facets. The associated set of adjacent face indices is The corresponding set of normal vectors is ; The curvature response at each vertex is defined as the sum of the standard deviations of its normal distribution, denoted as: ； in, Normal vector The One portion, Indicates standard deviation; Select curvature response from all vertices The largest front Each vertex is used as a corner sampling result to form a corner set. .

7. The adaptive point cloud reconstruction and ISAR simulation imaging method for space targets according to claim 4, characterized in that, Step S23 includes: For each triangular facet Its area is: ； Let the total number of sampling points be... Then the first The number of sampling points allocated to each patch is: ； Within each facet, a regular triangle resampling algorithm is performed to construct a regular mesh within the facet. For any sampling point, its centroid coordinates are: ； in, , ; Sampling points are generated through grid linear interpolation, using the following sampling strategy: ； in, , where is the grid sampling resolution within the patch.

8. The adaptive point cloud reconstruction and ISAR simulation imaging method for space targets according to claim 1, characterized in that, Step S3 specifically includes the following steps S31-S33: Step S31, Repetitive Point Removal and Point Cloud Simplification: Repeated points are detected and removed from the sampled point cloud to eliminate coordinate redundancy caused by sampling overlap or model error; Step S32, Voxel Mesh Generation and Small Cluster Removal: Divide the entire point cloud space into a three-dimensional mesh at a fixed voxel resolution, and perform cluster analysis on the point groups within each voxel; identify and remove isolated small clusters or low-density point sets in the voxels; Step S33, Outlier Removal and Statistical Filtering: Analyze the density, distance distribution, and statistical mean of points in the neighborhood of each point to remove outliers; and combine statistical filtering methods based on local mean and variance to further remove abnormal points that deviate significantly from the statistical characteristics of the neighborhood.

9. The adaptive point cloud reconstruction and ISAR simulation imaging method for space targets according to claim 8, characterized in that, Step S5 specifically includes the following steps S51-S53: Step S51, Reference profile selection: Perform energy statistics on the distance spectrum corresponding to all pulses, and select the pulse with the strongest energy as the reference pulse; The amplitude spectrum of the pulse is used as a reference distance profile; Let the original echo matrix be: ； in, For the number of pulses, The number of distance sampling points corresponds to each pulse; a range-to-fast Fourier transform (FFT) is performed on each line of signal to obtain the range spectrum matrix: ； Calculate the energy of the distance spectrum for each pulse: ； Select the one with the highest energy Line as reference outline : ； Step S52, Cross-correlation delay estimation: Perform normalized cross-correlation calculation between each pulse and the reference profile, search for the position of maximum cross-correlation, and estimate its relative distance drift. For each pulse Extract its distance spectrum amplitude and with reference contour Perform normalized cross-correlation calculation: ； Step S53: Obtain aligned distance compressed data: Apply a compensation factor to each element of the original signal matrix to eliminate the distance phase drift caused by motion in each pulse, and obtain aligned distance compressed data; For the estimated delay sequence Expand and perform polynomial fitting: ； Define the distance-oriented normalized frequency index as: ； Construct a two-dimensional phase compensation matrix: ； Element-by-element compensation is performed on the original echo signal matrix: ； Based on this, the range compression data after range compensation and phase compensation is obtained: 。

Images