SPR grating type sensor performance optimization method based on IMPA

Abstract

The invention relates to an SPR grating type sensor performance optimization method based on IMPA, in particular to a method for optimizing the structural performance of a Kretschmann type SPR grating sensor, and relates to the field of sensors and operational research. The SPR grating type sensor performance optimization method is characterized by comprising the following steps that (1) the sensor structure and parameters of the sensor structure are determined; (2) determining constraint conditions in sensor structure design; (3) determining an objective function of sensor structure optimization design; (4) solving and initializing model parameters; and (5) carrying out IMPA iterative solution on the optimization model. According to the method, the problem of multi-parameter joint optimization in the prior art is solved, the optimization theory in operational research and the design process of an actual sensor are deeply fused, multiple optimization targets such as sensitivity, response speed and stability of the sensor are comprehensively considered, the MPA algorithm is applied to sensor optimization model solution, and the optimization efficiency of the sensor is improved. The advantages of MPA in high-dimensional and nonlinear optimization problems are fully played, and a new method is provided for structural performance optimization of the Kretschmann type SPR grating sensor.

Description

Technical Field

[0001] The invention relates to the field of sensors and operations research, and mainly relates to a method for optimizing the structural performance of a Kretschmann type SPR grating sensor. Background Art

[0002] As an important implementation form of surface plasmon resonance sensing technology, Kretschmann-type SPR (Surface Plasmon Resonance, SPR) grating sensor has shown high sensitivity and real-time monitoring capabilities in the field of optical sensing. However, its detection performance is still limited by many factors. First, in terms of sensitivity, although the current technology can detect changes at the nanoscale, when dealing with low-concentration substance detection or weak molecular interactions, the signal response is weak and the problem of insufficient sensitivity is still difficult to overcome. Secondly, the performance of the sensor is closely related to the structural design, including multiple key parameters such as biofilm thickness, metal film thickness, grating tooth number, concave-convex spacing, incident angle, and incident light wavelength. The current mainstream design methods are mostly based on empirical adjustment or single variable optimization, lacking systematic multi-parameter joint optimization methods, resulting in the failure to fully exert the performance. In addition, the sensor is easily disturbed by external factors such as temperature changes, mechanical vibrations, and unstable light sources in complex environments, which reduces the signal accuracy and stability. Furthermore, although the currently widely used metal materials (such as gold and silver) have good plasma properties, they are limited in specific band ranges, chemical stability, and corrosion resistance, which limits their performance in some special applications. Therefore, although Kretschmann-type SPR sensors have important value, there is still room for improvement in structural performance optimization and sensitivity enhancement.

[0003] Kretschmann-type SPR grating sensors have shown extensive potential in practical applications, especially in the fields of biomedicine, environmental monitoring and food safety. In the biomedical field, the sensor is widely used to detect biomolecular interactions, such as antigen-antibody binding, DNA hybridization and drug screening. Such applications require the sensor to accurately capture weak signals of low-concentration biomolecules while having high selectivity and resolution. In environmental monitoring, SPR grating sensors are used to detect pollutants in the environment (such as heavy metal ions in water and toxic gases in the air). Such applications require the sensor to provide stable and anti-interference signal output in complex and changing environments. In the field of food safety, SPR sensors are used to quickly detect harmful substances in food (such as pesticide residues or pathogenic microorganisms), which puts higher requirements on detection speed and batch detection capabilities. In general, practical applications have put forward diverse requirements on the sensitivity, stability, response speed and environmental adaptability of Kretschmann-type SPR grating sensors, becoming a new direction for future technological improvements.

[0004] In the face of the current bottleneck of Kretschmann-type SPR grating sensors in performance optimization, operations research and intelligent optimization algorithms provide strong support. Operations research methods are good at dealing with complex optimization problems with multiple objectives and multiple constraints. They can jointly optimize multiple key parameters such as grating period, metal film thickness and incident angle to help find the design with the best performance. At the same time, intelligent optimization algorithms (such as particle swarm optimization, genetic algorithm and simulated annealing algorithm) have significant advantages in solving high-dimensional nonlinear problems. They can effectively avoid the local optimal trap in traditional optimization methods through global search and rapid iteration, and further improve the efficiency and accuracy of parameter optimization. In addition, intelligent optimization algorithms can combine the simulation model and experimental data of the sensor, and continuously improve the design through data-driven methods to ensure that the optimization scheme can adapt to the complex needs of practical applications. By selecting metal materials and regulating the grating structure, the environmental stability and working bandwidth of the sensor can be improved while maintaining sensitivity, thereby meeting the multiple requirements of modern optical sensing technology for high sensitivity, high stability and fast response. Combining operations research and intelligent optimization methods with the performance optimization of Kretschmann-type SPR grating sensors can not only break through the performance bottleneck of the current Kretschmann-type SPR grating sensors, but also promote the innovative application of this technology in more fields. Summary of the invention

[0005] In view of the problems existing in the above-mentioned prior art, the technical problem to be solved by the present invention is to provide a method for optimizing the performance of an SPR grating type sensor based on IMPA, the specific process of which is as follows: Figure 1 shown.

[0006] The technical solution implementation steps are as follows:

[0007] (1) Determine the sensor structure and its parameters:

[0008] The traditional Kretschmann type SPR grating sensor structure is usually composed of four parts: "substance to be measured-metal film-substrate-optical path system". However, with the in-depth study of SPR grating sensors, it is found that the surface of ordinary single metal film layer SPR grating sensors is relatively active and easily reacts with the solution to be measured during the detection process, thereby reducing the sensitivity of the sensor. In response to this problem, this paper attempts to improve the performance of the SPR sensor by adding other film materials as auxiliary film layers on the surface of the metal film in the structural design of the SPR grating sensor. The studied SPR grating sensor structure consists of five parts: "substance to be measured-biofilm-metal film-substrate-optical path system".

[0009] Based on the five-layer structure of "substance to be tested - biofilm - metal film - substrate layer - optical path system", the specific parameters we need to design are as follows:

[0010] Biofilm thickness d b ;

[0011] Metal film thickness d m ;

[0012] Grating teeth number N g ;

[0013] The convex-concave spacing of the grating is Λ;

[0014] ·Incident wavelength λ;

[0015] Angle of incidence θ;

[0016] (2) Determine the constraints in the sensor structure design:

[0017] 1) Physical restrictions

