A Method for Obtaining Wideband RCS of Electrically Large Targets Based on ACA and Improved AWE
By improving the AWE and ACA methods and combining the pseudospectral derivative method and the octree algorithm, the problems of low efficiency and high memory requirements in acquiring the broadband RCS of electrically large targets are solved, achieving efficient computation and memory saving.
Patent Information
- Application Number
- CN202411157636.7
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2024-08-22
- Publication Date
- 2025-11-11
- Estimated Expiration
- 2044-08-22
AI Technical Summary
Existing technologies are inefficient and require high computer memory when acquiring broadband radar cross sections of electrically large targets. The method of moments (MoM) is complex to calculate and the filling of high-order impedance matrices results in long computation time.
An improved asymptotic waveform estimation method (AWE) combined with an adaptive cross approximation method (ACA) is adopted. The impedance matrix is decomposed by the pseudospectral derivative method (PSDM) and the octree algorithm. The low-rank matrix is used to accelerate the calculation of the induced current vector, avoiding the direct filling of the high-order impedance matrix.
It effectively improves the efficiency of radar cross section acquisition, reduces computer memory requirements, simplifies the calculation process, and increases calculation speed.
Smart Images

Figure CN119026363B_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of electromagnetic simulation technology and relates to a method for obtaining the broadband RCS of electrically large targets, specifically a method for obtaining the broadband RCS of electrically large targets based on ACA and an improved AWE. Background Technology
[0002] Radar cross section (RCS) is a physical quantity that quantitatively characterizes the intensity of a target's scattering of incident waves. With the development of science and technology, radar detection systems and stealth technologies are constantly advancing, which places higher demands on the rapid and accurate acquisition of RCS. The Method of Moments (MoM), as one of the most commonly used numerical methods for electromagnetic field calculation, can accurately analyze the RCS of a target. However, electrically large targets are usually much larger than the wavelength of electromagnetic waves, requiring consideration of more physical effects and computational complexity when dealing with their electromagnetic properties. MoM, due to its large memory requirements and long computation time, is difficult to solve the electromagnetic problems of electrically large targets. For broadband problems, MoM requires repeated filling of the impedance matrix at each frequency point, and because the MoM generates a dense matrix, it wastes a lot of computation time during solution. Therefore, researching broadband RCS acquisition methods for electrically large targets is of great significance.
[0003] Xi'an University of Electronic Science and Technology, in its patent application "A Method for Obtaining Broadband RCS of Periodic Structures Based on CBFM and AWE" (application date: May 15, 2023, application number: 202310540021.1, publication number: CN 116484642A), proposed a method for obtaining broadband RCS of periodic structures based on CBFM and AWE. The implementation steps are: initializing parameters; calculating the excitation vector of each sub-region; calculating the induced current of each sub-region at each frequency corresponding to the wavenumber; and obtaining the broadband RCS of the periodic structure. This invention divides the entire target into many sub-regions, constructs basis functions on each sub-region using the characteristic basis function method CBFM, and reflects the interaction between sub-regions through secondary characteristic basis functions, thereby reducing the order of the impedance matrix and avoiding the defect of excessive impedance matrix due to too many unknowns, thus reducing the complexity of solving the impedance matrix. At the same time, it uses asymptotic waveform estimation AWE to calculate the current at other frequencies based on the induced current at the center frequency, avoiding the repeated filling of the impedance matrix at each frequency, effectively improving the acquisition efficiency. However, this invention requires iterative calculation of a large number of matrix-vector products, which affects the further improvement of the acquisition efficiency. In addition, when using asymptotic waveforms to estimate AWE, this invention also needs to fill a large number of high-order impedance matrices, which increases the requirements for computer memory. Summary of the Invention
[0004] The purpose of this invention is to overcome the shortcomings of the prior art and propose a method for acquiring the wideband RCS of electrically large targets based on ACA and improved AWE, which solves the technical problems of low acquisition efficiency and high computer memory requirements in the prior art.
[0005] To achieve the above objectives, the technical solution adopted by the present invention includes the following steps:
[0006] (1) Initialize parameters:
[0007] The electric field of a broadband uniform plane wave incident on the surface of an electrically large target is initialized as E. in The frequency band f includes frequencies starting from f1 and ending at f1. B Let B be the cutoff frequencies, and let f be the frequency of each frequency. b The corresponding wavelengths and wavenumbers are λ. b k b The center frequency is f0; the induced current vector I(k) on the electrically large target surface in the AWE asymptotic waveform estimation method. b The truncation order of the Taylor series of ) is Q; the RWG basis function of each of the W triangular facet pairs formed by partitioning the electrically large target surface is f. w (r), where B≥4, and f0 is the th (r) when B is even. The frequency, when B is odd, f0 is the frequency. There are frequencies, and r is the field point position vector;
[0008] (2) Sampling the incident wave on the surface of the electrically large target:
[0009] Based on the GLC distribution, the incident wave on the surface of an electrically large target is sampled N times with k0 as the center, resulting in N GLC sampling points, where N>(Q+1), and the nth GLC sampling point is σ. n , σ n ∈(k0-Δk,k0+Δk), Δk is the radius of wavenumber variation centered at k0;
[0010] (3) Construct the impedance matrix at each GLC sampling point:
[0011] Based on the octree algorithm, the RWG basis functions of W triangle facet pairs are divided into S groups Ω={Ω1,Ω2,…,Ω…} s ,…,Ω S}, and based on the method of moments (MoM) for σ n Two Ω points with a spatial distance less than λ1 μ and Ω ν The near-field impedance matrix formed between them is filled to obtain σ. nX near-field blocks at location X, and simultaneously based on the adaptive cross-approximation method ACA for σ n Two Ω points with a spatial distance greater than λ1 χ and Ω ξ The far-field impedance matrix formed between them is compressed and decomposed to obtain σ. n Y far-field blocks at location Y, then through σ n Each far-field block and near-field blocks Calculate σ n The impedance matrix Z(σ) at the point n ), where χ≠ξ, 1≤μ,ν,χ,ξ≤S and are integers, when μ≠ν The near-field block is a mutual impedance block, and when μ = ν it is a self-impedance near-field block;
[0012] (4) Obtain the electrically large target broadband RCS:
[0013] The asymptotic waveform estimation method AWE based on the pseudospectral derivative method PSDM, and through Z(σ n ) Calculate the induced current vector I(k) on the electrically large target surface. b The w-th element I w (k b ), and then through I w (k b ), f w (r) and E in Calculate the radar cross section (RCS) at each wavenumber kb. b ), to obtain the electrically large target bandwidth RCS = {RCS(k1), RCS(k2), ..., RCS(k b ),...,RCS(k B )}.
[0014] Compared with the prior art, the present invention has the following advantages:
[0015] (1) When the present invention uses the improved progressive waveform estimation method to calculate the induced current vector of electrically large target surface at each wave number, it uses the pseudo-spectral differential matrix and impedance matrix to represent the high-order impedance matrix of each order, which avoids the defect of directly filling the high-order impedance matrix with high complexity in the prior art and effectively reduces the demand for computer memory space.
[0016] (2) Based on the adaptive cross approximation method, the present invention compresses and decomposes the far-field impedance matrix into the product of two low-rank matrices, and uses the low-rank matrix to accelerate the calculation of a large number of matrix-vector products in the asymptotic waveform estimation method. This avoids the defect of complex calculation process caused by directly calculating a large number of matrix-vector products in the existing asymptotic waveform estimation method, and effectively improves the efficiency of radar cross section acquisition. Attached Figure Description
[0017] Figure 1 This is a flowchart illustrating the implementation of the present invention. Detailed Implementation
[0018] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.
[0019] Reference Figure 1 The present invention includes the following steps:
[0020] Step 1) Initialize parameters:
[0021] The electric field of a broadband uniform plane wave incident on the surface of an electrically large target is initialized as E. in The frequency band f includes frequencies starting from f1 and ending at f1. B Let B be the cutoff frequencies, and let f be the frequency of each frequency. b The corresponding wavelengths and wavenumbers are λ. b k b The center frequency is f0; the induced current vector I(k) on the electrically large target surface in the AWE asymptotic waveform estimation method. b The truncation order of the Taylor series of ) is Q; the RWG basis function of each of the W triangular facet pairs formed by partitioning the electrically large target surface is f. w (r), where B≥4, and f0 is the th (r) when B is even. The frequency, when B is odd, f0 is the frequency. There are frequencies, and r is the field point position vector; in this embodiment, f1 = 1 GHz, f b =1.4GHz, B=41, f0 is the 21st frequency, f0=1.2GHz, W=18175.
[0022] The method for subdividing the electrically large target surface is as follows: the electrically large target surface is divided into H triangular patches, and each triangular patch and the triangular patches that share a common edge with it are combined into a pair of triangular patches to obtain W pairs of triangular patches, where H ≥ 4 and is an even number; in this embodiment, H = 12650.
[0023] The side length of each triangular facet obtained by subdivision does not exceed Furthermore, every two triangular faces either do not intersect or share only one common edge.
[0024] RWG basis functions f w The expression for (r) is:
[0025]
[0026] Among them, l w It represents two triangular facets and The length of the common side of the w-th triangular facet pair. and They represent and area, express The vector pointing from the vertex to the field point. express The vector pointing from the field point to the vertex.
[0027] Step 2) Sample the incident wave on the electrically large target surface:
[0028] Based on the GLC distribution, the incident wave on the surface of an electrically large target is sampled N times with k0 as the center, resulting in N GLC sampling points, where N>(Q+1), and the nth GLC sampling point is σ. n , σ n ∈(k0-Δk,k0+Δk), Δk is the radius of wavenumber variation centered at k0.
[0029] GLC sampling points σ n The calculation formula is:
[0030]
[0031] Where π represents the value of a circle's circumference.
[0032] Step 3) Construct the impedance matrix at each GLC sampling point:
[0033] Before constructing the impedance matrix, the RWG basis functions need to be grouped.
[0034] Based on the octree algorithm, the RWG basis functions of W triangle facet pairs are divided into S groups Ω={Ω1,Ω2,…,Ω…} s ,...,Ω S The implementation method is as follows:
[0035] Based on the octree algorithm, objects containing electrically large targets are grouped into (x) min ,y min ,z min Starting from θ, the cube with side length θ is evenly divided into eight equal parts. Each of these eight parts is then evenly divided into eight equal parts again, and so on, until the side length of each resulting cube is less than θ. Then, select S cubes from all the cubes obtained from the last split that contain RWG basis functions. Next, group the RWG basis functions in each of these cubes into a group, resulting in the RWG basis function group Ω = {Ω1, Ω2, ..., Ω...}. s ,…,ΩS}. The formula for calculating θ is:
[0036] θ=max{|x max -x min |,|y max -y min |,|z max -z min |}
[0037] Where x max x min y max y min z max , z min Let max{·} represent the maximum and minimum values of the x, y, and z axis components of the coordinates of the vertices of the H triangular facets, respectively. max{·} represents the maximum value operation, and |·| represents the modulus operation.
[0038] After grouping the RWG basis functions, the impedance matrix formed between the two groups of RWG basis functions needs to be divided into a near-field impedance matrix and a far-field impedance matrix based on the spatial distance between them. Then, at each GLC sampling point σ n At this point, σ is obtained by filling the near-field impedance matrix based on the method of moments (MoM). n The near-field block at the location is used to compress and decompose the far-field impedance matrix using the adaptive cross-approximation method (ACA) to obtain σ. n The far-field block at the location, and finally the obtained σ n σ is calculated from all near-field and far-field blocks at that location. n The impedance matrix at that point.
[0039] MoM based on σ n Two Ω points with a spatial distance less than λ1 μ and Ω ν The near-field impedance matrix formed between them is filled to obtain σ. n The x-th near-field block is... Where 1≤μ, ν≤S and are integers, when μ≠ν It is a near-field block with mutual impedance, and when μ = ν it is a near-field block with self-impedance.
[0040] Based on the adaptive cross-approximation method ACA, σ n Two Ω points with a spatial distance greater than λ1 χ and Ω ξ The far-field impedance matrix formed between them is compressed and decomposed to obtain σ. n The y-th far-field block is located at Y far-field blocks. The specific steps are as follows:
[0041] First, set up two low-rank matrices. calculate Corresponding approximate far-field block The calculation formula is:
[0042]
[0043] Then, calculate and Error matrix R y (σ n And ensure that the following conditions are met:
[0044]
[0045] Finally, using replace
[0046]
[0047] Where 1≤χ, ξ≤S and are integers, χ≠ξ; express The corresponding dimension is u y ×v y Approximate far-field block, and They represent The two corresponding low-rank matrices, ι y express effective rank, R y (σ n )express and The error matrix is denoted by ||·||, where ||·|| represents the Frobenus norm of the matrix, and ε represents the given error iteration threshold.
[0048] Through σ n Each near-field block at the location and far-field blocks Calculate σ n The impedance matrix Z(σ) at the point n ):
[0049]
[0050] Step 4) Obtain the electrically large target bandwidth RCS:
[0051] The key to obtaining the broadband RCS of electrically large targets is to obtain the induced current vector I(k) on the surface of the electrically large target. b ).
[0052] The asymptotic waveform estimation method AWE based on the pseudospectral derivative method PSDM, and through Z(σ n) Calculate the induced current vector I(k) on the electrically large target surface. b The w-th element I w (k b The specific steps are as follows:
[0053] 4a) Based on MoM, and through f w (r) and E in Calculate the excitation vector Γ(k0) at k0 and its κ-th derivative Γ. (κ) The w-th element Γ of (k0) w (k0), Γ w (κ) (k0):
[0054] Γ w (k0)=∫f w (r)·E in dr
[0055] Γ w (κ) (k0)=∫E in(κ) ·f w (r)dr
[0056] Among them, E in(κ) For E in The κ-th derivative, κ≥1.
[0057] 4b) Calculate the element D in the α-th row and β-th column of the N×N dimensional pseudospectral differential matrix D based on the pseudospectral derivative method (PSDM). αβ :
[0058]
[0059] Where, σ α σ β For the αth and βth GLC sampling points σ α σ β The corresponding sampling points after coordinate transformation The coefficients are related to the values of α and β, where 1 ≤ α, β ≤ N.
[0060] 4c) Calculate the Taylor expansion coefficient vectors of order 0 and q using the results of 4b) and 4c). q :
[0061] m0 = Z -1 (k0)·Γ(k0)
[0062]
[0063]
[0064] Where Z(k0) is the impedance matrix at the center wavenumber k0, Z -1 (k0) is the inverse matrix of Z(k0), Z (i) (k0) is the i-th order higher-order impedance matrix at the center wavenumber k0, D i Let i be the power of matrix D. For matrix D i The The elements of a row, i≥1, 1≤q≤Q.
[0065] In step 3), σ has been... n Impedance matrix Z(σ) n The far-field block portion of the matrix is represented by two low-rank matrices based on ACA. The product of the products can be solved more quickly in 4c) by using a low-rank matrix. q The matrix-vector product Z(σ) at time n )m q-i This increases the computational speed, thereby improving the efficiency of radar cross section (RCS) acquisition.
[0066] Solve for m in 4c) q During the process, the pseudospectral differential matrix D and σ are constructed using AWE based on the pseudospectral derivative method PSDM. n Impedance matrix Z(σ) n The expression represents the higher-order impedance matrix Z of order i, which has higher complexity. (i) (k0), avoiding direct filling of Z (i) (k0), thereby reducing the need for computer memory space.
[0067] Obtain m0 and m through the above steps. q After that, I can be obtained through the Taylor expansion formula. w (k b However, due to the bandwidth limitation of Taylor expansion, the calculation results may deviate significantly outside the bandwidth. In order to obtain a larger bandwidth, Padé approximation is introduced in 4d) to transform the Taylor series into a rational function.
[0068] 4d) Through m0, m q The w-th element m w,0 m w,q Calculate each k b The induced current vector I(k) at the location b The w-th element I w (k b ):
[0069]
[0070]
[0071] L+O=Q
[0072] Among them, L, a w,l They are polynomials The number of times and corresponding to I w (k b The coefficients of the l-th term, O, d w,o They are polynomials The number of times and corresponding to I w (k b The coefficient of the o-th term of Q is L, O, where L and O are positive integers. If Q is even, L = O; if Q is odd, the absolute value of the difference between L and O is |LO| = 1.
[0073] Then through I w (k b ), f w (r) and E in Calculate each wave number k b Radar cross section RCS(k) at the location b ), thus obtaining the electrically large target bandwidth RCS = {RCS(k1), RCS(k2), ..., RCS(k... b ),...RCS(k B )},RCS(k b The formula for calculating ) is:
[0074]
[0075] Where R represents the distance between the electrically large target surface and the radar. This indicates that R approaches infinity, j represents the imaginary unit, and η represents the wave impedance in air. The unit vector representing the observation point position vector r', Λ represents the outer surface of the electrically large target, ∫∫ Λ This indicates that the surface integral operation is performed on Λ, e represents the natural constant, and |·| represents the modulo operation.
Claims
1. A method for acquiring the broadband RCS of electrically large targets based on ACA and an improved AWE, characterized in that, An improvement to the Asymptotic Waveform Estimation (AWE) method is achieved by introducing the pseudospectral derivative method (PSDM), specifically including the following steps: (1) Initialize parameters: The electric field of a broadband uniform plane wave incident on the surface of an electrically large target is initialized as E. in The frequency band f includes frequencies starting from f1 and ending at f1. B Let B be the cutoff frequencies, and let f be the frequency of each frequency. b The corresponding wavelengths and wavenumbers are λ. b k b The center frequency is f0; the induced current vector I(k) on the electrically large target surface in the AWE asymptotic waveform estimation method. b The truncation order of the Taylor series of ) is Q; The RWG basis function of each of the W pairs of triangular facets formed by partitioning the electrically large target surface is f. w (r), where B≥4, and f0 is the th (r) when B is even. The frequency, when B is odd, f0 is the frequency. There are frequencies, and r is the field point position vector; (2) Sampling the incident wave on the surface of the electrically large target: Based on the GLC distribution, the incident wave on the surface of an electrically large target is sampled N times with k0 as the center, resulting in N GLC sampling points, where N>(Q+1), and the nth GLC sampling point is σ. n , σ n ∈(k0-Δk,k0+Δk), Δk is the radius of wavenumber variation centered at k0; (3) Construct the impedance matrix at each GLC sampling point: Based on the octree algorithm, the RWG basis functions of W triangle facet pairs are divided into S groups Ω={Ω1,Ω2,…,Ω…} s ,…,Ω S }, and based on the method of moments (MoM) for σ n Two Ω points with a spatial distance less than λ1 μ and Ω v The near-field impedance matrix formed between them is filled to obtain σ. n X near-field blocks at location X, and simultaneously based on the adaptive cross-approximation method ACA for σ n Two Ω points with a spatial distance greater than λ1 χ and Ω ξ The far-field impedance matrix formed between them is compressed and decomposed to obtain σ. n Y far-field blocks at location Y, then through σ n Each far-field block and near-field blocks Calculate σ n The impedance matrix Z(σ) at the point n ), where χ≠ξ, 1≤μ,v,χ,ξ≤S and are integers, when μ≠v The near-field block is a mutual impedance block, and the near-field block is a self-impedance block when μ = v. (4) Obtain the electrically large target broadband RCS: The asymptotic waveform estimation method AWE based on the pseudospectral derivative method PSDM, and through Z(σ n ) Calculate the induced current vector I(k) on the electrically large target surface. b The w-th element I w (k b ), and then through I w (k b ), f w (r) and E in Calculate each wave number k b Radar cross section RCS(k) at the location b ), to obtain the electrically large target bandwidth RCS = {RCS(k1), RCS(k2), ..., RCS(k b ),...RCS(k B )}.
2. The method according to claim 1, characterized in that, The method for dividing the electrically large target surface described in step (1) is as follows: the electrically large target surface is divided into H triangular facets, and each triangular facet and the triangular facets that share a common edge with it are combined into a pair of triangular facets to obtain W pairs of triangular facets, where H≥4 and is an even number.
3. The method according to claim 1, characterized in that, The RWG basis function f mentioned in step (1) w (r), whose expression is: Among them l w It represents two triangular facets and The length of the common side of the w-th triangular facet pair. and They represent and area, express The vector pointing from the vertex to the field point. express The vector pointing from the field point to the vertex.
4. The method according to claim 1, characterized in that, The σ mentioned in step (2) n The calculation formula is: Where π represents the value of a circle's circumference.
5. The method according to claim 1, characterized in that, The steps in step (3) are as follows: The RWG basis functions of the W triangle facet pairs are divided into S groups based on the octree algorithm. Based on the octree algorithm, objects containing electrically large targets are grouped into (x) min ,y min ,z min Starting from θ, the cube with side length θ is evenly divided into eight equal parts. Each of these eight parts is then evenly divided into eight equal parts again, and so on, until the side length of each resulting cube is less than θ. Then, select S cubes from all the cubes obtained from the last split that contain RWG basis functions. Next, group the RWG basis functions in each of these cubes into a group, resulting in the RWG basis function group Ω = {Ω1, Ω2, ..., Ω...}. s ,…,Ω S }, where the formula for calculating θ is: θ=max{|x max -x min |,|y max -y min |,|z max -z min |} Where x max x min y max y min z max , z min Let max{·} represent the maximum and minimum values of the x, y, and z axis components of the coordinates of the vertices of the H triangular facets, respectively. max{·} represents the maximum value operation, and |·| represents the modulus operation.
6. The method according to claim 1, characterized in that, The σ mentioned in step (3) n Each far-field block The calculation formula is: in, express The corresponding dimension is u y ×v y Approximate far-field block, and They represent The two corresponding low-rank matrices, ι y express effective rank, R y (σ n )express and The error matrix is denoted by ||·||, where ||·|| represents the Frobenus norm of the matrix, and ε represents the given error iteration threshold.
7. The method according to claim 1, characterized in that, The σ mentioned in step (3) n The impedance matrix Z(σ) at the point n The calculation formula is:
8. The method according to claim 1, characterized in that, Step (4) describes the calculation of the induced current vector I(k) on the electrically large target surface. b The w-th element I w (k b The implementation steps are as follows: (4a) Based on MoM, and through f w (r) and E in Calculate the excitation vector Γ(k0) at k0 and its κ-th derivative Γ. (k) The w-th element Γ of (k0) w (k0), Γ w (k) (k0): C w (k0)=∫f w (r)·E in doctor C w (κ) (k0)=∫E in(κ) ·f w (r)dr Among them, E in(κ) For E in The κ-th derivative, κ≥1; (4b) Calculate the element D in the α-th row and β-th column of the N×N dimensional pseudospectral differential matrix D based on the pseudospectral derivative method PSDM. αβ : in, For the αth and βth GLC sampling points σ α σ β The corresponding sampling points after coordinate transformation The coefficients are related to the values of α and β, where 1 ≤ α, β ≤ N; (4c) Calculate the Taylor expansion coefficient vectors of order 0 and q using the results of steps (4a) and (4b), and the Taylor expansion coefficient vectors of order m0 and m0. q : m0=Z -1 (k0)·Γ(k0) Where Z(k0) is the impedance matrix at the center wavenumber k0, Z -1 (k0) is the inverse matrix of Z(k0), Z (i) (k0) is the i-th order higher-order impedance matrix at the center wavenumber k0, D i Let i be the power of matrix D. For matrix D i The The elements of the row, i≥1, 1≤q≤Q; (4d) via m0, m q The w-th element m w,0 m w,q Calculate each k b The induced current vector I(k) at the location b The w-th element I w (k b ): L+O=Q Among them, L, a w,l They are polynomials The number of times and corresponding to I w (k b The coefficients of the l-th term, O, d w,o They are polynomials The number of times and corresponding to I w (k b The coefficient of the o-th term of Q is L, O, where L and O are positive integers. If Q is even, L = O; if Q is odd, the absolute value of the difference between L and O is |LO| = 1.
9. The method according to claim 1, characterized in that, Each wavenumber k described in step (4) b Radar cross section RCS(k) at the location b The calculation formula is: Where R represents the distance between the electrically large target surface and the radar. This indicates that R approaches infinity, j represents the imaginary unit, and η represents the wave impedance in air. The unit vector representing the observation point position vector r', Λ represents the outer surface of the electrically large target, ∫∫ Λ This indicates that the surface integral operation is performed on Λ, e represents the natural constant, and |·| represents the modulo operation.
Citation Information
Patent Citations
CBFM and AWE-based periodic structure broadband RCS acquisition method
CN116484642A
Method for estimating electromagnetic scattering characteristics of inhomogeneous medium target based on PMCHWT integral equation
CN115859613A
Novel electromagnetic property calculation method for electrically large metal target with uncertain appearance based on Maehly approximation
CN117436273A