CN107145667B - Rapid calculation method of wave front structure function - Google Patents

Rapid calculation method of wave front structure function Download PDF

Info

Publication number
CN107145667B
CN107145667B CN201710316115.5A CN201710316115A CN107145667B CN 107145667 B CN107145667 B CN 107145667B CN 201710316115 A CN201710316115 A CN 201710316115A CN 107145667 B CN107145667 B CN 107145667B
Authority
CN
China
Prior art keywords
wavefront
surface shape
original surface
structure function
formula
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201710316115.5A
Other languages
Chinese (zh)
Other versions
CN107145667A (en
Inventor
安其昌
张景旭
杨飞
赵宏超
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Changchun Institute of Optics Fine Mechanics and Physics of CAS
Original Assignee
Changchun Institute of Optics Fine Mechanics and Physics of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Changchun Institute of Optics Fine Mechanics and Physics of CAS filed Critical Changchun Institute of Optics Fine Mechanics and Physics of CAS
Priority to CN201710316115.5A priority Critical patent/CN107145667B/en
Publication of CN107145667A publication Critical patent/CN107145667A/en
Application granted granted Critical
Publication of CN107145667B publication Critical patent/CN107145667B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Telescopes (AREA)
  • Complex Calculations (AREA)

Abstract

A quick calculation method for a wavefront structure function relates to the field of wavefront quality evaluation and solves the problem that the existing calculation method cannot effectively represent the full-frequency-domain characteristics of the surface shape of a primary mirror of a large-aperture telescope. The method comprises the following steps: calculating structural functions of two orthogonal directions of the original surface shape of a primary mirror of a large-aperture telescope; step two, rotating the original surface shape; step three, interpolating the rotated original surface shape into a regular grid, and repeating the step two until the original surface shape is rotated for multiple times; and step four, obtaining a final structure function by a square average method. The invention utilizes the mirror structure function to carry out frequency domain expression and combines the power spectrum to carry out calculation, thereby greatly improving the calculation efficiency.

Description

Rapid calculation method of wave front structure function
Technical Field
The invention relates to the technical field of wavefront quality evaluation, in particular to a rapid calculation method of a wavefront structure function.
Background
With the development of large-aperture telescopes, the wave front root mean square is a characteristic which cannot effectively represent the whole frequency range by only relying on the wave front root mean square. The power spectrum is short for power spectral density function, and is defined as signal power in a unit frequency band, and represents the variation of the signal power with frequency, i.e. the distribution of the signal power in the frequency domain. The power spectrum evaluation method was proposed by NIF at the end of the last century and has been a more mature standard (ISO 10110). However, the power spectrum evaluation method can only evaluate the fluctuation of a wavefront (a surface formed by a particle which starts to be displaced at a certain moment when a wave propagates in a medium, which is called a wavefront and represents a spatial position where wave energy arrives at a certain moment, and which is moving) in a certain direction, and for an ultra-smooth surface, the certain direction can be taken as a standard for evaluation. However, for large aperture telescopes, due to cost constraints and atmospheric effects, the surface profile does not necessarily reach sub-nanometer scale, and at this scale, the anisotropy of the surface profile becomes significant over larger scales. And then, a two-dimensional power spectrum is introduced for the research of the power spectrum, and finally, the surface shape is evaluated by collapsing to one dimension. However, the fourier transform must be processed on orthogonal data, and the statistical characteristics of the surface shapes at other angles cannot be evaluated. On the other hand, in the low frequency part, the power spectrum evaluation method lacks averaging to reduce noise, so the evaluation effect is also affected.
The structure function is proposed in 1965 to describe the atmospheric turbulence under different scales, fundamentally, the structure function represents the second-order center distance on different calculation scales, and all points of the wave front need to participate in calculation on different calculation scales.
Therefore, it is necessary to design a fast algorithm for the wavefront structure function.
Disclosure of Invention
The invention provides a rapid calculation method of a wavefront structure function, aiming at solving the problem that the existing calculation method cannot effectively represent the full-frequency-domain characteristics of the surface shape of a primary mirror of a large-aperture telescope.
The technical scheme adopted by the invention for solving the technical problem is as follows:
the invention discloses a method for quickly calculating a wave front structure function, which comprises the following steps:
calculating structural functions of two orthogonal directions of the original surface shape of a primary mirror of a large-aperture telescope;
step two, rotating the original surface shape;
step three, interpolating the rotated original surface shape into a regular grid, and repeating the step two until the original surface shape is rotated for multiple times;
and step four, obtaining a final structure function by a square average method.
Further, the specific process of the step one is as follows:
according to wave front structure function DwavefrontThe structural functions of the primary mirror of the large-aperture telescope in two orthogonal directions are defined and calculated; wave front structure function DwavefrontIs as defined in formula (1):
Figure BDA0001288493680000021
in formula (1): phi, which represents the wavefront obtained by detection with a wavefront sensor, is a two-dimensional matrix,
Figure BDA0001288493680000031
and
Figure BDA0001288493680000032
are all vectors of positions on the wave front,<·>represents the average over the wavefront.
Further, the specific process of the second step is as follows:
① determining the number M of the original surface shape rotation;
② let the space coordinate of the point on the original surface be (x)i,yi,zi) The spatial coordinate after the rotation by the angle θ becomes (x)b,yb,zb) As shown in formula 2:
Figure BDA0001288493680000033
③ calculating equation (2) by using the definition of the expected operation, and obtaining the convolution form as shown in equation (3):
Figure BDA0001288493680000034
in the formula (3), the window function
Figure BDA0001288493680000035
w denotes the conjugate of the window function w, Φ denotes the conjugate of the wavefront Φ, r is the position vector on the wavefront
Figure BDA0001288493680000036
The mold of (4);
④ fast Fourier transform
The fast fourier transform of equation (3) is performed by the convolved fast fourier transform being equal to the product of the fast fourier transform:
Figure BDA0001288493680000037
further, the specific process of determining the number M of rotations of the original surface shape is as follows:
one fast Fourier transform is equivalent to calculation for 4-neighborhood, one rotation is 8-neighborhood, and two rotations are equivalent to calculation of 12 points on the circumference, and so on; and taking the area enclosed by the average power spectrum obtained by rotating for different times and the horizontal axis as a characteristic quantity, and determining the rotation times M required by the original surface shape when the increment of the characteristic quantity is less than 10% after the rotation times are increased.
Further, the specific process of step three is as follows:
interpolating the rotated wavefront into a regular grid by interpolation to correspond to the grid before rotation, wherein the resulting interpolation error is represented by the difference in RMS values of the rotated wavefront, the rotated θ and re-interpolated wavefront becomes
Figure BDA0001288493680000041
Theta is the angle of rotation of the original profile,
Figure BDA0001288493680000042
the angle of the original surface shape before rotation is obtained, and meanwhile, the formula (5) is obtained according to the definition of the power spectrum PSD:
Figure BDA0001288493680000043
in the formula (5), the reaction mixture is,
Figure BDA0001288493680000044
is a spatial frequency vector.
Further, the specific process of step four is as follows:
① obtaining expression (6) by performing an inverse fast Fourier transform using expressions (3), (4) and (5):
Figure BDA0001288493680000045
in the formula (6), θiRepresenting the angle of the ith rotation of the original surface shape;
② the formula (6) is brought into the formula (1) to calculate the final wave front structure function, which is shown in the formula (7):
Figure BDA0001288493680000046
the invention has the beneficial effects that:
for the calculation of the wavefront structure function, all data within a certain radius can be evaluated, and in actual operation, due to the discreteness of the data, a specific embodiment needs to be considered carefully. The basic idea of the invention is as follows: firstly, the wavefront to be analyzed is rotated, the power spectrum of the wavefront is calculated, and finally, the obtained power spectrums are averaged.
1. Compared with the prior method for evaluating the surface shape of the system by using the power spectrum, the power spectrum essentially analyzes frequency information along a certain dimension, and researches information of two orthogonal dimensions as an improved two-dimensional power spectrum. The structural function describes information at a certain interval in any direction, and the representation is more comprehensive. The invention utilizes the mirror structure function to carry out frequency domain expression and combines the power spectrum to carry out calculation, thereby greatly improving the calculation efficiency.
2. The structural function can be better linked to the optical transfer function than the power spectrum. For the actual optical system of the large-aperture telescope, a structural function curve can be obtained according to the specific situation of each optical element, and finally, the final error curve can be synthesized through simple addition. The errors caused by the mechanical and tracking can be uniformly considered through the optical transfer function and error distribution is carried out.
3. By the invention, a system engineer can control the intermediate frequency error of the system, reasonably distribute limited precision indexes and predict the optical information transmission performance of the actual system; meanwhile, different functions of different frequency components in system errors can be more comprehensively understood, and limited resources can be more fully utilized.
Drawings
FIG. 1 is a schematic diagram illustrating a method for fast calculating a wavefront structure function according to the present invention.
Detailed Description
The present invention will be described in further detail with reference to the accompanying drawings.
The invention provides a method for quickly calculating a wave front structure function by utilizing Fourier transform and a power spectrum, which has the basic idea that the rotation of an original matrix is realized by utilizing a rotation matrix; after rotation, new data and original data are aligned by interpolation operation, and finally a final wave front structure function is obtained by a square average method.
As shown in fig. 1, the method for fast calculating a wavefront structure function of the present invention includes the following specific steps:
1. first according to the wave front structure function DwavefrontThe structural functions of the two orthogonal directions of the original surface shape of the primary mirror of the large-aperture telescope are calculated. Wave front structure function DwavefrontAs shown in formula (1):
Figure BDA0001288493680000051
in formula (1): phi, which represents the wavefront obtained by detection with a wavefront sensor, is a two-dimensional matrix,
Figure BDA0001288493680000052
and
Figure BDA0001288493680000061
are all vectors of positions on the wave front,<·>mean over wavefront, content in expression, the same as below, for example: in the formula (1), the reaction mixture is,
Figure BDA0001288493680000062
to represent
Figure BDA0001288493680000063
Averaged over the wavefront.
2. Determining the number M of required rotations of the original surface shape
As shown in fig. 1, one fft corresponds to the calculation for the 4-neighborhood, one rotation for the 8-neighborhood, two rotations for the 12 points on the circle, and so on. The area enclosed by the average power spectrum obtained by rotating for different times and the horizontal axis is taken as a characteristic quantity, when the number of times of rotation is increased, the increment of the characteristic quantity is less than 10 percent, the number of times of rotation M required by the original surface shape can be determined, and the proper average number of times of rotation can be directly selected according to experience.
3. Rotation of original profile
Let the spatial coordinates of the points on the original surface shape be (x)i,yi,zi) The space coordinate after rotating a certain angle theta becomes (x)b,yb,zb) That is, the relative rotation of the θ angle between the two is performed, the equation (2) is calculated by using the basic definition of the expected operation, and the obtained convolution form is very similar to the convolution, and the obtained convolution form is shown in equation (3):
Figure BDA0001288493680000064
Figure BDA0001288493680000065
in the formula (3), the window function: (
Figure BDA0001288493680000066
w represents: the conjugate of the window function w, Φ denotes: the conjugate of the wave front phi,
Figure BDA0001288493680000067
is a position vector on the wavefront, r is a position vector on the wavefront
Figure BDA0001288493680000068
The die of (1).
4. Fast Fourier Transform (FFT)
By the convolved fast fourier transform being equal to the product of the fast fourier transform, the fast fourier transform of equation (3) can be obtained:
Figure BDA0001288493680000069
5. after the original surface shape is rotated, interpolation operation is needed to be carried out on the data, the rotated wave front is interpolated into a regular grid by utilizing the interpolation operation, the regular grid is enabled to correspond to the grid before the rotation, the generated interpolation error is represented by the difference of RMS (root mean square) values of the wave front before and after the rotation, the wave front which is rotated by theta and re-interpolated is changed into the wave front which is rotated by theta
Figure BDA0001288493680000071
Theta is the angle of rotation required for the original profile,
Figure BDA0001288493680000072
the angle of the original surface shape before rotation is obtained, and meanwhile, the angle can be obtained according to the basic definition of a Power Spectrum (PSD):
Figure BDA0001288493680000073
in the formula (5), the reaction mixture is,
Figure BDA0001288493680000074
is a space frequency vector (circle/rad).
6. And repeating the steps 3 to 5 until the original surface shape is rotated for M times.
7. And (3) obtaining a final wave front structure function by a square average method, wherein the formula is shown as a formula (7). Specifically, the method comprises the following steps: first, Inverse Fast Fourier Transform (IFFT) is performed using equations (3), (4), and (5) to obtain:
Figure BDA0001288493680000075
in the formula (6), θiRepresents: angle of ith rotation of original surface shape.
Then, the formula (6) is taken into the formula (1) to calculate to obtain a final wave front structure function, which is shown in the formula (7):
Figure BDA0001288493680000076
it should be noted that: in the above method for fast calculation of the wave front structure function, the meaning of the same letter representation in all formulas is the same.
For the evaluation means of the power spectrum, an expression of the two-dimensional power spectrum can be obtained according to the basic definition of the power spectrum. By comparing the structural function with the calculation process of the power spectrum, the main differences can be found in that: 1. firstly, calculating the error with smaller scale by the structural function, and starting the power spectrum by the representation of the low-frequency error; 2. the structure function needs to calculate the error of each direction of a certain scale, and the power spectrum only calculates the error of two orthogonal directions. The comparison shows that the calculation precision and efficiency of the wave front structure function are greatly improved.
The foregoing is only a preferred embodiment of the present invention, and it should be noted that, for those skilled in the art, various modifications and decorations can be made without departing from the principle of the present invention, and these modifications and decorations should also be regarded as the protection scope of the present invention.

Claims (5)

1. A method for rapidly calculating a wave front structure function is characterized by comprising the following steps:
calculating structural functions of two orthogonal directions of the original surface shape of a primary mirror of a large-aperture telescope;
according to wave front structure function DwavefrontThe structural functions of the primary mirror of the large-aperture telescope in two orthogonal directions are defined and calculated; wave front structure function DwavefrontIs as defined in formula (1):
Figure FDA0002370137400000011
in formula (1): phi, which represents the wavefront obtained by detection with a wavefront sensor, is a two-dimensional matrix,
Figure FDA0002370137400000012
and
Figure FDA0002370137400000013
are all vectors of positions on the wave front,<·>represents the average over the wavefront;
step two, rotating the original surface shape;
step three, interpolating the rotated original surface shape into a regular grid, and repeating the step two until the original surface shape is rotated for multiple times;
and step four, obtaining a final structure function by a square average method.
2. The method for rapidly calculating the wavefront structure function according to claim 1, wherein the specific process of the second step is as follows:
① determining the number M of the original surface shape rotation;
② let the space coordinate of the point on the original surface be (x)i,yi,zi) Rotation of thetaThe spatial coordinate after the angle becomes (x)b,yb,zb) As shown in formula 2:
Figure FDA0002370137400000014
③ calculating equation (2) by using the definition of the expected operation, and obtaining the convolution form as shown in equation (3):
Figure FDA0002370137400000015
in the formula (3), the window function
Figure FDA0002370137400000016
w denotes the conjugate of the window function w, Φ denotes the conjugate of the wavefront Φ, r is the position vector on the wavefront
Figure FDA0002370137400000017
The mold of (4);
④ fast Fourier transform
The fast fourier transform of equation (3) is performed by the convolved fast fourier transform being equal to the product of the fast fourier transform:
Figure FDA0002370137400000021
3. the method for rapidly calculating the wavefront structure function according to claim 2, wherein the specific process for determining the number M of rotations of the original surface shape is as follows:
one fast Fourier transform is equivalent to calculation for 4-neighborhood, one rotation is 8-neighborhood, and two rotations are equivalent to calculation of 12 points on the circumference, and so on; and taking the area enclosed by the average power spectrum obtained by rotating for different times and the horizontal axis as a characteristic quantity, and determining the rotation times M required by the original surface shape when the increment of the characteristic quantity is less than 10% after the rotation times are increased.
4. The method for rapidly calculating the wavefront structure function according to claim 2, wherein the specific process of the third step is as follows:
interpolating the rotated wavefront into a regular grid by interpolation to correspond to the grid before rotation, wherein the resulting interpolation error is represented by the difference in RMS values of the rotated wavefront, the rotated θ and re-interpolated wavefront becomes
Figure FDA0002370137400000022
Theta is the angle of rotation of the original profile,
Figure FDA0002370137400000023
the angle of the original surface shape before rotation is obtained, and meanwhile, the formula (5) is obtained according to the definition of the power spectrum PSD:
Figure FDA0002370137400000024
in the formula (5), the reaction mixture is,
Figure FDA0002370137400000025
is a spatial frequency vector.
5. The method for rapidly calculating the wavefront structure function according to claim 4, wherein the specific process of the fourth step is as follows:
① obtaining expression (6) by performing an inverse fast Fourier transform using expressions (3), (4) and (5):
Figure FDA0002370137400000026
in the formula (6), θiRepresenting the angle of the ith rotation of the original surface shape;
② the formula (6) is brought into the formula (1) to calculate the final wave front structure function, which is shown in the formula (7):
Figure FDA0002370137400000031
CN201710316115.5A 2017-05-08 2017-05-08 Rapid calculation method of wave front structure function Active CN107145667B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710316115.5A CN107145667B (en) 2017-05-08 2017-05-08 Rapid calculation method of wave front structure function

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710316115.5A CN107145667B (en) 2017-05-08 2017-05-08 Rapid calculation method of wave front structure function