[0018] Biofilm thickness d b :The thickness of the biofilm has a certain reasonable range. If it is too thick, it may affect the molecular diffusion rate and sensor response time. If it is too thin, it may not be able to effectively fix the molecules to be tested or provide sufficient physical and chemical modification. Assume that the lower limit of the biofilm thickness is d bmin , the upper limit is d bmax ,but:

[0019] d bmin ≤d b ≤d bmax

[0020] Metal film thickness d m:The thickness of the metal film also has its limit. If it is too thin, it may not be able to effectively excite the surface plasmon resonance. If it is too thick, it will increase the cost and may cause signal attenuation. Let the lower limit of the metal film thickness be d mmin , the upper limit is d mmax ,but:

[0021] d mmin ≤d m ≤d mmax

[0022] Grating teeth number N g : The number of grating teeth must be a positive integer and is limited by factors such as the grating size. Assume that the lower limit of the number of grating teeth is N gmin , the upper limit is N gmax , and is an integer, then:

[0023] N gmin ≤N g ≤N gmax ,

[0024] The convex-concave spacing of the grating Λ: The convex-concave spacing must be within a reasonable range to ensure the optical performance of the grating, and is related to the number of grating teeth N. g There is an inverse relationship, and the lower limit of the concave-convex spacing is Λ min , the upper limit is Λ max ,but:

[0025]

[0026] Λ min ≤Λ≤Λ max

[0027] Incident wavelength λ: The incident light wavelength must be within the wavelength range that the sensor can detect, assuming the lower limit is λ min , with an upper limit of λ max ,but:

[0028] λ min ≤λ≤λ max

[0029] ·Incident angle θ: The incident angle usually has to be within a certain angle range to effectively stimulate the SPR effect, that is, to reach a specific exit angle θ at a specific wavelength λ m , can be adjusted by the grating equation:

[0030] (sin(θ)+sin(θ m ))Λ=mλ

[0031] Where m is the diffraction order, θ m is the diffraction angle, and Λ is the concave-convex spacing.

[0032] 2) Environmental restrictions

[0033] Effect of temperature T: Assume that the sensor is within the temperature variation range (T min ≤T≤T max ) when working, the signal change rate is required Cannot exceed a certain threshold ΔS Tmax ,but:

[0034]

[0035] Where S is the sensitivity, ΔT is the change in temperature, and ΔS is the change in sensitivity.

[0036] Effect of humidity H: Assume that the sensor is within the temperature variation range (H min ≤H≤H max ) when working, the signal change rate is required Cannot exceed a certain threshold ΔS Hmax ,but:

[0037]

[0038] Where S is the sensitivity, ΔH is the change in humidity, and Δs is the change in sensitivity.

[0039] (3) Determine the objective function of sensor structure optimization design:

[0040] 1) Maximize sensitivity s

[0041] Assume that the concentration of the substance to be measured is C, ΔS represents the change in signal response, ΔC represents the change in concentration, and define In the sensor setting, we often seek to maximize the sensitivity, that is, max(S). However, in order to avoid the problem of negative sensitivity caused by algorithm convergence, we perform a logarithmic transformation on the sensitivity, that is, max(ln(S+1)). At the same time, considering that there may be certain measurement errors in the signal response change and concentration change in actual measurement, the weight ω is introduced. s1 and ω s2 To adjust the importance of the two, the objective function can be further optimized as follows:

[0042] max(ω s1 ×ln(ΔS+1)-ω s2 ×ln(ΔC+1))

[0043] 2) Stability optimization

[0044] Stability can be measured by the degree of signal fluctuation when the sensor is working under different environmental conditions or for a long time. Let σ be the standard deviation of signal fluctuation. In order to make the stability index have a good balance with the sensitivity and other indicators in the objective function, the standard deviation is normalized, that is, (σ max is the maximum standard deviation determined based on prior knowledge or experimental data), and the objective function is min(σ n ). At the same time, in order to avoid the indicator having too large a weight in the objective function due to a small standard deviation in the early stage of optimization, an exponential form can be used, namely:

[0045]

[0046] 3) Optimize response speed

[0047] Response speed r It is defined as the time required from the time when the substance to be tested contacts the sensor to the time when the signal reaches a stable value. In order to make the response speed index and other indicators have a unified trend in the objective function (the smaller the better), the reciprocal of the response time can be taken. The objective function is set as:

[0048]

[0049] 4) Establishment of comprehensive objective function

[0050] Taking the above factors into consideration, the objective function of multi-objective optimization is constructed as follows:

[0051]

[0052] Where ω1, ω2, and ω3 are the weights corresponding to sensitivity, stability, and response speed, respectively.

[0053] (4) Model parameter solution initialization:

[0054] After constructing the above complete "decision variables-constraints-objective function" operations research model, we consider using the Marine Predators Algorithm (MPA) to solve the optimal sensor parameters of the model. MPA is a new meta-heuristic optimization algorithm proposed by Afshin Faramarzi et al. in 2020. The algorithm is inspired by the foraging strategies of marine predators, especially Lévy and Brownian motion, and the optimal encounter rate strategy between predators and prey. MPA simulates the theory of survival of the fittest in the ocean and has strong evolutionary ability, fast search speed and strong optimization ability. The optimization process of MPA includes three main stages to simulate the movement behavior of predators and prey throughout their life cycle. Each stage considers a different speed ratio, which enables MPA to more accurately simulate the movement relationship between predators and prey and provide better performance in the optimization process.

[0055] MPA first uses a set of random solutions as the initial search space, and generates the initial solution according to the formula based on the range of the search space limit. The optimal solution is selected as the top predator, and a matrix called Elite is constructed. This matrix array monitors the search process and searches for prey based on the location information of the prey. The second matrix is ​​the Prey matrix, which has the same dimension as the Elite, and the predator is updated based on this matrix. At the end of each iteration, if a predator with a higher fitness value appears, the current top predator will be replaced and the elite matrix will be updated accordingly.

[0056] Similar to most metaheuristics, MPA is a population-based method where the initial solutions are randomly distributed over the search space. In the first step:

[0057] X0=X min +rand(X max -X min )

[0058] Where X max and X min are the lower and upper bounds of the variable, and rand is a uniform random vector in the range of 0 to 1.

[0059] However, the standard MPA random initialization may cause the initial position of the population to be relatively scattered, so the individuals in the population may need to go through multiple iterations to gradually gather and optimize, thereby increasing the number of iterations and time complexity of the algorithm. Therefore, we improve it and use Logistics-Tent chaos initialization to ensure that the initial solution is generated more uniformly. The specific improvements are as follows:

[0060]

[0061] In the formula, r∈(0,4), P0∈(0,1), the chaotic sequence {P n}, and then use this chaotic sequence to initialize the position of the population as follows:

[0062] X0=X min +{P n}(X max -X min )

[0063] During the initialization process, the best solution is selected as the top predator and the elite matrix is ​​constructed. During the search process, the solutions after the survival of the fittest are named top predators, and they form a matrix, also called the elite matrix. The array of this matrix supervises the search and finds the prey based on the information of the prey's location.

[0064]

[0065] Where X is the top predator vector, which is replicated n times to construct the elite matrix. n is the number of search agents and d is the dimension. It is worth noting that both predators and prey are considered search agents. At the end of each iteration, if the top predator is replaced by a better predator, the elite will be updated.

[0066] Another matrix with the same dimension as the "elite matrix" is called the prey matrix, and the predator updates its position based on this matrix. In short, initialization creates the initial prey, and the predator constructs the elite. The prey is shown in the following formula:

[0067]

[0068] X i,j represents the jth dimension of the i-th prey, n is the number of search agents, and d is the dimension. The optimization process is mainly directly related to these two matrices.

[0069] (5) IMPA iterative solution optimization model:

[0070] Before introducing the steps of the proposed IMPA (Improved Marine Predators Algorithm, IMPA) algorithm, it is necessary to know the mathematical models of the two random motions of Lévy and Brownian motion. Then the optimal foraging, movement and memory management strategies of marine predators are introduced one by one and the corresponding mathematical models are established, and the overall process of the IMPA algorithm is gradually given.

[0071] 1) Brownian motion

[0072] Standard Brownian motion is a random process in which their step sizes are drawn from a probability function defined by a normal (Gaussian) distribution with zero mean (μ = 0) and unit variance (σ 2 =1). The probability density function (PDF) of the motion at point x is as follows:

[0073]

[0074] 2) Lévy Flight

[0075] A Lévy flight is a random walk whose step size is determined by the probability function defined by the Lévy distribution (power-law tail):

[0076] L(x j )≈|x j | 1-α

[0077] where x j is the step size, 1<α≤2 is the power law exponent, and the probability density function of the Levy stable process in integral form is as follows:

[0078]

[0079] Where α defines the distribution index and the scaling properties of the control process, and γ is the selected scaling unit. The integral in the formula has an analytical solution only in a few cases. When α = 2, it represents a Gaussian distribution, and when α = 1, it represents a Cauchy distribution. The solution of the integral in the formula generally requires the use of a series expansion method only when the x value is large, as shown in the following formula:

[0080]

[0081] Wherein Γ is the gamma function, α is an integer, and Γ(1+α)=α.

[0082] Mantegna proposed an exact and fast algorithm for generating Levy stable processes under exponential distribution (α) between 0.3 and 1.99. The Mantegna method is used to generate random numbers based on Levy distribution as follows:

[0083]

[0084] In Lévy flight, the area of ​​interest is traced primarily by small steps associated with long hops, whereas Brownian motion, by contrast, can cover the area of ​​interest with more uniform and controlled steps.

[0085] 3) Early stage of iteration

[0086] At high speed ratios or when the predator moves faster than the prey, the predator does not move at all. This happens in the initial iterations of the optimization, which is the key to exploration. The mathematical model for this stage is:

[0087]

[0088]

[0089]

[0090] in is a vector containing random numbers from a normal distribution, representing Brownian motion, represents moment multiplication, P=0.5 is a constant, R is a standard random number, Iter represents the current number of iterations, Max_ Iter represents the maximum number of iterations, Represents the step length of each movement of prey i.

[0091] 4) Mid-Iteration

[0092] At unit speed ratio or when both predator and prey are moving at the same speed. At this point both the predator and prey are simulating looking for prey. At this stage, both exploration and exploitation are important. Therefore, half of the population is designated for exploration and the other half is designated for exploitation. At this stage, the prey is responsible for exploitation and the predator is responsible for exploration. According to this rule, when the prey is doing Lévy motion, the best strategy for the predator is Brownian motion.

[0093] For the first half of the population:

[0094]

[0095]

[0096]

[0097] in is a random number vector based on the Levy distribution, representing the Levy motion. The product of the prey and the prey simulates the movement of the prey in the Levy way, Iter represents the current number of iterations, Max_ Iter represents the maximum number of iterations, Represents the step length of each movement of prey i.

[0098] For the second half of the population:

[0099]

[0100]

[0101] in is an adaptive parameter that controls the predator's movement step size.

[0102] 5) Late iteration

[0103] At low speed ratios or when the predator moves faster than the prey, the optimal strategy for the predator is Lévy motion. This occurs at the end of the optimization process, which is mainly associated with high exploitation capabilities. At low speed ratios, this stage can be described as:

[0104]

[0105]

[0106]

[0107] and The multiplication of simulates the Lévy motion of the predator, while The step length is increased in the position to simulate the movement of the predator. Iter represents the current number of iterations, Max_ Iterrepresents the maximum number of iterations, Represents the step length of each movement of prey i.

[0108] 6) Vortex formation and the influence of FADs

[0109] Another point that can cause changes in the behavior of marine predators is environmental issues, such as vortex formation or the fish aggregating device (FADs) effect. FADs can be thought of as local optima that capture these points in the search space. Considering these longer jumps during the simulation can avoid stagnation in local optima. Therefore, the FADs effect can be expressed mathematically as:

[0110]

[0111] Where FADs = 0.2 is the probability that FADs affects the optimization process. U is a binary vector, and the array includes 0 and 1. This generates a random vector in [0,1], and if the array is less than 0.2, it is changed to 0, and if the array is greater than 0.2, it is changed to 1. r is a uniform random number in [0,1]. X max and X min are the dimension vectors of the lower and upper bounds, respectively. The subscripts r1 and r2 represent the random indices of the prey matrix, respectively.

[0112] 7) Memories in the Ocean

[0113] Marine predators have a good memory to remind them of places where they have successfully foraged before, and this ability is simulated by storing it in the MPA. After updating the prey and realizing the effect of FADs, the fitness of this matrix is ​​evaluated to update the elite position. The fitness of each solution in the current iteration is compared to its solution in the previous iteration, and if the fitness of the current solution is greater, the current solution will replace the old solution. This process improves the quality of the solution over time and also simulates the predator returning to areas where prey is abundant and successfully foraging.

