A method for inverse design of computational spectral metasurface based on iterative optimization algorithm

Abstract

The present application relates to the technical field of nanophotonics and computational optical imaging, and particularly relates to a computational spectral metasurface reverse design method based on an iterative optimization algorithm, comprising: S1, physical structure parameterization, defining the geometric parameters, arrangement period and material refractive index of the nanostructure in the unit, initializing the physical structure parameters and the corresponding electromagnetic response operator; S2, constructing a differentiable simulation model, using a rigorous coupled wave analysis algorithm based on an automatic differentiation framework to calculate the complex amplitude transmission spectrum data under the current parameters, and keeping the electromagnetic response gradient flow continuous in the calculation graph. In the present application, through end-to-end soft and hard collaborative optimization of the metasurface physical structure parameters and the neural network reconstruction algorithm weight, the deep coupling and collaborative evolution of the hardware coding characteristics and the software decoding logic are realized, the performance mismatch problem caused by the traditional step-by-step design is fundamentally solved, and the restoration accuracy and recognition sensitivity of the system to complex spectral characteristics are significantly improved.

Description

Technical Field

[0001] This invention relates to the field of interdisciplinary technology of nanophotonics and computational optical imaging, and in particular to a computational spectral metasurface inverse design method based on iterative optimization algorithm. Background Technology

[0002] With the deep integration of micro-nano fabrication technology and computational optics theory, metasurface-based computational spectral detection systems have become a cutting-edge research direction in micro-spectrometers. Metasurfaces consist of arrays of subwavelength-scale nanostructures, enabling multidimensional fine modulation of the phase, amplitude, and polarization dimensions of the light field within extremely thin spatial scales. In computational spectral imaging systems, the metasurface acts as a physical encoding layer, responsible for mapping high-dimensional incident spectral information into a specific spatial intensity distribution. This distribution is then captured by a photoelectric sensor and accurately reconstructed using a backend digital unwinding algorithm. This architecture overcomes the size limitations of traditional spectroscopic elements, laying the physical foundation for chip-based, high-throughput spectral sensing applications.

[0003] However, existing technologies generally suffer from a disconnect between hardware response design and software reconstruction algorithms. Because hardware design and algorithm development are often treated as independent modules during the R&D phase, the physical coding characteristics of the metasurface and the backend decoding logic cannot achieve a deep match. As a result, the overall system performance is difficult to reach its theoretical limit, and both reconstruction accuracy and system robustness are insufficient when dealing with low signal-to-noise ratio or complex spectral environments. Summary of the Invention

[0004] To overcome the above shortcomings, this invention provides a computational spectral metasurface reverse design method based on iterative optimization algorithms, aiming to improve the performance mismatch caused by the lack of coordination between hardware design and reconstruction algorithms.

[0005] This invention provides the following technical solution: a method for reverse design of computational spectral metasurfaces based on iterative optimization algorithms, comprising:

[0006] S1. Physical structure parameterization: Define the geometric parameters, arrangement period and material refractive index of the nanostructure within the unit, and initialize the physical structure parameters and their corresponding electromagnetic response operators.

[0007] S2. Construct a differentiable simulation model, use a rigorous coupled-wave analysis algorithm based on an automatic differentiation framework to calculate the complex amplitude transmission spectrum data under the current parameters, and keep the electromagnetic response gradient flow continuous in the calculation graph.

[0008] S3. Establish a photoelectric conversion model, map the transmission spectrum data into a system transmission matrix, and combine the preset spectral data source and noise model to simulate and generate the coded spectral response signal received by the detector.

[0009] S4. Design a neural network reconstruction model, input the encoded spectral response signal into the neural network reconstruction model, and obtain the reconstructed spectral data through nonlinear feature mapping;

[0010] S5. Construct a joint loss function by weighting the mean square error between the reconstructed spectrum and the original spectrum, as well as the processability morphology constraint term for the minimum physical size limit of the nanostructure, to quantify the overall bias of the system.

[0011] S6. Perform soft and hard collaborative iterative optimization, use the backpropagation algorithm to synchronously calculate the gradient of the joint loss function with respect to physical parameters and network weights, and collaboratively update the physical structure and reconstruction algorithm until the preset convergence criterion is met.

[0012] Preferably, in step S1, the step of initializing the physical structure parameters and their corresponding electromagnetic response operators includes:

[0013] The unit is divided into The pixel grid region is defined, and an initial material density attribute value is assigned to each pixel, wherein the material density attribute value is limited to a continuous physical range of 0 to 1 to characterize the probability of the material's existence.

[0014] The pixel grid region is convolved using a dual smoothing filter operator. By controlling the filter radius and projection slope, the step discontinuities at the pixel edges are eliminated, generating a spatially continuous and gradient-differentiable material density distribution map.

[0015] Based on the effective medium theory, a refractive index mapping function is established to transform the material density distribution map into the corresponding complex permittivity matrix. The matrix is ​​then projected onto the reciprocal space using a two-dimensional fast Fourier transform to complete the initialization of the electromagnetic response operator.

[0016] Preferably, in step S2, the step of calculating the complex amplitude transmission spectrum data under the current parameters includes:

[0017] Construct the characteristic eigenvalue matrix of Maxwell's equations in the frequency domain, and use the electromagnetic response operator to fill the scattering coupling elements in the matrix to describe the momentum exchange of the light field between different diffraction orders;

[0018] The eigenvalue solving algorithm is used to perform full-space mode decomposition on the eigenmode matrix to obtain the propagation constants of each eigenmode in the metasurface layer and the normal mode vectors of the transverse electromagnetic field.

[0019] Based on the transfer matrix method, the continuity conditions of the electromagnetic field tangential components at the boundaries of each layer are matched, and the complex amplitude transmission coefficients at each wavelength point in the target band are obtained through interlayer iterative cascade calculation, thus forming the complex amplitude transmission spectrum data.

[0020] Preferably, in step S3, the step of simulating and generating the encoded spectral response signal received by the detector includes:

[0021] Perform a modulus square operation on the complex amplitude transmission spectrum data to extract the intensity transmittance of the metasurface at different wavelengths, and arrange them in wavelength sequence to construct the system transmission matrix;

[0022] A linear superposition operation is performed between the preset spectral data source and the system transmission matrix to calculate and generate a photocurrent response sequence under ideal conditions, which is used to characterize the expected value of photon count at the detector end.

[0023] A hybrid probability model incorporating Poisson statistical distribution and Gaussian white noise is introduced to randomly perturb the photocurrent response sequence under the ideal state, generating an encoded spectral response signal that conforms to the physical noise characteristics of the hardware.

[0024] Preferably, in step S4, the step of designing the neural network reconstruction model includes:

[0025] Construct a deep convolutional architecture with residual connections, and use an initial mapping layer to project the encoded spectral response signal into a high-dimensional latent space;

[0026] By using multiple cascaded residual blocks to perform nonlinear evolution on latent features, and by using a skip connection mechanism to preserve the original components of the signal during the depth mapping process, global morphological features of the spectral signal are extracted.

[0027] By standardizing the neuron outputs using a layer normalization operator, the propagation efficiency of the network gradient is optimized and the ability to suppress probe noise is enhanced.

[0028] Preferably, in step S5, the step of performing weighted mapping and quantifying the overall system deviation includes:

[0029] Calculate the root mean square error between the reconstructed spectral data and the original spectral data, and generate a spectral fidelity loss term to characterize the deconvolution accuracy of the system;

[0030] A total variational regularization operation is performed on the material density distribution map of the nanostructure, and a smoothing constraint term is constructed to suppress structural fragmentation by limiting the integral magnitude of the spatial gradient.

[0031] The material density distribution is projected using the Heaviside projection function. By controlling the slope of the spatial transition zone, the minimum size features in the structure are captured and constrained. The processability morphological constraint term is constructed, and it is linearly combined with the fidelity loss term and the smoothness constraint term to obtain the system comprehensive deviation.

[0032] Preferably, in step S6, the step of collaboratively updating the physical structure and reconstruction algorithm includes:

[0033] Initiate a chain-like differentiation process based on an automatic differentiation framework to simultaneously obtain the model gradient of the system's overall deviation relative to the internal weights of the neural network, as well as the structural gradient relative to the material density attribute value.

[0034] Using an adaptive optimizer with first-order and second-order momentum prediction capabilities, backpropagation updates are performed on the model gradient and the structural gradient respectively according to a dynamic learning rate strategy.

[0035] In each iteration, the geometric distribution of the hardware and the nonlinear mapping operator of the software are adjusted synchronously, so that the spectral coding characteristics of the hardware and the decoding characteristics of the software approach the global optimum in the joint search space.

[0036] Preferably, in step S6, the step of continuing until a preset convergence criterion is met includes:

[0037] The system monitors the change in the overall deviation of the system within a preset number of consecutive iterations in real time. When the rate of change of deviation is lower than the preset minimum threshold, an iteration termination command is triggered.

[0038] The current material density distribution map is subjected to binarization hard thresholding segmentation, which forces the continuous probability distribution to be transformed into discrete geometric boundaries representing solid materials or air gaps.

[0039] Perform the final system performance verification simulation to verify whether the binarized metasurface structure still meets the preset spectral reconstruction accuracy requirements after the superposition of manufacturing error disturbances, and output the final physical structure design model after the requirements are met.

[0040] The present invention has the following beneficial effects:

[0041] 1. In this invention, by performing end-to-end hardware and software co-optimization of the metasurface physical structure parameters and the neural network reconstruction algorithm weights, the deep coupling and co-evolution of hardware encoding characteristics and software decoding logic are realized, fundamentally solving the performance mismatch problem caused by traditional step-by-step design, and significantly improving the system's accuracy in restoring complex spectral features and its recognition sensitivity.

[0042] 2. In this invention, a hybrid noise operator containing Poisson statistics and Gaussian distribution is introduced into the photoelectric conversion model, and residual connection and layer normalization techniques are used in the reconstruction architecture. This enables the system to learn spontaneously and suppress the influence of detector dark current and photon shot noise during the iteration process, ensuring that the reconstructed data still has extremely high robustness and fidelity in high background noise environments.

[0043] 3. In this invention, by embedding a total variational regularization term and Heaviside projection constraints in the joint loss function, the spatial gradient of the nanostructure is effectively limited and the fuzzy intermediate state of the physical boundary is eliminated. This ensures that the geometry designed in reverse has good spatial coherence and meets the minimum size requirements of the actual processing technology, which greatly reduces the manufacturing difficulty of the device and improves the mass production yield. Attached Figure Description

[0044] Figure 1 This is a flowchart of a computational spectral metasurface reverse design method based on an iterative optimization algorithm proposed in this invention. Detailed Implementation

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

[0046] In embodiments of the present invention, the present invention provides a computational spectral metasurface reverse design method based on an iterative optimization algorithm, such as... Figure 1 As shown, it includes the following steps:

[0047] S1. Physical structure parameterization: Define the geometric parameters, arrangement period and material refractive index of the nanostructure within the unit, and initialize the physical structure parameters and their corresponding electromagnetic response operators.

[0048] Further, in step S1, the step of initializing the physical structure parameters and their corresponding electromagnetic response operators includes: dividing the unit into... The pixel grid region is defined, and an initial material density attribute value is assigned to each pixel. The material density attribute value is limited to a continuous physical range of 0 to 1 to characterize the probability of material existence. The pixel grid region is convolved using a dual smoothing filter operator. By controlling the filter radius and projection slope, the step discontinuity at the pixel edge is eliminated, generating a spatially continuous and gradient-differentiable material density distribution map. Based on the effective medium theory, a refractive index mapping function is established to transform the material density distribution map into the corresponding complex permittivity matrix. The matrix is ​​then projected onto the reciprocal space through a two-dimensional fast Fourier transform to complete the initialization of the electromagnetic response operator.

[0049] Specifically, firstly, within the physical space of the metasurface unit, according to the preset design cycle... Establish a two-dimensional Cartesian coordinate system, and discretize the space into a coordinate system with dimensions of 1. The pixel grid area. For each pixel... Assign a continuously varying initial material density property value. This constitutes the initial matter density matrix. ,in This density value corresponds to the material filling state of the nanostructure at a specific coordinate position, where 1 represents the filling medium material and 0 represents the filling background medium.

[0050] To ensure continuous gradient propagation during the reverse design process and eliminate boundary abrupt changes caused by discretization, the initial material density matrix... Perform convolution filtering. Use a filter with a specific radius. Conical kernel function or Gaussian kernel function Convolution with the matter density matrix yields the filtered density matrix. The calculation formula is as follows:

[0051] ;

[0052] In the above formula, Represented by pixels Centered on, with radius The neighborhood range, Represents pixels and neighboring points The Euclidean distance between them This is a preset weight allocation function.

[0053] Subsequently, a projection operator is used to control the binarization tendency of the filtered density, generating a spatially continuous and gradient-differentiable density distribution map. The calculation process is as follows:

[0054] ;

[0055] In the above formula, This is the projection slope coefficient, used to control the steepness of the density distribution evolution towards 0 or 1. This is the preset projection threshold.

[0056] Obtain the density distribution map of the material Subsequently, based on the effective medium theory, the geometric distribution of each point in space is mapped to physical optical properties. A linear mapping relationship between density values ​​and complex permittivity is established, yielding the permittivity matrix of the metasurface unit in real space. Its mapping formula is expressed as:

[0057] ;

[0058] In the above formula, The complex permittivity of the metasurface background medium is given by [insert value here]. The complex permittivity of the semiconductor or metallic materials used in nanostructures.

[0059] Finally, in order to initialize the electromagnetic response operator required for subsequent rigorous coupled-wave analysis, the dielectric constant matrix in real space is obtained by two-dimensional fast Fourier transform. Projecting onto the reciprocal lattice points of momentum space, the scattering coefficients of electromagnetic waves at different diffraction orders are obtained. The transformation formula is:

[0060] ;

[0061] In the above formula, Let be the area of ​​the metasurface unit. and They are respectively in direction and Reciprocal vector coordinates in the direction, This is the generated frequency-domain electromagnetic response operator, which serves as the input parameter for subsequent calculations of transmission spectrum data at various wavelengths.

[0062] This step, through parameterization and smoothing filtering, transforms the discrete geometric structure into a mathematically differentiable continuous field model, ensuring the lossless transfer of the optimization gradient between the physical structure and the electromagnetic response, and improving the convergence stability and physical distribution rationality of the reverse design algorithm.

[0063] S2. Construct a differentiable simulation model, use a rigorous coupled-wave analysis algorithm based on an automatic differentiation framework to calculate the complex amplitude transmission spectrum data under the current parameters, and keep the electromagnetic response gradient flow continuous in the calculation graph.

[0064] Further, in step S2, the step of calculating the complex amplitude transmission spectrum data under the current parameters includes: constructing the characteristic eigenvalue matrix of Maxwell's equations in the frequency domain, filling the scattering coupling elements in the matrix with electromagnetic response operators to describe the momentum exchange between different diffraction orders of the light field; performing full-space mode decomposition on the characteristic eigenvalue matrix using an eigenvalue solving algorithm to obtain the propagation constants of each eigenmode within the metasurface layer and the normal mode vector of the transverse electromagnetic field; matching the continuity conditions of the tangential components of the electromagnetic field at the boundaries of each layer based on the transfer matrix method, and obtaining the complex amplitude transmission coefficients at each wavelength point in the target band through interlayer iterative cascade calculations to constitute the complex amplitude transmission spectrum data.

[0065] Specifically, in the differentiable simulation stage, the frequency domain electromagnetic response operator generated in the previous physical structure parameterization stage is first retrieved. Construct the Teplitz matrix based on these discrete Fourier dilation coefficients. It is used to characterize the spatial distribution of dielectric constant within a metasurface unit.

[0066] Combined with the preset incident wavelength Incident angle and metasurface period Define a diagonal matrix and The diagonal elements of these matrices correspond to the normalized wavenumber components in the transverse direction for different diffraction orders. Subsequently, these matrix components are used to construct the eigenvalue matrix of Maxwell's equations in the frequency domain. To describe the momentum exchange and mode coupling of the light field at different diffraction orders, its construction formula is as follows:

[0067] ;

[0068] In the above formula, It is the identity matrix. It is the inverse of the dielectric constant matrix. This is the system matrix that couples the electric and magnetic field components.

[0069] After constructing the eigenvalue matrix, the differentiable eigenvalue solver is called within the automatic differentiation framework to solve the matrix. Perform full-space mode decomposition. This operator must ensure accurate computation of the gradients of eigenvalues ​​and eigenvectors with respect to the elements of the input matrix during backpropagation. Through decomposition computation, obtain a diagonal matrix composed of eigenvalues. and the matrix composed of eigenvectors Among them, matrix Each element in With the metasurface layer The propagation constant of each eigenmode The following relationship must be satisfied:

[0070] ;

[0071] The above propagation constant With eigenvector matrix Together, they constitute the physical description of the various normal modes within the metasurface.

[0072] Next, the continuity condition of the tangential component of the electromagnetic field at the boundary of each layer is matched based on the transfer matrix method or the scattering matrix method. For a thickness of... The metasurface layer defines the internal mode evolution matrix. Its diagonal elements are represented as ,in The free-space wavenumber is given. The coefficients of each diffraction order are solved by simultaneously applying the boundary condition equations for the incident region, the metasurface region, and the substrate region. For computational spectroscopy, the complex solution of the zeroth-order diffraction is extracted to obtain the wavelengths within the target band. Complex amplitude transmission coefficient at the location .

[0073] The final generated complex amplitude transmission spectrum data The definition is as follows:

[0074] ;

[0075] By discretizing the wavelengths within the target spectral range and repeating the above process, a complete complex amplitude transmission spectrum dataset is constructed.

[0076] This step, by embedding the numerical simulation process into the automatic differential calculation graph, achieves a precise positive mapping from physical parameters to electromagnetic response and ensures the continuity of gradient flow in the electromagnetic field evolution logic under complex geometric structures, providing physically real derivative information for subsequent hardware-software co-optimization.

[0077] S3. Establish a photoelectric conversion model, map the transmission spectrum data into a system transmission matrix, and combine the preset spectral data source and noise model to simulate and generate the coded spectral response signal received by the detector.

[0078] Further, in step S3, the step of simulating the generation of the encoded spectral response signal received by the detector includes: performing a modulus-square operation on the complex amplitude transmission spectrum data, extracting the intensity transmittance of the metasurface at different wavelengths, and arranging them in a wavelength sequence to construct a system transmission matrix; performing a linear superposition operation on the preset spectral data source and the system transmission matrix to calculate and generate a photocurrent response sequence under ideal conditions, which is used to characterize the expected value of photon count at the detector end; introducing a hybrid probability model containing Poisson statistical distribution and Gaussian white noise, randomly perturbing and sampling the photocurrent response sequence under ideal conditions, and generating an encoded spectral response signal that conforms to the hardware physical noise characteristics.

[0079] Specifically, in the photoelectric conversion modeling stage, the complex amplitude transmission spectrum data calculated in the previous step S2 is retrieved first. ,in The incident wavelength, This is a set of physical parameters characterizing the geometry of a metasurface. Modulus-square operations are performed on this complex data to calculate the intensity transmittance of the metasurface unit at a specific wavelength. Its calculation formula is expressed as:

[0080] ;

[0081] In the above formula, It is the conjugate complex number of the complex amplitude.

[0082] Subsequently, by performing traversal calculations on metasurface units under different structural parameters, the obtained multiple sets of intensity transmittance were arranged according to wavelength sequence. Spatial arrangement to construct the system transmission matrix In this matrix, the first... The line represents the first Spectral response curves of a type of metasurface structural unit, the first Columns represent the first Transmission efficiency at discrete wavelengths, matrix elements That is, the first Each detection channel corresponds to wavelength The intensity modulation response coefficient.

[0083] After establishing the system transfer matrix, obtain the preset raw spectral data source. This vector characterizes the incident light power distribution of the target light source at each discrete wavelength. The system transmission matrix... With the original spectral vector By performing linear superposition operations, the physical process of photons passing through the metasurface modulation layer and being received by detector pixels is simulated, and the photocurrent response sequence under ideal conditions is calculated and generated. The calculation formula is:

[0084] ;

[0085] In the above formula, the output vector Include Each component The corresponding detector array The total photon energy received by each pixel. This response sequence represents the expected photon count at the detector end, and is the original measurement benchmark unaffected by environmental and circuit noise contamination.

[0086] To simulate the detector's operating characteristics in a real physical environment, a hybrid probability model incorporating Poisson statistical distribution and Gaussian white noise is introduced. First, quantum shot noise during the detection process is considered. The expected photon count is randomly perturbed and sampled using a Poisson distribution operator, and its probability distribution formula is expressed as:

[0087] ;

[0088] In the above formula, The first one after adding Poisson noise Each channel signal value, The preset photoelectric conversion gain coefficient for the system.

[0089] Subsequently, a Gaussian white noise term representing the circuit's dark current and readout noise is superimposed on the Poisson noise. Generates the final encoded spectral response signal that conforms to the physical characteristics of the hardware. The expression is as follows:

[0090] ;

[0091] in, The photoelectric response signal after Poisson noise is added, and the noise term is... Follows the pattern of mean 0 and variance 0 The generated signal follows a normal distribution. As input data for subsequent reconstruction layers, it fully characterizes the actual output features of the metasurface detection system in complex background environments.

[0092] This step enables cross-scale mapping from electromagnetic field simulation data to digital electrical signals. By introducing a physically realistic hybrid noise sampling mechanism, it significantly enhances the robustness of the design scheme and the reliability of spectral decoding in real hardware environments.

[0093] S4. Design a neural network reconstruction model, input the encoded spectral response signal into the neural network reconstruction model, and obtain the reconstructed spectral data through nonlinear feature mapping;

[0094] Further, in step S4, the step of designing the neural network reconstruction model includes: constructing a deep convolutional architecture containing residual connections; projecting the encoded spectral response signal to a high-dimensional latent space using an initial mapping layer; performing nonlinear evolution on the latent features using multiple cascaded residual blocks; retaining the original components of the signal during the deep mapping process through a skip connection mechanism; and extracting the global morphological features of the spectral signal; and standardizing the neuron outputs through a layer normalization operator to optimize the propagation efficiency of the network gradient and enhance the ability to suppress detection noise.

[0095] Specifically, in the construction phase of the neural network reconstruction model, the noisy encoded spectral response signal from step S3 is first received. The signal has a length equal to the number of pixel channels of the detector. The signal is a one-dimensional vector. An initial mapping layer projects this low-dimensional measurement signal into a high-dimensional latent space to establish a preliminary characterization of the signal features. This mapping process is implemented through a fully connected linear operator combined with a nonlinear activation function, and its mathematical expression is:

[0096] ;

[0097] In the above formula, For size The projection weight matrix, As a preset latent spatial feature dimension, For bias vectors, For linear rectified functions, the output is This is the initial high-dimensional feature vector that will be used in subsequent calculations.

[0098] After initial spatial mapping, multiple cascaded residual blocks are used to nonlinearly evolve the latent features to extract complex global morphological features from the spectral signal. Each residual block contains two layers of one-dimensional convolution operators and corresponding nonlinear mapping units, and the block input is directly accumulated to the convolution output through a skip connection mechanism. For the first... For each residual block, the feature evolution logic follows the following formula:

[0099] ;

[0100] In the above formula, For the first Input features of each residual block, This represents a residual mapping operator composed of convolution and activation functions. This is the set of learnable parameters for this block. This skip connection mechanism ensures that the original linear components of the photoelectric signal are preserved during transmission in the deep network, effectively solving the gradient vanishing problem in deep networks, thereby accurately extracting the peak position and envelope contour of the spectral curve.

[0101] To optimize the propagation efficiency of network gradients and enhance the system's ability to suppress detector noise, a layer normalization operator is introduced after the convolution operation to standardize the neuron output. This operator, by calculating the mean and variance of the current feature layer, forces the feature distribution to be constrained within a stable numerical range. Its calculation formula is expressed as:

[0102] ;

[0103] In the above formula, and respectively, feature vectors Mean and variance along the channel dimension To prevent extremely small constants with a denominator of zero, and These are learnable scaling and translation parameters. This operation can effectively smooth numerical fluctuations caused by noise and improve the reconstruction stability of the model under different signal-to-noise ratio environments.

[0104] Finally, after The depth features processed by each residual block are transformed back into the spectral wavelength domain through an output mapping layer to generate the final reconstructed spectral data. The mapping relationship is defined as follows:

[0105] ;

[0106] In the above formula, To output the bias vector of the mapping layer, The output layer weight matrix has dimensions derived from the latent space dimensions. Number of target spectral sampling points The output is determined jointly. This refers to the spectral power distribution data corresponding to the wavelength coordinate system defined in step S1.

[0107] This step achieves high-fidelity decoding of the encoded signal through a residual learning architecture, which significantly improves the ability to restore details of complex multi-peak spectral structures while effectively suppressing detection noise.

[0108] S5. Construct a joint loss function by weighting the mean square error between the reconstructed spectrum and the original spectrum, as well as the processability morphology constraint term for the minimum physical size limit of the nanostructure, to quantify the overall bias of the system.

[0109] Further, in step S5, the step of performing weighted mapping and quantifying the overall system bias includes: calculating the root mean square error between the reconstructed spectral data and the original spectral data, and generating a spectral fidelity loss term to characterize the deconvolution accuracy of the system; performing total variational regularization on the material density distribution map of the nanostructure, and constructing a smoothing constraint term to suppress structural fragmentation by limiting the integral magnitude of the spatial gradient; using the Heaviside projection function to perform a projection transformation on the material density distribution, capturing and limiting the minimal size features in the structure by controlling the slope of the spatial transition zone, constructing a manufacturability morphology constraint term, and linearly combining it with the fidelity loss term and the smoothing constraint term to obtain the overall system bias.

[0110] Specifically, in the implementation phase of constructing the joint loss function, the reconstructed spectral data output from the preceding step S4 is first retrieved. and the preset raw spectral data source The root mean square error between the two is calculated to quantitatively evaluate the system's deconvolution accuracy of the encoded signal. Spectral fidelity loss term. The calculation formula is defined as follows:

[0111] ;

[0112] In the above formula, This represents the total number of spectral sampling points. To reconstruct the spectrum in the first Power value at each wavelength This corresponds to the original reference spectral value. The magnitude of this value directly reflects the degree of matching between the metasurface hardware coding features and the backend neural network decoding algorithm.

[0113] Physical constraints are applied to the topological morphology of the metasurface nanostructure to obtain the continuous mass density distribution map defined in step S1. A total variational regularization operation is performed on the distribution to suppress potential random fragmentation or isolated tiny fragments in the structure by limiting the integral magnitude of the spatial gradient. Smoothing constraint term. The construction formula is as follows:

[0114] ;

[0115] In the above formula, The computational domain for the metasurface element. The preset smoothing factor, and and represent the spatial rate of change of the density of the nanostructure material in the horizontal and vertical directions, respectively. By minimizing these terms, an optimization algorithm can be induced to generate geometric figures with clear physical boundaries and spatial continuity.

[0116] To further ensure that the generated nanostructures meet the minimum size requirements of photolithography or etching processes, the density distribution of the material is projected using the Heaviside projection function. This is achieved by adjusting the slope parameter in the projection operator. Force the density value to evolve towards either 0 or 1, and construct a processability morphological constraint term. This term is used to quantify the width of the transition zone at the structural boundary. It is achieved through the dispersion of the statistical density distribution, and its expression is defined as:

[0117] ;

[0118] In the above formula, This represents the pixel density value at the spatial coordinates. This value reaches its maximum when the pixel value is in a debinarized state near 0.5. By applying this constraint during iteration, blurred areas at the edges of nanostructures can be effectively suppressed, thereby capturing and limiting extremely small sharp corners or slits that are difficult to process within the structure.

[0119] Finally, the spectral fidelity loss term, smoothness constraint term, and manufacturability morphological constraint term are linearly weighted and summed to obtain the system comprehensive bias used to guide the gradient descent algorithm in performing global optimization. Its linear combination form is expressed as:

[0120] ;

[0121] In the above formula, , as well as These are the weighting coefficients for each loss term, used to balance the relative importance of spectral reconstruction accuracy, structural smoothness, and physical manufacturing feasibility. The output scalar value... It will be used as input to the backpropagation algorithm to drive the system parameters to evolve toward the global optimum.

[0122] This step achieves a deep mathematical unification of performance indicators and manufacturing constraints through weighted mapping of multi-dimensional objectives, ensuring that the reverse-designed structure has high reconstruction accuracy while possessing good mechanical strength and process reliability.

[0123] S6. Perform soft and hard collaborative iterative optimization, use the backpropagation algorithm to synchronously calculate the gradient of the joint loss function with respect to physical parameters and network weights, and collaboratively update the physical structure and reconstruction algorithm until the preset convergence criterion is met.

[0124] Furthermore, in step S6, the steps of collaboratively updating the physical structure and reconstruction algorithm include: initiating a chain-like differentiation process based on an automatic differentiation framework to simultaneously obtain the model gradient of the system's overall deviation relative to the internal weights of the neural network, and the structural gradient relative to the material density attribute value; using an adaptive optimizer with first-order momentum and second-order momentum prediction capabilities to perform backpropagation updates on the model gradient and structural gradient respectively according to a dynamic learning rate strategy; and simultaneously adjusting the geometric distribution of the hardware and the nonlinear mapping operator of the software in each iteration so that the spectral encoding characteristics of the hardware and the decoding characteristics of the software approach the global optimum in the joint search space.

[0125] Furthermore, in step S6, the steps until the preset convergence criterion is met include: real-time monitoring of the change in the overall system deviation within a preset number of consecutive iterations; triggering an iteration termination command when the deviation change rate is lower than a preset minimum threshold; performing binarization hard thresholding on the current material density distribution map to forcibly transform the continuous probability distribution into discrete geometric boundaries representing solid materials or air gaps; performing the final system performance verification simulation to verify whether the binarized metasurface structure still meets the preset spectral reconstruction accuracy requirements after being superimposed with manufacturing error disturbances, and outputting the final physical structure design model after meeting the requirements.

[0126] Specifically, in the execution phase of the hardware-software co-optimization, the joint loss function constructed in step S5 is first retrieved. Based on the automatic differentiation framework, a complete computational graph from output error to input parameters is established. A chained differentiation process is initiated to simultaneously calculate the joint loss function relative to the set of weight parameters within the neural network. Model gradient and the material density attribute value defined in step S1. structural gradient Its gradient calculation process follows the following derivative mapping relationship:

[0127] ;

[0128] In the above formula, the weight gradient Evolution of nonlinear mapping capability used to guide spectral reconstruction algorithms, structural gradient It is used to indicate the evolution direction of the metasurface geometry, ensuring that the physical modulation characteristics approach the optimal coding target.

[0129] After obtaining the gradient vector, an adaptive moment estimator optimizer with first-order and second-order momentum prediction capabilities is used to update the parameters. For each iteration step... Calculate the first moment estimate based on the current gradient. With second-order moment estimation The physical structure parameters and network weights are updated according to the dynamic learning rate strategy. The core update operator is expressed as follows:

[0130] ;

[0131] ;

[0132] ;

[0133] In the above formula, Let be the gradient vector at the current time. Representative includes and The set of parameters to be optimized For dynamic learning rate, and The attenuation coefficient is estimated based on the preset moment. This is a smoothing factor. In each iteration, the process synchronously adjusts the hardware's nanostructure arrangement and the software's decoding weights, forcing the hardware's spectral coding characteristics and the software's reconstruction algorithm to converge collaboratively within the joint search space.

[0134] During the iterative loop execution, continuous real-time monitoring is performed. Change in joint loss function within each iteration step When the rate of change of deviation meets the condition When the iteration terminates, a termination instruction is triggered, where This is the preset minimum convergence threshold. At this point, the current mass density distribution map is extracted. The binary hard thresholding process is performed, forcibly transforming the probability distribution within continuous intervals into discrete geometric boundaries representing solid media materials or air gaps. Its binary logic definition is as follows:

[0135] ;

[0136] In the above formula, the output result is... This constitutes the final metasurface physical layout data.

[0137] After binarization, the final system performance verification simulation is performed. Perturbation operators representing manufacturing tolerances are artificially injected into the generated binarized structure. This simulates critical dimensional deviations that may occur in actual photolithography or etching processes. The structure containing the perturbation is then resubmitted into steps S2 and S3 for forward electromagnetic verification. This verifies whether the encoded signal output by the system can still recover the spectral data meeting the preset accuracy requirements through the neural network in step S4 after the manufacturing error perturbation is superimposed. If the verification passes, the final physical structure design model is output.

[0138] This step enables the synchronous evolution of hardware parameters and algorithm weights. Through binarization processing and robustness verification, it ensures the determinism of the design scheme at the physical implementation level and maintains stable spectral reconstruction performance under complex process error environments.

[0139] Finally, it should be noted that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A computational spectral metasurface reverse design method based on an iterative optimization algorithm, characterized in that, include: S1. Physical structure parameterization: Define the geometric parameters, arrangement period and material refractive index of the nanostructure within the unit, and initialize the physical structure parameters and their corresponding electromagnetic response operators. S2. Construct a differentiable simulation model, use a rigorous coupled-wave analysis algorithm based on an automatic differentiation framework to calculate the complex amplitude transmission spectrum data under the current parameters, and keep the electromagnetic response gradient flow continuous in the calculation graph. S3. Establish a photoelectric conversion model, map the transmission spectrum data into a system transmission matrix, and combine the preset spectral data source and noise model to simulate and generate the coded spectral response signal received by the detector. S4. Design a neural network reconstruction model, input the encoded spectral response signal into the neural network reconstruction model, and obtain the reconstructed spectral data through nonlinear feature mapping; S5. Construct a joint loss function by weighting the mean square error between the reconstructed spectrum and the original spectrum, as well as the processability morphology constraint term for the minimum physical size limit of the nanostructure, to quantify the overall bias of the system. S6. Perform soft and hard collaborative iterative optimization, use the backpropagation algorithm to synchronously calculate the gradient of the joint loss function with respect to physical parameters and network weights, and collaboratively update the physical structure and reconstruction algorithm until the preset convergence criterion is met.

2. The method for reverse design of computational spectral metasurfaces based on iterative optimization algorithm according to claim 1, characterized in that, Step S1, the step of initializing the physical structure parameters and their corresponding electromagnetic response operators, includes: The unit is divided into The pixel grid region is defined, and an initial material density attribute value is assigned to each pixel, wherein the material density attribute value is limited to a continuous physical range of 0 to 1 to characterize the probability of the material's existence. The pixel grid region is convolved using a dual smoothing filter operator. By controlling the filter radius and projection slope, the step discontinuities at the pixel edges are eliminated, generating a spatially continuous and gradient-differentiable material density distribution map. Based on the effective medium theory, a refractive index mapping function is established to transform the material density distribution map into the corresponding complex permittivity matrix. The matrix is ​​then projected onto the reciprocal space using a two-dimensional fast Fourier transform to complete the initialization of the electromagnetic response operator.

3. The method for reverse design of computational spectral metasurfaces based on iterative optimization algorithm according to claim 1, characterized in that, Step S2, the step of calculating the complex amplitude transmission spectrum data under the current parameters, includes: Construct the characteristic eigenvalue matrix of Maxwell's equations in the frequency domain, and use the electromagnetic response operator to fill the scattering coupling elements in the matrix to describe the momentum exchange of the light field between different diffraction orders; The eigenvalue solving algorithm is used to perform full-space mode decomposition on the eigenmode matrix to obtain the propagation constants of each eigenmode in the metasurface layer and the normal mode vectors of the transverse electromagnetic field. Based on the transfer matrix method, the continuity conditions of the electromagnetic field tangential components at the boundaries of each layer are matched, and the complex amplitude transmission coefficients at each wavelength point in the target band are obtained through interlayer iterative cascade calculation, thus forming the complex amplitude transmission spectrum data.

4. The method for reverse design of computational spectral metasurfaces based on iterative optimization algorithm according to claim 1, characterized in that, Step S3, the step of simulating and generating the encoded spectral response signal received by the detector, includes: Perform a modulus square operation on the complex amplitude transmission spectrum data to extract the intensity transmittance of the metasurface at different wavelengths, and arrange them in wavelength sequence to construct the system transmission matrix; A linear superposition operation is performed between the preset spectral data source and the system transmission matrix to calculate and generate a photocurrent response sequence under ideal conditions, which is used to characterize the expected value of photon count at the detector end. A hybrid probability model incorporating Poisson statistical distribution and Gaussian white noise is introduced to randomly perturb the photocurrent response sequence under the ideal state, generating an encoded spectral response signal that conforms to the physical noise characteristics of the hardware.

5. The method for reverse design of computational spectral metasurfaces based on iterative optimization algorithm according to claim 1, characterized in that, Step S4, the step of designing the neural network reconstruction model, includes: Construct a deep convolutional architecture with residual connections, and use an initial mapping layer to project the encoded spectral response signal into a high-dimensional latent space; By using multiple cascaded residual blocks to perform nonlinear evolution on latent features, and by using a skip connection mechanism to preserve the original components of the signal during the depth mapping process, global morphological features of the spectral signal are extracted. By standardizing the neuron outputs using a layer normalization operator, the propagation efficiency of the network gradient is optimized and the ability to suppress probe noise is enhanced.

6. The method for reverse design of computational spectral metasurfaces based on iterative optimization algorithm according to claim 1, characterized in that, Step S5, the step of performing weighted mapping and quantifying the overall system deviation, includes: Calculate the root mean square error between the reconstructed spectral data and the original spectral data, and generate a spectral fidelity loss term to characterize the deconvolution accuracy of the system; A total variational regularization operation is performed on the material density distribution map of the nanostructure, and a smoothing constraint term is constructed to suppress structural fragmentation by limiting the integral magnitude of the spatial gradient. The material density distribution is projected using the Heaviside projection function. By controlling the slope of the spatial transition zone, the minimum size features in the structure are captured and constrained. The processability morphological constraint term is constructed, and it is linearly combined with the fidelity loss term and the smoothness constraint term to obtain the system comprehensive deviation.

7. The method for reverse design of computational spectral metasurfaces based on iterative optimization algorithm according to claim 1, characterized in that, Step S6, the step of collaboratively updating the physical structure and reconstruction algorithm, includes: Initiate a chain-like differentiation process based on an automatic differentiation framework to simultaneously obtain the model gradient of the system's overall deviation relative to the internal weights of the neural network, as well as the structural gradient relative to the material density attribute value. Using an adaptive optimizer with first-order and second-order momentum prediction capabilities, backpropagation updates are performed on the model gradient and the structural gradient respectively according to a dynamic learning rate strategy. In each iteration, the geometric distribution of the hardware and the nonlinear mapping operator of the software are adjusted synchronously, so that the spectral coding characteristics of the hardware and the decoding characteristics of the software approach the global optimum in the joint search space.

8. The method for reverse design of computational spectral metasurfaces based on iterative optimization algorithm according to claim 1, characterized in that, In step S6, the step of continuing until the preset convergence criterion is met includes: The system monitors the change in the overall deviation of the system within a preset number of consecutive iterations in real time. When the rate of change of deviation is lower than the preset minimum threshold, an iteration termination command is triggered. The current material density distribution map is subjected to binarization hard thresholding segmentation, which forces the continuous probability distribution to be transformed into discrete geometric boundaries representing solid materials or air gaps. Perform the final system performance verification simulation to verify whether the binarized metasurface structure still meets the preset spectral reconstruction accuracy requirements after the superposition of manufacturing error disturbances, and output the final physical structure design model after the requirements are met.

Images