Publications (2)

Publication Number Publication Date
CN107145667A CN107145667A (en) 2017-09-08
CN107145667B true CN107145667B (en) 2020-07-03

Family

ID=59776916

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710316115.5A Active CN107145667B (en) 2017-05-08 2017-05-08 Rapid calculation method of wave front structure function

Country Status (1)

Country Link
CN (1) CN107145667B (en)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114034470B (en) * 2021-11-10 2022-09-20 中国科学院长春光学精密机械与物理研究所 Telescope wavefront rotation angle calculation method and device and telescope

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101017178A (en) * 2006-11-15 2007-08-15 中国科学院安徽光学精密机械研究所 Image drifting velocity method for measuring average wind speed and direction of atmosphere
CN104614083A (en) * 2014-12-20 2015-05-13 佛山市南海区欧谱曼迪科技有限责任公司 Method for recovering phase distribution of phase shift interference figures and method for obtaining phase shift between two figures
CN105425378A (en) * 2015-12-31 2016-03-23 中国科学院光电技术研究所 Virtual-aperture complex-amplitude splicing super resolution astronomical telescope system
CN106529104A (en) * 2016-12-28 2017-03-22 哈尔滨工业大学 Phase screen simulation method of short distance transmission of light in underwater turbulent flow

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8103045B2 (en) * 2005-01-04 2012-01-24 Stc.Unm Structure function monitor

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101017178A (en) * 2006-11-15 2007-08-15 中国科学院安徽光学精密机械研究所 Image drifting velocity method for measuring average wind speed and direction of atmosphere
CN104614083A (en) * 2014-12-20 2015-05-13 佛山市南海区欧谱曼迪科技有限责任公司 Method for recovering phase distribution of phase shift interference figures and method for obtaining phase shift between two figures
CN105425378A (en) * 2015-12-31 2016-03-23 中国科学院光电技术研究所 Virtual-aperture complex-amplitude splicing super resolution astronomical telescope system
CN106529104A (en) * 2016-12-28 2017-03-22 哈尔滨工业大学 Phase screen simulation method of short distance transmission of light in underwater turbulent flow

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
"基于结构函数的子孔径拼接算法研究";安其昌等;《红外与激光工程》;20150331;第44卷(第3期);第929-932页 *