[0114] 8) Multi-subgroup improvement strategy

[0115] However, for the problems of low solution accuracy and easy to fall into local optimum in the standard MPA algorithm, a multi-subpopulation improvement strategy is proposed. After the population of marine predators is iterated, each population is sorted according to the size of the obtained population fitness value, and the population is divided into three subpopulations. The individuals with the fitness values ​​ranked in the top N / 3 are high-quality populations, those with fitness values ​​between N / 3 and 2N / 3 are classified as follower populations, and the rest are classified as low-quality populations.

[0116] In a high-quality population, excellent individuals influence each other in two different dimensions by learning from each other. This learning process helps excellent individuals to conduct reinforcement learning based on their own experience, making them more focused on areas with a higher probability of optimal solutions, thereby accelerating the convergence speed of the algorithm and improving the accuracy of optimization.

[0117]

[0118] is the new position of the individual in the high-quality population, P1 is a standard random number, and are any two dimensions in an individual.

[0119] For the follower population, reverse learning strategy update is used to improve the algorithm's ability to escape from local optimality.

[0120]

[0121] is the new position of the individual in the follower population. min and X max are the upper and lower limits of the problem dimension.

[0122]

[0123] is the new position of the individual in the low-quality population, and They are the global optimal solution and the solution with average fitness respectively. r1 and r2 are standard random numbers. The global optimal solution participates in the update of the low-quality population, which speeds up the convergence of the algorithm.

[0124] In summary, after multiple updates of the fitness function, the model will output the optimal sensor parameter design values.

[0125] The present invention has the following advantages over the prior art:

[0126] (1) Combination of operations research model and sensor optimization design: This invention innovatively integrates the optimization theory in operations research with the actual sensor design process. Traditional sensor design often focuses on the optimization of a single parameter, while this invention comprehensively considers multiple optimization objectives, such as key performance indicators such as sensor sensitivity, response speed, and stability, through the application of operations research models. This multi-objective optimization method can maximize the efficiency and cost-effectiveness of the design while ensuring that the sensor performance meets the requirements.

[0127] (2) Innovative application of IMPA in solving models: This invention innovatively applies the marine predator algorithm to the solution of complex sensor optimization models, giving full play to the advantages of IMPA in high-dimensional, nonlinear optimization problems. IMPA can better avoid the trouble of local optimal solutions in the solution process and achieve rapid convergence of global optimal solutions. IMPA simulates the predation behavior of marine predators, maintains a balance between global search and local refined exploration, and has strong global search capabilities and high computational efficiency. Therefore, the use of the IMPA algorithm to optimize sensor design not only improves the accuracy of the solution, but also speeds up the optimization process, greatly improving the design performance and efficiency of the sensor. BRIEF DESCRIPTION OF THE DRAWINGS

[0128] In order to better understand the present invention, further description is given below in conjunction with the accompanying drawings.

[0129] Figure 1 It is a flowchart of the steps for establishing a performance optimization method for SPR grating type sensors based on IMPA;

[0130] Figure 2 It is the algorithm flow chart of IMPA;

[0131] Figure 3 It is the result of optimizing the design of 4 sets of sensors by using the present invention; DETAILED DESCRIPTION

[0132] The present invention is further described in detail below through implementation cases.

[0133] This implementation case uses four different Kretschmann-type SPR grating sensor design structures. For different structures, the model is used for specific optimization. The overall process of the SPR grating sensor performance optimization method based on IMPA provided by the present invention is as follows: Figure 1 As shown, the specific steps are as follows:

[0134] (1) Determine the sensor structure and its parameters:

[0135] The traditional Kretschmann type SPR grating sensor structure is usually composed of four parts: "substance to be measured-metal film-substrate-optical path system". However, with the in-depth study of SPR grating sensors, it is found that the surface of ordinary single metal film layer SPR grating sensors is relatively active and easily reacts with the solution to be measured during the detection process, thereby reducing the sensitivity of the sensor. In response to this problem, this paper attempts to improve the performance of the SPR sensor by adding other film materials as auxiliary film layers on the surface of the metal film in the structural design of the SPR grating sensor. The studied SPR grating sensor structure consists of five parts: "substance to be measured-biofilm-metal film-substrate-optical path system".

[0136] Based on the five-layer structure of "substance to be tested - biofilm - metal film - substrate layer - optical path system", the specific parameters we need to design are as follows:

[0137] Biofilm thickness d b , assuming that in this experiment, d bmin =10nm,d bmax =50nm.

[0138] Metal film thickness d m , assuming that in this experiment, d mmin =20nm,d mmax =80nm.

[0139] Grating teeth number N g , assuming that in this experiment, N gmin =10, N gmax =50, and is an integer.

[0140] The convex-concave spacing of the grating is Λ. Assume that in this experiment, Λ min =0.02μm, Λ max = 0.1 μm, and

[0141] ·Incident wavelength λ, assuming that in this experiment, λ min =400nm,λ max =800nm.

[0142] The incident angle θ is determined by the grating diffraction equation.

[0143] (2) Determine the constraints in the sensor structure design:

[0144] 1) Physical restrictions

[0145] Biofilm thickness d b :

[0146] 10nm≤d b ≤50nm

[0147] Metal film thickness d m :

[0148] 20nm≤d m ≤80nm

[0149] Grating teeth number N g :

[0150] 10≤N g ≤50,

[0151] The convex-concave spacing of the grating:

[0152]

[0153] 0.02μm≤Λ≤0.1μm

[0154] Incident wavelength λ:

[0155] 400nm≤λ≤800nm

[0156] ·Incident angle θ: Achieve a specific exit angle θ at different wavelengths λ m , can be adjusted by the grating equation:

[0157] (sin(θ)+sin(θ m ))Λ=mλ

[0158] Where m is the diffraction order, θ m is the diffraction angle, and Λ is the concave-convex spacing.

[0159] 2) Environmental restrictions

[0160] Influence of temperature T: When the sensor works within the temperature range (-10°≤T≤50°), the signal change rate is required to be Cannot exceed a certain threshold ΔS Tmax =0.1, then:

[0161]

[0162] Where S is the sensitivity, ΔT is the change in temperature, and ΔS is the change in sensitivity.

[0163] Effect of humidity H: When the sensor works within the temperature variation range (20% ≤ H ≤ 80%), the signal change rate is required to be Cannot exceed a certain threshold ΔS Hmax =0.25, then:

[0164]

[0165] Where S is the sensitivity, ΔH is the change in humidity, and ΔS is the change in sensitivity.

[0166] (3) Determine the objective function of sensor structure optimization design:

[0167] 1) Maximize sensitivity S

[0168] Assume that the concentration of the substance to be measured is C, ΔS represents the change in signal response, ΔC represents the change in concentration, and define In the sensor setting, we often seek to maximize the sensitivity, that is, max(S). However, in order to avoid the problem of negative sensitivity caused by algorithm convergence, we perform a logarithmic transformation on the sensitivity, that is, max(ln(S+1)). At the same time, considering that there may be certain measurement errors in the signal response change and concentration change in actual measurement, the weight ω is introduced. s1 =0.6 and ω s2 = 0.4 to adjust the importance of the two, and the objective function can be further optimized as follows:

[0169] max(0.6×ln(ΔS+1)-0.4×ln(ΔC+1))

[0170] 2) Stability optimization

[0171] Stability can be measured by the degree of signal fluctuation when the sensor is working under different environmental conditions or for a long time. Let σ be the standard deviation of signal fluctuation. In order to make the stability index have a good balance with the sensitivity and other indicators in the objective function, the standard deviation is normalized, that is, (σ max is the maximum standard deviation determined based on prior knowledge or experimental data, σ max =0.5), the objective function is min(σ n ). At the same time, in order to avoid the indicator having too large a weight in the objective function due to a small standard deviation in the early stage of optimization, an exponential form can be used, namely:

[0172]

[0173] 3) Optimize response speed

[0174] Response speed r It is defined as the time required from the time when the substance to be tested contacts the sensor to the time when the signal reaches a stable value. In order to make the response speed index and other indicators have a unified trend in the objective function (the smaller the better), the reciprocal of the response time can be taken. The objective function is set as:

[0175]

[0176] 4) Establishment of comprehensive objective function

[0177] Taking the above factors into consideration, let ω1 = 0.5, ω2 = 0.2, ω3 = 0.3, and construct the objective function of multi-objective optimization as follows:

[0178]

[0179] Where ω1, ω2, and ω3 are the weights corresponding to sensitivity, stability, and response speed, respectively.

[0180] (4) Model parameter solution initialization:

[0181] The Logistics-Tent chaos initialization is used to ensure that the initial solution is generated more evenly. The formula is:

[0182]

[0183] In the formula, r = 3, P0 = 0.3, and the chaotic sequence {P n}, and then use this chaotic sequence to initialize the position of the population as follows:

[0184] X0=X min +{P n}(X max -X min )

[0185] Where, X min =0,X max =10;

[0186] During the initialization process, the best solution is selected as the top predator and the elite matrix is ​​constructed. During the search process, the solutions after the survival of the fittest are named top predators, and they form a matrix, also called the elite matrix. The array of this matrix supervises the search and finds the prey based on the information of the prey's location.

[0187]

[0188] Where X is the top predator vector, which is replicated n times to construct the elite matrix. n = 100, d = 6. It is worth noting that both predators and prey are considered search agents. At the end of each iteration, if the top predator is replaced by a better predator, the elite will be updated.

[0189] Another matrix with the same dimension as the "elite matrix" is called the prey matrix, and the predator updates its position based on this matrix. In short, initialization creates the initial prey, and the predator constructs the elite. The prey is shown in the following formula:

[0190]

[0191] X i,j represents the jth dimension of the i-th prey, n=100, d=6, and the optimization process is mainly directly related to these two matrices.

[0192] (5) IMPA iterative solution optimization model:

[0193] Before introducing the steps of the proposed IMPA algorithm, it is necessary to know the mathematical models of two random motions: Lévy and Brownian motion. Then, the optimal foraging, movement, and memory management strategies of marine predators are introduced one by one and the corresponding mathematical models are established, and the overall process of the IMPA algorithm is gradually given.

[0194] 1) Brownian motion

[0195] Standard Brownian motion is a random process in which their step sizes are drawn from a probability function defined by a normal (Gaussian) distribution with zero mean (μ = 0) and unit variance (σ 2 =1). The probability density function (PDF) of the motion at point x is as follows:

[0196]

[0197] 2) Lévy Flight

[0198] A Lévy flight is a random walk whose step size is determined by the probability function defined by the Lévy distribution (power-law tail):

[0199] L(x j )≈|x j | 1-α

[0200] where x j is the step size, 1<α≤2 is the power law exponent, and the probability density function of the Levy stable process in integral form is as follows:

[0201]

[0202] Where α defines the distribution index and the scaling properties of the control process, and γ is the selected scaling unit. The integral in the formula has an analytical solution only in a few cases. When α = 2, it represents a Gaussian distribution, and when α = 1, it represents a Cauchy distribution. The solution of the integral in the formula generally requires the use of a series expansion method only when the x value is large, as shown in the following formula:

[0203]

[0204] Wherein Γ is the gamma function, α is an integer, and Γ(1+α)=α.

[0205] Mantegna proposed an exact and fast algorithm for generating Levy stable processes under exponential distribution (α) between 0.3 and 1.99. The Mantegna method is used to generate random numbers based on Levy distribution as follows:

[0206]

[0207] In Lévy flight, the area of ​​interest is traced primarily by small steps associated with long hops, whereas Brownian motion, by contrast, can cover the area of ​​interest with more uniform and controlled steps.

[0208] 3) Early stage of iteration

[0209] Let Max_ Iter = 1000, at high speeds or when the predator moves faster than the prey, the predator does not move at all. This happens in the initial iteration of the optimization, which is the key to exploration. The mathematical model for this stage is:

[0210]

[0211]

[0212]

[0213] in is a vector containing random numbers from a normal distribution, representing Brownian motion, represents moment multiplication, P=0.5 is a constant, and R is a standard random number.

[0214] 4) Mid-Iteration

[0215] At unit speed ratio or when both predator and prey are moving at the same speed. At this point both the predator and prey are simulating looking for prey. At this stage, both exploration and exploitation are important. Therefore, half of the population is designated for exploration and the other half is designated for exploitation. At this stage, the prey is responsible for exploitation and the predator is responsible for exploration. According to this rule, when the prey is doing Lévy motion, the best strategy for the predator is Brownian motion.

[0216] For the first half of the population:

[0217]

[0218]

[0219]

[0220] in is a random number vector based on the Levy distribution, representing the Levy motion. The product with the prey simulates the movement of the prey in the manner of Levy.

[0221] For the second half of the population:

[0222]

[0223]

[0224] in is an adaptive parameter that controls the predator's movement step size.

[0225] 5) Late iteration

[0226] At low speed ratios or when the predator moves faster than the prey, the optimal strategy for the predator is Lévy motion. This occurs at the end of the optimization process, which is mainly associated with high exploitation capabilities. At low speed ratios, this stage can be described as:

[0227]

[0228]

[0229]

[0230] and The multiplication of simulates the Lévy motion of the predator, while The step length is increased in position to simulate the movement of the predator.

[0231] 6) Vortex formation and the influence of FADs

[0232] Another point that can cause changes in the behavior of marine predators is environmental issues, such as vortex formation or the fish aggregating device (FADs) effect. FADs can be thought of as local optima that capture these points in the search space. Considering these longer jumps during the simulation can avoid stagnation in local optima. Therefore, the FADs effect can be expressed mathematically as:

[0233]

[0234] Where FADs = 0.2 is the probability that FADs affects the optimization process. U is a binary vector, and the array includes 0 and 1. This generates a random vector in [0,1], and if the array is less than 0.2, it is changed to 0, and if the array is greater than 0.2, it is changed to 1. r is a uniform random number in [0,1]. X max and X min are the dimension vectors of the lower and upper bounds, respectively. The subscripts r1 and r2 represent the random indices of the prey matrix, respectively.

[0235] 7) Memories in the Ocean

[0236] Marine predators have a good memory to remind them of places where they have successfully foraged before, and this ability is simulated by storing it in IMPA. After updating the prey and implementing the FADs effect, a fitness evaluation is performed on this matrix to update the elite positions. The fitness of each solution in the current iteration is compared to its solution in the previous iteration, and if the fitness of the current solution is greater, the current solution will replace the old solution. This process improves the quality of the solution over time and also simulates the predator returning to areas where prey is abundant and successfully foraging.

[0237] 8) Multi-subgroup improvement strategy

[0238] However, for the problems of low solution accuracy and easy to fall into local optimum in the standard MPA algorithm, a multi-subpopulation improvement strategy is proposed. After the population of marine predators is iterated, each population is sorted according to the size of the obtained population fitness value, and the population is divided into three subpopulations. The individuals with the fitness values ​​ranked in the top N / 3 are high-quality populations, those with fitness values ​​between N / 3 and 2N / 3 are classified as follower populations, and the rest are classified as low-quality populations.

[0239] In a high-quality population, excellent individuals influence each other in two different dimensions by learning from each other. This learning process helps excellent individuals to conduct reinforcement learning based on their own experience, making them more focused on areas with a higher probability of optimal solutions, thereby accelerating the convergence speed of the algorithm and improving the accuracy of optimization.

[0240]

[0241] is the new position of the individual in the high-quality population, P1 is a standard random number, and are any two dimensions in an individual.

[0242] For the follower population, reverse learning strategy update is used to improve the algorithm's ability to escape from local optimality.

[0243]

[0244] is the new position of the individual in the follower population. min and X max are the upper and lower limits of the problem dimension.

[0245]

[0246] is the new position of the individual in the low-quality population, and They are the global optimal solution and the solution with average fitness respectively. r1 and r2 are standard random numbers. The global optimal solution participates in the update of the low-quality population, which speeds up the convergence of the algorithm.

[0247] In summary, after multiple updates of the fitness function, the model will output the optimal sensor parameter design values.

[0248] In order to verify the effectiveness of the present invention in optimizing the structure of the Kretschmann type SPR grating sensor, four sets of sensor optimization design experiments were carried out. The experimental results are shown in the figure. Figure 3 As shown. Figure 3 It can be seen that the Kretschmann grating type SPR sensor design method established in the present invention improves the sensor performance under a specific structure by more than 95%. Specifically, for the optimization design of six parameters, namely, biofilm thickness, metal film thickness, grating tooth number, concave-convex spacing, incident angle, and incident light wavelength, we use the quantitative index Goos- Displacement distance means, Figure 3 The result data further indicate that our design method can achieve higher stability and faster response performance while ensuring sensitivity, which indicates that the Kretschmann grating type SPR sensor structure optimization design method established in the present invention is effective and has certain practicality.

Claims

1. A method for optimizing the performance of an SPR grating sensor based on IMPA, characterized in that: (1) Determine the sensor structure and its parameters; (2) Determine the constraints in the sensor structure design; (3) Determine the objective function of the sensor structure optimization design; (4) Initialize the model parameter solution; (5) IMPA iteratively solve the optimization model. The specific steps include the following five steps: Step 1: Determine the sensor structure and its parameters; The traditional Kretschmann type SPR grating sensor structure is usually composed of four parts: "substance to be measured-metal film-base layer-optical path system". However, with the in-depth study of SPR grating sensors, it is found that the surface of ordinary single metal film layer SPR grating sensors is relatively active and easily reacts with the solution to be measured during the detection process, thereby reducing the sensitivity of the sensor. To address this problem, this paper attempts to improve the performance of the SPR sensor by adding other film materials as auxiliary film layers on the surface of the metal film in the structural design of the SPR grating sensor. The structure of the SPR grating sensor studied consists of five parts: "substance to be measured-biofilm-metal film-base layer-optical path system"; Based on the five-layer structure of "substance to be tested - biofilm - metal film - substrate layer - optical path system", the specific parameters we need to design are as follows: Biofilm thickness d b ; Metal film thickness d m ; Grating teeth number N g ; The convex-concave spacing of the grating is Λ; ·Incident wavelength λ; Angle of incidence θ; Step 2: Determine the constraints in the sensor structure design; 1) Physical restrictions Biofilm thickness d b :The thickness of the biofilm has a certain reasonable range. If it is too thick, it may affect the molecular diffusion rate and sensor response time. If it is too thin, it may not be able to effectively fix the molecules to be tested or provide sufficient physical and chemical modification. Assuming that the lower limit of the biofilm thickness is d bmin , the upper limit is d bmax ,but: d bmin ≤d b ≤d bmax ; Metal film thickness d m :The thickness of the metal film also has its limit. If it is too thin, it may not be able to effectively excite the surface plasma resonance. If it is too thick, it will increase the cost and may cause signal attenuation. Let the lower limit of the thickness of the metal film be d mmin , the upper limit is d mmax ,but: d mmin ≤d m ≤d mmax ; Grating teeth number N g : The number of grating teeth must be a positive integer and is limited by factors such as the grating size. Assume that the lower limit of the number of grating teeth is N gmin , The upper limit is N gmax , and is an integer, then: N gmin ≤N g ≤N gmax , The convex-concave spacing of the grating Λ: The convex-concave spacing must be within a reasonable range to ensure the optical performance of the grating, and is related to the number of grating teeth N. g There is an inverse relationship, and the lower limit of the concave-convex spacing is Λ min , the upper limit is Λ max ,but: L min ≤Λ≤Λ max ; Incident wavelength λ: The incident light wavelength must be within the wavelength range that the sensor can detect, assuming the lower limit is λ min , with an upper limit of λ max ,but: l min ≤λ≤λ max ; ·Incident angle θ: The incident angle usually has to be within a certain angle range to effectively stimulate the SPR effect, that is, to reach a specific exit angle θ at a specific wavelength λ m , can be adjusted by the grating equation: (sin(θ)+sin(θ) m ))Λ=mλ; Where m is the diffraction order, θ m is the diffraction angle, Λ is the concave-convex spacing; 2) Environmental restrictions Effect of temperature T: Assume that the sensor is within the temperature variation range (T min ≤T≤T max ) when working, the signal change rate is required Cannot exceed a certain threshold ΔS Tmax ,but: In the formula, S is the sensitivity, ΔT represents the change in temperature, and ΔS represents the change in sensitivity; Effect of humidity H: Assume that the sensor is within the temperature variation range (H min ≤H≤H max ) when working, the signal change rate is required Cannot exceed a certain threshold ΔS Hmax ,but: In the formula, S is the sensitivity, ΔH represents the change in humidity, and ΔS represents the change in sensitivity; Step 3: Determine the objective function of sensor structure optimization design; 1) Maximize sensitivity S Assume that the concentration of the substance to be measured is C, ΔS represents the change in signal response, ΔC represents the change in concentration, and define In the sensor setting, we often seek to maximize the sensitivity, that is, max(S). However, in order to avoid the problem of negative sensitivity caused by algorithm convergence, we perform a logarithmic transformation on the sensitivity, that is, max(ln(S+1)). At the same time, considering that there may be certain measurement errors in the signal response change and concentration change in actual measurement, the weight ω is introduced. s1 and ω s2 To adjust the importance of the two, the objective function can be further optimized as follows: max(ω s1 ×ln(ΔS+1)-ω s2 ×ln(ΔC+1)); 2) Stability optimization Stability can be measured by the fluctuation degree of the signal when the sensor works in different environmental conditions or for a long time; let σ be the standard deviation of signal fluctuation. In order to make the stability index have a good balance with the sensitivity and other indicators in the objective function, the standard deviation is normalized, that is, (σ max is the maximum standard deviation determined based on prior knowledge or experimental data), and the objective function is min(σ n ); At the same time, in order to avoid the indicator having too large a weight in the objective function due to a small standard deviation in the early stage of optimization, an exponential form can be used, namely: 3) Optimize response speed Response speed r It is defined as the time required from the time when the substance to be tested contacts the sensor to the time when the signal reaches a stable value. In order to make the response speed index and other indicators have a unified trend in the objective function (the smaller the better), the reciprocal of the response time can be taken. The objective function is set as: 4) Establishment of comprehensive objective function Taking the above factors into consideration, the objective function of multi-objective optimization is constructed as follows: In the formula, ω1, ω2, and ω3 are the weights corresponding to sensitivity, stability, and response speed, respectively; Step 4: Initialization of model parameter solution; After constructing the above complete "decision variables-constraints-objective function" operations research model, we consider using the Marine Predators Algorithm (MPA) to solve the optimal sensor parameters of the model. MPA is a new meta-heuristic optimization algorithm proposed by Afshin Faramarzi et al. in 2020; the algorithm is inspired by the foraging strategies of marine predators, especially Lévy and Brownian motions, and the optimal encounter rate strategy between predators and prey; MPA simulates the theory of survival of the fittest in the ocean and has strong evolutionary ability, fast search speed and strong optimization ability; the optimization process of MPA includes three main stages to simulate the movement behavior of predators and prey throughout their life cycle; each stage considers different speed ratios, so that MPA can more accurately simulate the movement relationship between predators and prey and provide better performance in the optimization process; MPA first uses a set of random solutions as the initial search space, and generates the initial solution by the formula according to the range of the search space restriction; The optimal solution is selected as the top predator, and a matrix named Elite is constructed. This matrix array monitors the search process and searches for prey based on the location information of the prey. The second matrix is ​​the Prey matrix, which has the same dimension as the Elite matrix. The predator is updated based on this matrix. At the end of each iteration, if a predator with a higher fitness value appears, the current top predator will be replaced and the elite matrix will be updated accordingly. Similar to most metaheuristics, MPA is a population-based method where the initial solutions are randomly distributed over the search space. In the first step: X0=X min +rand(X max -X min ); Where X max and X min are the lower and upper bounds of the variable, and rand is a uniform random vector in the range of 0 to 1; However, the standard MPA random initialization may cause the initial position of the population to be relatively scattered, so the individuals in the population may need to go through multiple iterations to gradually gather and optimize, thereby increasing the number of iterations and time complexity of the algorithm; the use of Logistics-Tent chaos initialization ensures that the initial solution is generated more uniformly. The specific improvements are as follows: In the formula, r∈(0,4), p0∈(0,1), the chaotic sequence {P n }, and then use this chaotic sequence to initialize the position of the population as follows: X0=X min +{P n }(X max -X min ); During the initialization process, the optimal solution is selected as the top predator to construct an elite matrix; during the search process, the solutions after the survival of the fittest are named top predators, and they form a matrix, also called the elite matrix; the array of this matrix supervises the search and finds the prey based on the information of the prey's location; Where X is the top predator vector, which is replicated n times to construct the elite matrix; n is the number of search agents, and d is the dimension; it is worth noting that both predators and prey are considered search agents; At the end of each iteration, if the top predator is replaced by a better predator, the elite will be updated; Another matrix with the same dimension as the "elite matrix" is called the prey matrix, and the predator updates its position based on this matrix; in short, initialization creates the initial prey, and the predator constructs the elite; the prey is shown in the following formula: X i,j represents the jth dimension of the i-th prey, n is the number of search agents, and d is the dimension. The optimization process is directly related to these two matrices; Step 5: IMPA iteratively solves the optimization model; Before introducing the steps of the proposed IMPA (Improved Marine Predators Algorithm, IMPA) algorithm, it is necessary to know the mathematical models of the two random motions of Lévy and Brownian motion; then the optimal foraging, movement and memory management strategies of marine predators are introduced one by one and the corresponding mathematical models are established, and the overall process of the MPA algorithm is gradually given; 1) Brownian motion Standard Brownian motion is a random process in which their step sizes are drawn from a probability function defined by a normal (Gaussian) distribution with zero mean (μ = 0) and unit variance (σ 2 =1); the probability density function (PDF) of the motion at point x is as follows: 2) Lévy Flight A Lévy flight is a random walk whose step size is determined by the probability function defined by the Lévy distribution (power-law tail): L(x j )≈|x j | 1-α ; where x j is the step size, 1<α≤2 is the power law exponent, and the probability density function of the Levy stable process in integral form is as follows: Where α defines the distribution index and the scaling properties of the control process, and γ is the selected scaling unit; the integral in the formula has an analytical solution only in a few cases; when α = 2, it represents a Gaussian distribution, and when α = 1, it represents a Cauchy distribution; the solution of the integral in the formula generally requires the use of a series expansion method only when the x value is large, as shown in the following formula: Where Γ is the gamma function, α is an integer, and Γ(1+α)=α; Mantegna proposed an exact and fast algorithm for generating Levy stable processes under exponential distribution (α) in the range of 0.3 and 1.99; the Mantegna method is used to generate random numbers based on the Levy distribution as follows: In Lévy flight, the area of ​​interest is traced primarily by small steps associated with long jumps, whereas in contrast, Brownian motion can cover the area of ​​interest with more uniform and controlled steps; 3) Early stage of iteration At high speeds or when the predator moves faster than the prey, the predator does not move at all. This happens in the initial iterations of the optimization, which is the key to exploration; the mathematical model for this stage is: in is a vector containing random numbers from a normal distribution, representing Brownian motion, represents moment multiplication; P = 0.5 is a constant, R is a standard random number, Iter represents the current number of iterations, Max_ Iter represents the maximum number of iterations, represents the step length of each movement of prey i; 4) Mid-term iteration At unit speed ratio or when both predator and prey move at the same speed, the simulation is that both predator and prey are searching for prey. At this stage, both exploration and exploitation are important. Therefore, half of the group is designated for exploration and the other half is designated for exploitation. At this stage, the prey is responsible for exploitation and the predator is responsible for exploration. According to this rule, when the prey is doing Levy motion, the best strategy for the predator is Brownian motion. For the first half of the population: in is a random number vector based on Levy distribution, representing Levy motion; The product of the prey and the prey simulates the movement of the prey in the Levy way; Iter represents the current number of iterations, Max_ Iter represents the maximum number of iterations, represents the step length of each movement of prey i; For the second half of the population: in is an adaptive parameter that controls the predator's moving step size; 5) Late iteration At low speed ratios or when the predator moves faster than the prey; the predator's optimal strategy is Levy movement; This happens at the end of the optimization process, which is mainly related to high development capabilities and can be described as follows at low speed ratios: and The multiplication of simulates the Lévy motion of the predator, while The step length is increased in the position to simulate the movement of the predator; Iter represents the current number of iterations, Max_ Iter represents the maximum number of iterations, represents the step length of each movement of prey i; 6) Vortex formation and the influence of FADs Another point that causes changes in the behavior of marine predators is environmental issues, such as vortex formation or fish aggregating devices (FADs) effects. FADs can be considered as local optima that capture these points in the search space. Considering these longer jumps during the simulation process can avoid stagnation in local optima; therefore, the FADs effect is expressed in mathematical form as: Where FADs = 0.2 is the probability that FADs affects the optimization process, U is a binary vector, the array includes 0 and 1, this is to generate a random vector in [0,1], if the array is less than 0.2, it is changed to 0, if the array is greater than 0.2, its array is changed to 1, r is a uniform random number in [0,1], X max and X min are the dimensional vectors of the lower and upper bounds, respectively, and the subscripts r1 and r2 represent the random index of the prey matrix; 7) Memories in the Ocean Marine predators have a good memory to remind them of places where they have successfully foraged before, and this ability is simulated by storing it in the MPA; after updating the prey and realizing the FADs effect, the fitness of this matrix is ​​evaluated to update the elite position; the fitness of each solution in the current iteration is compared with its solution in the previous iteration, and if the fitness of the current solution is greater, the current solution will replace the old solution; as the iterations go by, this process improves the quality of the solution and also simulates the predator returning to the area where the prey is abundant and successfully foraging; 8) Multi-subgroup improvement strategy However, for the problems of low solution accuracy and easy to fall into local optimum in the standard MPA algorithm, a multi-subpopulation improvement strategy is proposed; after the population of marine predators is iterated, each population is sorted according to the size of the obtained population fitness value, and the population is divided into 3 sub-populations. The individuals with fitness values ​​ranked in the top N / 3 are high-quality populations, those with fitness values ​​between N / 3 and 2N / 3 are classified as follower populations, and the rest are classified as low-quality populations; In a high-quality population, excellent individuals influence each other in two different dimensions by learning from each other. This learning process helps excellent individuals to conduct reinforcement learning based on their own experience, making them more focused on areas with a higher probability of optimal solutions, thereby accelerating the convergence speed of the algorithm and improving the accuracy of optimization. is the new position of the individual in the high-quality population, P1 is a standard random number, and are any two dimensions in an individual; For the follower population, reverse learning strategy is used to update: improve the algorithm's ability to escape from local optimality; is the new position of the individual in the follower population; X min and X max are the upper and lower limits of the problem dimension; is the new position of the individual in the low-quality population, and They are the global optimal solution and the solution with average fitness, r1 and r2 are standard random numbers. The global optimal solution participates in the update of the low-quality population, which speeds up the convergence of the algorithm. In summary, after multiple updates of the fitness function, the model will output the optimal sensor parameter design values.

Images