Also Published As

Publication number Publication date
CN107145667A (en) 2017-09-08

Similar Documents

Publication Publication Date Title
CN107315160B (en) Relatively prime array Wave arrival direction estimating method based on interpolation virtual array signal atom norm minimum
CN105181121B (en) Using the high-precision near field acoustic holography method of weighted iteration equivalent source method
CN110045323B (en) Matrix filling-based co-prime matrix robust adaptive beamforming algorithm
CN106443587B (en) A kind of high-resolution quick deconvolution sound source imaging algorithm
CN107450047B (en) Compressed sensing DOA estimation method based on unknown mutual coupling information under nested array
CN107561484B (en) Direction-of-arrival estimation method based on interpolation co-prime array covariance matrix reconstruction
CN110197112B (en) Beam domain Root-MUSIC method based on covariance correction
KR101524793B1 (en) Apparatus and method for estimating direction of arrival using array antenna
CN108008355B (en) Steady sound source positioning method based on quaternary orthogonal microphone array
CN111220222B (en) Measuring algorithm for flow of ultrasonic gas meter
CN104391183A (en) Near-field-measurement-based rapid calculation method of antenna far-field region characteristic
CN109085556B (en) High-frequency ground wave radar wave field forming method based on first-order and second-order peak ratios
CN109343003B (en) Method for identifying sound source formed by fast iterative shrinking wave beams
CN104535987A (en) Amplitude phase error self-correcting method applicable to uniform circular array acoustic susceptance system
CN109709510A (en) A kind of estimation method and system of coherent 2-d direction finding
EP2201564A1 (en) Enhanced sound source localization system and method by using a movable microphone array
CN107145667B (en) Rapid calculation method of wave front structure function
CN109870669A (en) How soon a kind of two dimension claps mesh free compression Wave beam forming identification of sound source method
CN112285647A (en) Signal orientation high-resolution estimation method based on sparse representation and reconstruction
CN107894581A (en) A kind of wideband array Wave arrival direction estimating method
CN108614234B (en) Direction-of-arrival estimation method based on multi-sampling snapshot co-prime array received signal fast Fourier inverse transformation
CN104105049A (en) Room impulse response function measuring method allowing using quantity of microphones to be reduced
CN112087235A (en) Sparsity self-adaptive DOA estimation method and system based on pseudo-inverse perception dictionary
CN109613474B (en) Angle measurement compensation method suitable for short-distance vehicle-mounted radar
CN116299150B (en) Two-dimensional DOA estimation method of dimension-reduction propagation operator in uniform area array

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant