Extended Kantorovich method for local stresses in composite laminates upon polynomial stre

A regularizing L-curve Lanczos method forunderdetermined linear systemsS.Morigi,F.Sgallari *,1Department of Mathematics,CIRAM,University of Bologna,Piazza di Porta San Donato 5,I-40127Bologna,ItalyAbstractMany real applications give rise to the solution of underdetermined linear systems of equations with a very ill conditioned matrix A ,whose dimensions are so large as to make solution by direct methods impractical or infeasible.Image reconstruction from projections is a well-known example of such systems.In order to facilitate the com-putation of a meaningful approximate solution,we regularize the linear system,i.e.,we replace it by a nearby system that is better conditioned.The amount of regularization is determined by a regularization parameter.Its optimal value is,in most applications,not known a priori.A well-known method to determine it is given by the L-curve approach.We present an iterative method based on the Lanczos algorithm for inexpensively evaluating an approximation of the points on the L-curve and then determine the value of the optimal regularization parameter which lets us compute an approximate solution of the regularized system of equations.Ó2001Elsevier Science Inc.All rights reserved.Keywords:Ill posed problems;Regularization;L-curve criterion;Lanczos algorithm1.IntroductionIn this article,we describe a new regularizing iterative method for the so-lution of ill conditioned underdetermined linear systems of equationsA x b ;A P R m Ân ;x P R n ;b P R m ;m <n : 1/locate/amcApplied Mathematics and Computation 121(2001)55±73*Corresponding author.E-mail address:sgallari@dm.unibo.it (F.Sgallari).1This research was supported by the University of Bologna,funds for research topics.0096-3003/01/$-see front matter Ó2001Elsevier Science Inc.All rights reserved.PII:S 0096-3003(99)00262-356S.Morigi,F.Sgallari/put.121(2001)55±73Linear systems of this kind arise,for example,in the reconstruction of images from projections.In this application,each column corresponds to one pixel of the image to be reconstructed,and each rowto a``ray''from an emitter to a detector.The matrices that arise are typically large and have many more col-umns than rows.Further,they are generally very ill conditioned,because the linear system is obtained by discretizing an ill posed problem,the inversion of a Radon transform.The ill conditioning makes a solution of(1)very sensitive to perturbations in the right-hand side vector.Such perturbations arise,e.g.,from measurement errors[1,2].In order to be able to compute a meaningful approximate solution of(1),the system has to be regularized,i.e.,the linear system(1)has to be replaced by a linear system of equations with a not very ill conditioned matrix,such that the unique solution of this system yields an acceptable approximate solution of(1). We introduce the matrixB: AA T:The solution of the minimal norm of(1)can be written as x: A T y,where y P R m solves the linear systemB y b: 2 We regularize this system by adding a multiple of the identity matrix to B and obtainB k I y b; 3 where k>0is a regularization parameter.We denote the solution of(3)by y k. The smaller the value of k,the closer the solution y k of(3)is to the minimal norm solution of(2),and the more ill conditioned the matrix B k I.We wish to®nd a value of k in(3),such that the matrix B k I is not very ill conditioned, and the solution of(3)is an acceptable approximate solution of(2).Generally, it is not known a priori for which values of k these requirements are satis®ed. The approximate solution x k A T y k of(1),obtained by®rst solving(3)for y k,can also be computed by standard Tikhonov regularization of system(1). This approach,however,as has been shown in[3],requires more computer storage and arithmetic work when the linear system has many more columns than rows.Let us consider the minimization problemfk B yÀb k2 k k A T y k2g; 4 miny P R mwhere kÁk denotes the Euclidean norm here and in the sequel.The solution y k of(4)satis®es the equationB T B k AA T y B T b: 5S.Morigi,F.Sgallari/put.121(2001)55±7357 The comparison of(3)±(5)shows the equivalence of problem(3)and(5) when B is a nonsingular matrix and this justi®es what follows.Let components of the right-hand side vector b be contaminated by errors, e.g.,due to inaccurate measurements,and b exact and e denote the unperturbed right-hand side and the error,respectively.If the norm of the error e is known, then the value of k can be chosen so that the associated solution satis®es the Morozov discrepancy principle[3].However,in many applications the norm of e is not explicitly known,and other approaches to determine a value of the regularization parameter k have to be employed.We will show that the plot of the curve(k B y kÀb k;k A T y k k)is shaped roughly like an``L'',and the optimal value of the regularization parameter k corresponds to the``vertex''point of the L-curve.This idea of L-curve was observed by Lawson and Hanson[4]to determine the regularizing parameter k in the Tikhonov approach,while the choice of the optimal regularization pa-rameter as the vertex point was suggested recently in[5,6],and successively considered in[7].The computation of the L-curve is quite costly for large problems;in fact, the determination of a point on the L-curve requires that both the norm of the regularized approximate solution and the norm of the corresponding residual vector be available.Therefore,usually only a fewpoints on the L-curve are computed,and these values,which we will refer to as the discrete L-curve,are used to determine a value of the regularization parameter.In[8],a newmethod for approximately determining the points on the L-curve is proposed.They showed how rectangular regions that contain points on the L-curve(L-ribbon)can be determined fairly inexpensively for several values of the regularization parameter without having to solve the corresponding regularized minimization problem but using Gauss and Gauss±Radau quadr-ature rules.When the L-ribbon is narrow,it is possible to infer the approximate location of the vertex of the L-curve from the shape of the L-ribbon.In order to solve problem(4),or equivalently(3),we cannot follow the Calvetti et al.[8]approach because their assumptions are not satis®ed due to the term A T in the regularizing functional.We present an iterative procedure that determines successive encapsulated intervals containing the optimal regularization parameter corresponding to the vertex of the L-curve.At each step,some newpoints on the L-curve are computed by a fewsteps of the Lanczos bidiagonalization algorithm,while a new interval for k is de-termined by a minimum distance from the local origin criterion.The®nal value obtained for the regularization parameter is used to compute an approximation of the solution of the regularized system(3).This article is organized as follows.In Section2,we give a brief descrip-tion of properties of the L-curve and we show how to use it to solve prob-lem(4).Section3explains the relation between the bidiagonalization Lanczosalgorithm and the approximated L-curve and howto update the information extracted from the Lanczos process to compute a newpoint.Section4de-scribes the regularizing algorithm for the solution of the linear system(3). Section5contains several computed examples which illustrate the performance of the algorithm.2.L-curve for the choice of regularization parameter valueA recent method to choose the regularization parameter is the so-called L-curve.The L-curve is a parametric plot of q k ;g k ,where g k and q k measure the size of the regularized solution and the corresponding residual. The underlying idea is that a good method for choosing the regularization parameter for discrete ill posed problems must incorporate information about the solution size in addition to using information about the residual size.This is indeed quite natural,because we are seeking a fair balance in keeping both of these values small.The L-curve has a distinct L-shaped corner located exactly where the solution y k changes in nature from being dominated by regulariza-tion errors(i.e.,by over smoothing)to being dominated by the errors in the right-hand side.Hence,the corner of the L-curve corresponds to a good bal-ance between minimization of the two quantities,and the corresponding reg-ularization parameter k is a good choice.In particular the continuous L-curve for our problem consists of all the points q k ;g k for k P 0;I ,where q k k B y kÀb k and g k k A T y k k. In what follows,we present some important properties of the continuous L-curve and a method to compute a good approximation of the discrete L-curve. We®rst need to introduce the following singular value decomposition for matrix A U R V T.Here U P R mÂm and V P R nÂn are orthogonal matrices, R R m0 P R mÂn,and R m diag f r i g,where the singular values r i are ordered as0<r1T r2TÁÁÁT r m.It is straightforward to see that the solution of problem(4)can be written asy kX mi 1u T i bk r2iu i: 6Using these notations,the L-curve has the following properties:Proposition1.Let y k denote the solution to(4).Then,k A T y k k is a monotonically decreasing function of k B y kÀb k.Proof.The fact that the k A T y k k is a decreasing function of k B y kÀb k follows immediately from the following expressions:k A T y k k2X mi 1r i u T i br2i k2; 758S.Morigi,F.Sgallari/put.121(2001)55±73k B y k Àb k 2X m i 1kr 2i k u T ib 2; 8which are easily derived from (6).ÃIf a fewreasonable assumptions are made about the ill posed problem,then it is possible to characterize the L-curve behaviour analogously to that done in [5,6]for the Tikhonov regularization.The three assumptions made are:1.The coe cients j u T i b exact j on average decay to zero faster than the r i ;2.The perturbation vector e is essentially ``white noise''with standard devia-tion r ;3.The norm of e is bounded by the norm of b exact .Under Assumptions 1±3,the L-curve q k ;g k is shown to exhibit a corner behaviouras a function of k ,and the corner appears approximately at r m p ;k A T y 0k .Here,y 0is the unregularized solution to the unperturbed problem (2).The larger the di erence between the decay rates of j u T i b exact j and j u Ti e j ,the more distinct the corner will appear.Following [6],we summarize the reasoning that leads to the previous characterization.We consider the behaviour of the L-curve for an unperturbedproblem with e 0.For k (r 21,w e have r i = r 2i k %1=r i ,so that y k %y 0,andk B y k Àb exact k 2%X m i 1k r 2i u T i b exact 2T k r 21k b exact k 2: 9 Hence,for small values of k ,the L-curve is approximately a horizontal line atk A T y k k k A T y 0k .As k increases,it follows from (7)that k A T y k k starts to de-crease,while k B y k Àb k still grows.Now,consider the L-curve associated with the perturbation e of the right-hand side.The corresponding solution y k given by (6)with b replaced by e ,satis®es (from Assumption 2)k A T y k k 2 X m i 1r i u T i e r 2i k 2%r 2X m i 1r i r 2i k 2; 10 k B y k Àb k 2X m i 1kr 2i k u T ie 2%r 2X m i 1kr 2i k2:11For k (r 21,this L-curve is approximately a horizontal line at k A T y k k %m p r =r 1,and it starts to decrease towards the abscissa axis for much smaller values of k than the L-curve for b exact (see Assumption 3).Moreover,we see that as k increases,k B y k Àb k becomes almost independent of k ,while k A T y k k isS.Morigi,F.Sgallari /put.121(2001)55±7359dominated by a fewterms and becomes about r2=kp.Hence,this L-curve soonbecomes almost a vertical line at k B y kÀb k%rmpas k3I.The actual L-curve for a given problem,with a perturbed right-hand side b b exact e,is a combination of the above two special L-curves.For small k, the behaviour is completely dominated by contributions from e,whereas,for large k it is entirely dominated by contributions from b exact.In between there is a small region where both b exact and e contribute,and this region contains the L-shaped corner of the L-curve.In actual computations,the values q k and g k are usually not available; hence,we search for good approximations.Let y k k be an approximation of the solution y k of(3),and introduce the residual vectorr k k: bÀ B k I y k k 12 and the quantitiesq k k : k B y k kÀb k; 13g k k : k A T y k k k: 14 The following proposition is concerned with how well the functions q k and g k approximate q and g,respectively.Proposition2.Let q k k and g k k be defined by(13)and(14),respectively. Thenj q k k Àq k j T k r k k k; 15 j g k k Àg k j T k r k k k: 16 Proof.It follows from the de®nition of q and q k thatj q k k Àq k j T k B y k kÀy k k k B y k kÀ B k I À1b kk B B k I À1r k k k T k B B k I À1kk r k k k;and in viewof k B B k I À1k T1,the®rst part of the proposition follows.For the second part,we have from the de®nition of g and g kj g k k Àg k j T k A T y k kÀy k k k A T y k kÀ B k I À1b kk A T B k I À1r k k k T k A T B k I À1kk r k k k:Using the singular value decomposition of the matrix A U R V T and B U R2U T,we obtain that k A T B k I À1k T1,then the last part of the proposition follows.Ã60S.Morigi,F.Sgallari/put.121(2001)55±73Thus,we can obtain a good approximation of the discrete L-curve by computing q k k ;g k k ,and quantifying the errors by the length of the re-sidual vector.3.The Lanczos algorithm and regularizationIn this section,we brie¯y review the bidiagonalization Lanczos process[9], and,we explore its relation with the functions q k k and g k k introduced in Section2.Given a symmetric matrix B AA T and a nonvanishing vector b,if we ex-ecute k iterations of the Lanczos bidiagonalization Algorithm1in exact arithmetic,we get the orthogonal mÂk matrixV k v0;v1;...;v kÀ1 17 and the orthogonal nÂk matrixW k w0;w1;...;w kÀ1 : 18 The matrices V k and W k are connected by the two equationsAW k V k L k d k v k e T k; 19A T V k W k L T k; 20 where L k denotes the kÂk lower bidiagonal matrixL k:c0d1c1......d kÀ2c kÀ2d kÀ1c kÀ12666664377777521and e k is the k th column of the identity matrix.Algorithm1(Lanczos bidiagonalization algorithm).Input:b P R m n f0g,B P R mÂm,0<k<m;Output:f d j g k j 1;f c j g kÀ1j 0;f v j g kÀ1j 0;v0: b=k b k;s0: A T v0;c0: k s0k;w0: s0=c0;for j 2;3;...;k doS.Morigi,F.Sgallari/put.121(2001)55±7361r jÀ1: A w jÀ2Àc jÀ2v jÀ2;d jÀ1: k r jÀ1k;v jÀ1: r jÀ1=d jÀ1;s jÀ1: A T v jÀ1Àd jÀ1w jÀ2;c jÀ1: k s jÀ1k;w jÀ1: s jÀ1=c jÀ1;end jLet T k 1;k denote a tridiagonal matrixT k 1;k:a0b1b1a1b2.........b kÀ2a kÀ2b kÀ1b kÀ1a kÀ1b k2666666666437777777775P R k 1 Âk: 22It is nowstraightforw ard to identify T k;k L k L T k if a j d2j c2j, b j d j c jÀ1;j 1;...;kÀ1,and a0 c20.By combining the two relations(19) and(20),we getBV k V k L k L T k d k c kÀ1v k e T k V k 1T k 1;k; 23 with V k 1 v0;v2;...;v k ,such that V T k 1V k 1 I k 1and V k 1e1 b=k b k. Throughout this article,I k 1denotes the identity matrix of order k 1,I k 1;k its leading principal k 1 Âk submatrix.It follows from Proposition2that in order to have a good approximation of the L-curve points q k ;g k ,we have to choose k su ciently large and in-crease it until k r k k k is under a given tolerance.Now,we want to show how,by means of Algorithm1,it is possible to compute the residual k r k k k for increasing values of k in a straightforward manner.For notational simplicity,we will consider in the sequel the T k 1;k matrix even if,from a computational point of view,only the necessary elements are determined starting from L k.We determine an approximate solution of(3)of the form y k k V k z by solving the minimization problemmin z P R k k B k I V k zÀb k minz P R kk V k 1 T k 1;k k I k 1;k zÀV k 1e1k b kkminz P R kk T k 1;k k I k 1;k zÀe1k b kk: 24Let us introduce the QR-factorization of the tridiagonal matrix~Tk 1;k : T k 1;k k I k 1;k ~Q k 1~R k 1;k; 2562S.Morigi,F.Sgallari/put.121(2001)55±73where ~Q k 1P R k 1 Â k 1 ,~Q T k 1~Q k 1 I k 1and ~R k 1;k P R k 1 Âk has an upper triangular leading principal k Âk submatrix,denoted by ~R k ,and a vanishinglast row.It follows thatmin z P R kk B k I V k z Àb k min z P R kk ~R k 1;k z À~Q T k 1e 1k b kk j e T k 1~Q T k 1;k e 1jk b khas the solutionz k k: ~RÀ1kI T k 1;k ~Q Tk 1e 1k b k :26Thus,the associated approximate solution of (3)is given byy k k : V k ~R À1k I T k 1;k ~Q T k 1e 1k b k ;27and the corresponding residual vector (12)has the normk r k k k j e T k 1~Q Tk 1e 1jk b k :28We note for future reference that,when y k k is given by (27),the quantities q kand g k de®ned by (13)and (14),respectively,can be evaluated asq k k k T k 1;k z k k Àe 1k b kk ;29 g k k k L T k z k k k ;30where z k k is given by (26).When increasing k ,the QR-factorizations of the matrices ~Tk 1;k are updated.We outline the computations required for this.As the matrix ~Tk 1;k is tridiag-onal,it is convenient to compute its QR-factorization by Givens rotations [10,Chapter 5].Let G ik 1be a Givens rotation of order k 1that rotates the i ;i 1 coordinate planes so that the i th subdiagonal element of the matrix G i À1 k 1;...;G 1 k 1~T k 1;k vanishes.Then,the factors in the QR-factorization (25)of ~Tk 1;k are given by ~Q k 1: G k k 1;...;G 1 k 1T ;~R k 1;k : G k k 1;...;G 1 k 1~T k 1;k :31Taking an additional step with the Lanczos algorithm yields the matrixwhose QR-factorization can be determined by updating factors (31)as follows.For 1T i T k ,let G ik 2denote the Givens rotation of order k 2whose leadingprincipal submatrix of order k 1is G ik 1.Then,32where Ãdenotes a matrix element that may be nonvanishing.The upper tri-angularmatrixin the QR-factorization ~Tk 2;k 1 ~Q k 2~R k 2;k 1is obtained by multiplying matrix (32)from the left by a Givens rotation G k 1k 2that annihilates the k 2;k 1 entry b k 2.Thus,andQ T k 2e 1 G k 1 k 2~Q T k 1e 102435:64S.Morigi,F.Sgallari /put.121(2001)55±73We point out that in order to compute ~Rk 2;k 1,given ~R k 1;k ,we only need to apply the Givens rotations G k À1 k 2,G k k 2and G k 1k 2,in this order,to the lastcolumn of ~Tk 2;k 1.We note that q k k and g k k ,for ®xed k ,can be evaluated for di erent values of k without recomputing the Lanczos decomposition (23).Each newvalue of k yields a newtridiagonal matrix ~Tk 1;k ,whose QR-factorization,see (25)and (31),has to be computed.The computation of this factorization is inexpensive;it requires only O k arithmetic operations,as it can be carried with k Givens rotations.4.A regularizing algorithm based on the Lanczos processIn this section,our method for the solution of underdetermined ill condi-tioned linear systems (1)is described.The evaluation of points on the L-curve requires that the exact solution of the linear system (3)is available.However,as already mentioned,computation of the exact or a highly accurate solution of (3)for each value of k is expensive,and therefore,we replace the equations for q k and g k by (29)and (30),respectively.It follows from Proposition 2and (28)that the norm of the re-sidual vector r k k is a computable inexpensive bound for the error introduced when approximating q k and g k by q k k and g k k .Algorithm 2computes the regularization parameter and the corresponding approximate solution of the regularized system (3).Algorithm 2(Regularizing L-curve Lanczos iteration ).Input :A P R m Ân ,b P R m , d , r P R ,[k i ,k f ],k i ;k f P R ,N ;k val P N ;Output :k opt ,regularized approximate solution x kopt k P R n ;%Initialization %.k : k val;j 0;k 1i k i ,k 1f k f ;k 0opt : 0;repeat.j j 1;.k j s : k j f Àk ji = N À1 ,%L-curve computation phase %.Perform k Lanczos steps by applying Algorithm 1;for k `: k j i to k j f step k js do.Compute QR-factorization of ~Tk 1;k ;.Evaluate k r k k ;while (k r k k > r )do .k k 1;.Apply one more Lanczos step;end ;S.Morigi,F.Sgallari /put.121(2001)55±736566S.Morigi,F.Sgallari/put.121(2001)55±73 .Compute point P` q k k` ;g k k` ;end;%Detection corner phase%.Detect the smallest interval that contains the closest point P j to the origin; .Compute the associated k j opt;until(j k jÀ1optÀk j opt j T d j k jÀ1opt j);%Resolution phase%.Evaluate the®nal solution of A x b by solving x k opt k A T y k opt k.The initialization step needs the de®nition of an initial interval k i;k f for the regularizing parameter,the number of required points(N)on the discrete L-curve,and the initial dimension k val of the Krylov space.We will refer to each step j of the repeat loop as re®nement step j.Each of these steps consists of an L-curve computation phase,which calculates N even spaced points on the discrete L-curve belonging to the interval k j i;k j f ,and of a detection phase for determining the corner of the L-curve at step j.Let us consider the L-curve computation phase in more detail.Given an initial k value(k val),k Lanczos steps are preliminarily performed. Using the relations between d and c entries of L k matrix and a and b entries of T k 1;k the~Q~R factorization(25)is then computed and further Lanczos steps are performed until the residual norm,evaluated by relation(28),does not satisfy the given precision r.The approximate point P` q k k` ;g k k` of the discrete L-curve is then given by(29)and(30).To determine the next point P`of the L-curve the algorithm uses for k and L k the previous values,then,using(25),a newQR-factorization is computed and consequently the residual is available.The corner detection phase considers the L-curve resulting from the previ-ous phase and determines the smallest interval containing the k value corre-sponding to the closest point P j opt to a local origin to a given tolerance d.The local origin coordinates are computed in the®rst re®nement step as min`q k k` and min`g k k` .Note that during each re®nement step,we can also verify whether the given j th interval contains the corner point P j opt.In fact,if the k value cor-responding to the returned point P opt j is equal to one of the extreme points of the interval k j i;k j f ,then the k value associated with the corner point is not contained in the present interval.Thus,it will be necessary to expand the interval and consequently to compute some additional points to complete the L-curve.The last step of Algorithm2evaluates the approximate solution y k opt k of(3) by(27),and the regularized approximate solution of the original problem(1) by x k opt k A T y k opt k.As a®nal remark,we observe that each re®nement step can take advantage of the previous Lanczos step without recomputing all the values.5.Numerical resultsThe purpose of this section is to illustrate the performance of Algorithm2in some numerical examples obtained on a Sun Ultra workstation with a unit roundo of 2À52%2:2204Â10À16.We implemented our method to com-pute an approximation of the L-curve in FORTRAN77language,and we used MATLAB commercial software package[11]to obtain the exact computation of points on the L-curve.In order to obtain small and large scale test problems,we generated an m-by-n matrix A de®ned by its singular value decompositionA U R V T; 33 whereU IÀ2uu Tk u k22;V IÀ2vv Tk v k22; 34and the entries of u and v are random numbers.R is a diagonal matrix with the singular values r i eÀ0:2 iÀm .On the right-hand side of the linear system(1)to be solved b exact is a random vector,and e is the noise vector generated by®rst determining a vector with normally distributed random numbers with zero mean and variance one as components,and then scaling it to have the desired norm.We will refer to the quantity k e k=k b k as the noise level.All computations here reported correspond to small ill conditioned test problems m 100;n 200;r m 1:0;r1 2:5175Â10À9 in order to be able to perform qualitative comparisons with the MATLAB results.The results of each re®nement step of Algorithm2are summarized in Table 1which consider,respectively,the three di erent noise levels0:1;0:01and 0:001.The®rst column in the tables indicates the level of re®nement(LR),the second column labelled`` k i;k f ''shows the interval of the regularization pa-rameter k,the columns labelled``k opt e''and``k opt a''contain the optimal regular-ization parameter determined in that re®nement interval,respectively,by MATLAB(exact)and by Algorithm2(approximated).The last two columns report the distance of the corner from the local origin,chosen for a given noise level(distance),and the number of evaluation points in the re®nement step (NP).The precision`` r'',used in Lanczos procedure,is reported in the captions. The dimension k of the Krylov subspace is not reported in the tables because it maintains a value of®ve for all the computations.The only exception is for some extremely small k values(e.g.,1:0Â10À9;1:0Â10À10)when a high pre-cision r is used in the Lanczos procedure.To avoid the use of k greater than 100,we suggest either to consider a smaller interval for k,or to use,for the computation related to the critical extreme value,a higher precision r.S.Morigi,F.Sgallari/put.121(2001)55±7367Table1Test problemLR k i;k f k opt e k opt a Distance NP(a)Noise level0:1, r 0:11 1:0Â10À9;2:0 3:0Â10À13:0Â10À15:963200Â10À1202 4:0Â10À1;2:0Â10À1 2:5Â10À12:5Â10À15:956456Â10À1203 2:4Â10À1;2:6Â10À1 2:5Â10À12:5Â10À15:956456Â10À1201 1:0Â10À9;2:0 2:5Â10À12:5Â10À15:956456Â10À1200(b)Noise level0:01, r 0:051 1:0Â10À9;2:0 1:0Â10À11:1Â10À1 1.211257202 1:0Â10À9;2:0Â10À1 1:0Â10À21:0Â10À28.901137Â10À1203 1:0Â10À9;2:0Â10À2 1:0Â10À31:0Â10À37.939300Â10À1204 1:0Â10À9;2:0Â10À3 1:3Â10À31:0Â10À37.930785Â10À1201 1:0Â10À9;2:0 1:0Â10À21:0Â10À28.901137Â10À12002 1:0Â10À9;2:0Â10À2 9:0Â10À41:0Â10À37:930785Â10À1200(c)Noise level0:001, r 0:0051 1:0Â10À10;2:0 1:0Â10À11:0Â10À1 2.232340202 1:0Â10À10;2:0Â10À1 1:0Â10À21:0Â10À2 1.715200203 1:0Â10À10;2:0Â10À2 1:0Â10À31:0Â10À3 1.306169204 1:0Â10À10;2:0Â10À3 1:0Â10À41:0Â10À49.563796Â10À1205 1:0Â10À10;2:0Â10À4 1:0Â10À51:0Â10À58.587058Â10À1206 1:0Â10À10;2:0Â10À5 7:0Â10À66:0Â10À68.572574Â10À1207 6:0Â10À6;8:0Â10À6 7:2Â10À65:8Â10À68.570727Â10À1201 1:0Â10À10;2:0 1:0Â10À21:0Â10À29.806748Â10À12002 1:0Â10À10;2:0Â10À2 1:0Â10À41:0Â10À48.646498Â10À12003 1:0Â10À10;2:0Â10À4 5:0Â10À67:0Â10À68.572574Â10À120068S. Morigi, F. Sgallari / Appl. Math. Comput. 121 (2001) 55±73。

DECface:An Automatic Lip-Synchronization Algorithm for Synthetic FacesKeith Waters and Thomas M.Levergood Digital Equipment CorporationCambridge Research LabDigital Equipment Corporation has four research facilities: the Systems Research Center and the Western Research Laboratory, both in Palo Alto, California; the Paris Research Laboratory, in Paris; and the Cambridge Research Laboratory, in Cambridge, Massachusetts.The Cambridge laboratory became operational in 1988 and is located at One Kendall Square, near MIT. CRL engages in computing research to extend the state of the computing art in areas likely to be important to Digital and its customers in future years. CRL’s main focus is applications technology; that is, the creation of knowledge and tools useful for the preparation of important classes of applications.CRL Technical Reports can be ordered by electronic mail. To receive instructions, send a mes-sage to one of the following addresses, with the word help in the Subject line:On Digital’s EASYnet:CRL::TECHREPORTSOn the Internet:techreports@This work may not be copied or reproduced for any commercial purpose. Permission to copy without payment is granted for non-profit educational and research purposes provided all such copies include a notice that such copy-ing is by permission of the Cambridge Research Lab of Digital Equipment Corporation, an acknowledgment of the authors to the work, and all applicable portions of the copyright notice.The Digital logo is a trademark of Digital Equipment Corporation.Cambridge Research LaboratoryTMOne Kendall SquareCambridge, Massachusetts 02139DECface:An AutomaticLip-Synchronization Algorithmfor Synthetic FacesKeith Waters and Thomas M.LevergoodDigital Equipment CorporationCambridge Research LabCRL93/4September23,1993AbstractThis report addresses the problem of automatically synchronizing computer generated faces with synthetic speech.The complete process is called DECface which provides a novel form of face-to-face communication and the ability to create a new range of talking personable synthetic characters.Based on plain ASCII text input,a synthetic speech segment is generated and synchronized in real-time to a graphical display of an articulating mouth and face.The key component of DECface is the run-time facility that adaptively synchronizes the graphical display of the face to the audio.Keywords:face-to-face,key-framing,anatomically-based,text-to-speech,viseme, phoneticsc Digital Equipment Corporation1993.All rights reserved.1 1IntroductionFrom an early age we are sensitive to the bimodal nature of speech,using cues from both visual and auditory modalities for comprehension.The visual stimuli are so tightly integrated into our perception of speech that we commonly believe that only the hearing-impaired lip-read[McG85].In fact,people with normal hearing use all available visual information that accompanies speech,especially if there is a degradation in the acoustical signal[MM86].Fluent speech is also emphasized and punctuated by facial expressions,thereby increasing our desire to observe the face of the speaker.Our goal is to use the expressive bandwidth of the human face in real-time synthetic facial models capable of interacting with and eliciting responses from the user.Although an ideal interface might eventually be a natural dialogue between humans and computers,we consider a subset of this larger goal:a technique to automatically synchronize mouth shapes to a real-time speech synthesizer.A computer-generated face has distinct advantages over images of real people, primarily because it is possible to create and control precise,repeatable facial actions.These faces suggest some unique and novel scenarios for presenting information,particularly where two-way interaction can be enhanced by listening rather than reading information.Examples of this type of interaction can be found in walk-by kiosks,ATM tellers,office environments,and videophones of tomorrow. In the near future we are going to see man-machine interfaces that mimic the way we interact face-to-face.A few years ago Apple Computer,Inc.produced a futuristic video called the The Knowledge Navigator,popularizing and advertising the notion of active agents.In the video,Phil,an active agent,performs a variety of tasks at the request of the user.In reality,no such system or environment exists. However it is worth noting that Phil,a head and shoulders image of a real actor,was the primary interface to the computer.Synthetic facial images also have potential in situations where information has to be presented in a controlled,unbiased fashion, such as interviewing.Furthermore,if the synthetic speech/face generator were combined with systems that perform basic facial analysis by tracking the focus of the user and analyzing the user’s speech,it would be possible to transform the computer from an inert box into a personable computer[CHK92].2Background and Previous WorkSome of thefirst images of animated speech were created by G.Demeny in1892 with a device he called the Phonoscope[Des66].The device mounted images of22BACKGROUND AND PREVIOUS WORK the lower face on a disk that rotated fast enough to exploit persistence of vision. Although the Phonoscope is nearly a hundred years old,it is the underlying process employed in animation today.Rather than using photographs,traditional animation relies on hand drawn images of the lips[TJ81].A sound track isfirst recorded,then an exposure sheet is marked with timing and phonetic information.The animator then draws a tailor-made mouth shape corresponding to a precise frame time.As one might expect,the whole process is extremely labor intensive.Thefirst attempts in computer-based facial animation involved key-framing, where two or more complete facial expressions are captured and in between frames computed by interpolation[Par72].The immense variety of facial expressions makes this approach extremely data intensive,and prompted the development of parametric models for facial animation.Parametric facial models[Par74,Par82] create expressions by specifying sets of parameter value sequences;for instance, by interpolating the parameters rather than direct key-framing.The parameters control facial features,such as the mouth opening,height,width,and protrusion. The limitations of ad hoc parameterized models prompted a movement towards models whose parameters are based on the anatomy of the face[Pla80,PB81, Wat87,Wai89,MPT88].Such models operate with parameters based on facial muscle structures.When anatomically based models incorporate facial action coding schemes as control procedures[EF77],it becomes relatively straightforward to synthesize a range of recognizable expressions.The geometric surface of the face is typically described as a collection of polygons and displayed using standard lighting models[Gou71,Pho76].Texture mapping of reflectance data acquired from photographs of faces or from laser digitizers provide another valuable technique for further enhancing the realism of facial modeling and animation[NHS88,Wil90,CT91].Even when the geometry is coarse,striking images can be generated[OTO87].To date,non-automated techniques are most commonly used to achieve lip-synchronization.This process involves recording the sound track from which a series of controlfiles are manually created.Thesefiles specify jaw and lip positions with key timing information,so that when the graphics and audio are recorded at exactly the same time,the face appears to speak.Key-framing requires a complete face posture for each key position[BL85].Parametric models harness fudical nodes that are moved in synchronization with timings in a script[Par74]. Likewise,anatomical models coordinate muscle contractions to synchronize with the timing information in the scriptfile[Wat87].In essence these techniques are obvious extensions to traditional hand animation.Unfortunately they have the same inherent problem:their lack offlexibility.If the audio is modified,even slightly,ambiguities in the synchronization are created,and the whole process has2.1Lip-reading and the Phonetics of Speech3to be repeated.Automatic lip synchronization can be tackled from two directions:synchro-nization to synthetic speech and synchronization to real speech.In the latter case, the objective is to use speech analysis to automatically identify lip shapes for a given speech sequence.In brief,the speech sequence is transferred into a represen-tation in which the formant frequencies are emphasized and the pitch information largely removed.This representation is then parsed into words.This latter task is difficult,and the former acoustic pre-processing is reasonably effective for driving phonetic scripts[Wei82,LP87].Lip synchronization to a synthesizer is the converse problem,where all the governing speech parameters are known.Hill,Pearce,and Wyvill extended a rule-based synthesizer to incorporate parameters for a3D facial model[HPW88].When the synthesizer scripts were created,facial parameters could also be modified.Once a speech sequence had been generated,it was recorded to the audio channel of a video tape.The facial model was then generated frame-by-frame and recorded in non-real-time to the video section of the tape.Consequently,when the sequence was played back,the face appeared to speak.2.1Lip-reading and the Phonetics of SpeechIn English lip-reading is based on the observation of forty-five phonemes and associated visemes[Wal82].Traditionally,lip-reading has been considered to be a completely visual process developed by the small number of people who are completely deaf.There are,however,three mechanisms employed in visual speech perception:auditory,visual,and audio-visual.Those with hearing impair-ment concern themselves with the audio-visual,placing emphasis on observing the context in which words are spoken,such as posture and facial expression.Speech comprises a mixture of audio frequencies,and every speech sound belongs to one of the two main classes known as vowels and consonants.V owels and consonants belong to basic linguistic units known as phonemes which can be mapped into visible mouth shapes known as visemes.Each vowel has a distin-tive mouth shape,and viseme groups such as p,m,b and f,v can be reliably observed like the vowels,although confusion among individual consonants within each viseme group is more common[McG85].Despite the low threshold between understanding and misunderstanding,the discrete phonetics provide a useful ab-straction,because they group together speech sounds that have common acoustic or articulatory features.We use phonemes and visemes as the basic units of visible articulatory mouth shapes.42BACKGROUND AND PREVIOUS WORK 2.2Speech SynthesisAutomatic text-to-speech synthesis describes the process of generating synthetic speech from arbitrary text input.In most cases this is achieved with a letter-to-sound component and a synthesizer component.For a complete review of automatic speech synthesis,see[Kla87].In general,a letter-to-sound system(LTS)accepts arbitrary ASCII text as input and produces a phonemic transcription as output.The LTS component uses a pronunciation lexicon and possibly a collection of letter-to-sound rules to convert text to phonemes with lexical stress markings.These letter-to-sound rule sets are used to predict the correct pronunciation when a dictionary match is not found.Synthesizers typically accept the input phonemes from the letter-to-sound com-ponent to produce synthetic audible speech.Three classes of segmental synthesis-by-rule techniques have been identified by Klatt[Kla87].They are(1)formant-based rules programs,(2)articulation-based rule programs,and(3)concatenation systems.Each technique attempts to create natural and intelligible synthetic speech from phonemes.A formant synthesizer recreates the speech spectrum using a collection of rules and heuristics to control a digitalfilter model of the vocal tract.Klattalk[Kla87] and DECtalk[BMT83]are examples of formant-based synthesizers.Concatena-tion systems,as the name suggests,synthesize speech by splicing together short segments of parameterized stored speech.For example,speech segments might be stored as sequences of LPC(linear predictive coding)parameters which are used to resynthesize speech.Olive’s Utter system is an example of diphone synthesis in a concatenative system[Oli90].Formant-based rule programs and concatenation systems are both capable of producing intelligible synthetic speech.However,concatenation systems tend to introduce artifacts into the speech as a result of discontinuities at the boundaries between the stored acoustic units.Furthermore,concatenation systems have to store an inventory of sound units that may be on the order of1M bytes per voice for a diphone system.Formant-based synthesizers require far less storage than concatenative systems,but they are restricted in the number and quality of the voices(speakers)produced.In particular,synthesizing a voice with a truly feminine quality has been found to be difficult[Kla87].For the purpose of synchronizing speech to synthetic faces,either the formant-based or concatenative approach could have been used.We use a formant synthesizer for DECface.2.3DECface52.3DECfaceThe previous work described in Section2produces animation in a single frame mode that is subsequently recorded to video tape.In contrast to this previous work this report describes a fundamentally different approach to automatically synchronizing computer generated faces with synthetic speech.The complete process is called DECface which provides a novel form of real-time face-to-face communication.The unique feature of DECface is the ability to generate speech and graphics at real-time rates,where the audio and the graphics are tightly coupled to generate expressive synthetic facial characters.This demands a fundamentally different approach to traditional techniques.Furthermore,to compute synthetic faces and synchronize the audio in real-time requires a powerful computational resource such as an Alpha AXP workstation.3The AlgorithmThis section presents an algorithm to automate the process of synchronizing lip motion to a formant-based speech synthesizer.The key component of the algorithm is the run-time facility that adaptively synchronizes the graphical display of the face to the audio.DECface executes the following sequence of operations:1.Input ASCII text2.Create phonetic transcription from the text3.Generate synthesized speech samples from the text4.Query the audio server and determine the current phoneme from the speechplaybackpute the current mouth shape from nodal trajectories6.Play synthesized speech samples and synchronize the graphical displayStep1,2,and3are an integral component of the algorithm and can be viewed as a pre-processing stage.Therefore both the phonetic transcription and the audio samples can be generated in advance and stored to disk if desired.Steps4,5,and6 are concerned with adaptive synchronization and are repeated until no more speech samples are left to be played.63THE ALGORITHMDisplayFigure1:The synchronization model TTS=text-to-speech.The system requires a Digital-to-Analog Converter(DAC)to produce the ana-log speech waveform.Associated with the audio playback device is a sample clock maintaining a representation of time that can be queried by an application program. We use the term audio server to describe the system software supporting the audio device.Since loss of synchronization is more noticeable in the audio domain,we used the audio server’s clock for the global time base.The audio server’s device clock is sampled during initialization,and thereafter the speech samples are played relative to this initial device time.3.1Text-to-SpeechThe DECface algorithm uses DECtalk.DECtalk is an algorithm that has many implementations;in our case,it is a software implementation.With sufficient computing power available,there is no special hardware or coprocessor needed by DECtalk to synthesize real-time speech.DECtalk comprises three major algo-rithms:(1)the letter-to-sound system,(2)the phonemic synthesizer,and the(3) vocal tract model.The letter-to-sound system accepts arbitrary ASCII text as input and produces a phonemic transcription as output by using a pronunciation lexicon and a collection of letter-to-sound rules for English.As part of this process,the input text may be reformatted to convert numbers and abbreviations into words and punctuation.3.2Time Base7The phonemic synthesizer accepts the phonemic transcription output from the letter-to-sound system and produces parameter control records for the vocal tract model.This component applies intonation,duration,and stress rules to modify the phonemic representation based on phrase-level context.Once these rules have been applied,the phonemic synthesizer calculates parameters for the vocal tract model. The resulting phonetic sequence is also provided to the DECface synchronization component.The vocal tract model accepts the control records from the phonemic synthe-sizer and updates its internal state in order to produce the next frame of synthesized samples.The vocal tract model is a formant synthesizer based on the model de-scribed by Klatt[Kla80].The vocal tract model consists of voiced and unvoiced waveform generators,cascade and parallel resonators,and a summing stage.Fre-quency,bandwidth,and gain parameters in the control record are used to compute thefilter coefficients for the resonators.3.2Time BaseTiming is the most critical component of the DECface algorithm since the au-dio/graphics synchronization can only be achieved when a common global time base exists.The initial time0is recorded from the audio server as thefirst sample of speech is output.By re-sampling the audio server at time1,an exact corre-spondence to a phoneme or phoneme pair can be determined.The relative time since the start of speech output is10and is used to calculate the current viseme mouth shape.The audio server’s sample rate is8000and ideally,the graphical display frame rate is30.Therefore to avoid aliasing artifacts in the interpolation,all values are computed in server time.3.3Mouth DeformationsOnce,the current time relative to0,has been determined from the audio device, the displacement of the mouth nodes can be calculated.Each viseme mouth node is defined with position,where are sequences of nodes defining the geometry and topology of the mouth.To permit a complete mouth shape interpolation,the topology must remainfixed and the nodes in each prototype mouth shape must be in correspondence.An intermediate interpolation position can be calculated between viseme nodes0and1by:0 0111101(1)83THE ALGORITHM where1.The parameter is usually described by a linear or non-linear transformation of where0 1.However,motions based on linear interpolation fail to exhibit the inertial motion characteristics of acceleration and deceleration.A closer interpolated approximation to acceleration and deceleration uses a cosine function to ease in and out of the motion:1cos02(2) The cosine interpolant is an efficient solution and provides acceptable results. However,duringfluent speech,mouth shape rarely converges to discrete viseme targets due to the short interval between positions and the physical properties of the mouth.To emulatefluent speech we need to calculate co-articulated visemes.Piecewise linear interpolation can be used to smooth through node trajectories, and various splines can provide approximations to node positions[Far90].The addition of parameters(such as tension,continuity,and bias)to splines begins to emulate physical motion trajectories[KB82].However,the most significantflaw of splines for animation is that they were originally developed as a means for defining static geometry and not motion trajectories.A dynamic system of nodal displacements,based on the physical behavior of the mouth,can be developed if the initial configuration of the nodes are specified with positions,masses,and velocities;. Once the geometric configuration has been specified,Newtonian physics can be applied:(4)Because we know the state positions0and1,we are in fact calculating the trajectory along the vector01.To resolve the forces that are applied to nodes, it is assumed that0when1.Forces can then be applied to the nodes where is a velocity dependent damping coefficient.It should be noted that forces are not generated from node neighbors as in a mesh,but rather from target node to target node.The mouth shape deformations use a Hookean elastic force based on the separation of the nodes,thereby providing an approximation to elastic behavior of facial tissue:(5)3.4Synchronization9 where the vector separation of the nodes is1.The equations of motion can then be integrated using the projected time t,de-rived from the audio server time,to calculate the new velocities and node positions:104THE FACEMODELv3IScreen Space uv Texture Space(a)(b)Figure2:(a)Polygonal representation of the face.(b)Texture mapping.4.1Face RenderingWhile the wire frame provides a suitable model to observe the motion characteristics of the face,the visual representation is clearly unrealistic.The face polygons can be shaded using Gouraud or Phong shading models,but this often results in faces that look plastic.However,texture mapping is a powerful technique that adds visual detail to synthetic images[Hec86].In particular,texture mapping greatly enhances the realism of synthetic faces.We use an incremental scanline texture mapping technique to achieve realistic faces.Incremental texture mapping is based on the scanline Gouraud shading algo-rithm[Gou71].Instead of interpolating intensity values at polygon vertices,it interpolates texture puting texture coordinates for a poly-gon provides fast mapping into texture space.A color value can then be extracted and applied to the current scanline in screen space(Figure2(b)).For each step along the current scanline in screen space,and are incremented by a constant amount in texture space.Effectively the increment has a slope of and can be used to rapidly index through texture space.For efficiency the coordinate samples texture space and returns a value for the current scanline.11f s ax yu r k581475110977141ix k eh ch n775812811583563 Table1:Phonemes and durations(milliseconds)for the sample sentence.5ExampleThe text“First,a few practice questions.”is used to demonstrate the DECface algorithm in this section.The words are spoken at165words-per-minute in a female voice.These words will be referred to as the sample sentence.The synthesizer produces a sequence of phoneme/duration pairs(Table1)for the sample sentence. The phonemes belong to the two-character arpabet used by DECtalk and have a direct correlation to visemes depicted in Table3.To create Table3,we used shapshots of a persons mouth while uttering CvC and VcV strings(Table4).Figure3is a time displacement graph illustrating three different computed trajectories,while Figure4illustrates every third frame from each of the three trajectories.Trajectory(a)is the cosine displacement that peaks at the viseme mouth shapes.Trajectories(b)and(c)are two dynamic trajectories computed from equation(4).The two trajectories(b)and(c)are controlled by two variables and representing the velocity damping coefficient and the spring constant respectively.The mass remains constant between the two examples at025. For(b)0650and0500and for(c)0150and0850.Figure5 illustrates a sequence of frames of the complete texture mapped facial model speaking the sample sentence.The physical model of facial tissue motion provides acceleration and deceler-ation emulating inertial characteristics as the mouth changes shape.This is most evident during rapid speech,where the mouth does not make complete mouth shapes,but instead produces a blend between shapes under muscular forces.The final result is a more natural looking mouth motion.125EXAMPLE 500100015002000250030003500Time (ms)(a)(b)(c)167.5170172.5175177.5180182.5185Y Disp Figure 3:The mid top lip vertical displacement trajectories of the sample sen-tence.Trajectory (a )is the cosine activity peaking at each phonetic mouth shape.Trajectory (b )is the physical model with a small damping coefficient and a large spring constant,while trajectory (c )is the physical model with a larger damping coefficient and lower spring constant.13Trajectory(a)Trajectory(b)Trajectory(c)58ms466ms773ms1066ms1305ms1635ms1863ms2148ms2722ms Figure4:Every third viseme produced by the sample sentence for each trajectory.147PERFORMANCE 6ImplementationTo implement a real-time version of DECface,we incorporated several existing hardware and software components on an Alpha AXP workstation running the DEC OSF/1Alpha V1.2operating system.The synthesized speech output by the DECtalk vocal tract model can be played on any D/A converter hardware supporting16bit linear PCM(pulse code modula-tion)or8bit-law(log-PCM)encoded samples operating at an8KHz sampling rate.For example,the baseboard CODEC1on Alpha AXP workstations or the DECaudio module[Lev93]may be used.DECaudio is a TURBOchannel I/O peripheral containing highlyflexible audio A/D and D/A converter capabilities.The AudioFile audio server and client library[LPG93]were used to interface to the audio hardware and provide an applications programming interface for the audio portion of DECface.The client library provides the application programming interface for accessing the audio server and audio hardware.DECface uses the X Window System to display the rendered facial images.A Tk[Ous90,Ous91]widget-based user interface facilitates the interactive use of DECface(Figure6).A variety of commands can be piped to DECface.For the speech synthesizer,arbitary text can be created in the text widget and spoken using one of eight internal voices.In addition the user can specify the number of words per minute,as well as the comma and period pause durations.For the face synthesizer,sliders are associated with six linear muscles[Wat87]allowing simple facial expressions to be created.Finally several graphical display characteristics can be selected,including texture or wireframe,a physical simulation,muscles, and an SMP(Software Motion Pictures)clip generator.7PerformanceDECface performance data for the DEC Alpha AXP3000/500workstation(150MHz clock)were collected to provide an overall performance indicator.Both the wire-frame and the texture mapped versions were timed on images of512x320pixels. Table7illustrates frames rates for the wireframe and texture modes.The timing data includes the cost to generate the synthetic speech.The performances of the wireframe and the texture mapped models are both about15frames per second.At these frame rates,convincing lip-synchonization can be created.This is because the texture mapping code had been highly optimized,whereas little effort was spent in improving the performance of the wireframe model.15....silence(f)(rr)(s)(t)ma pause(ax)(f)(yu)Figure5:Frames1through18of the test sequence:“First,a few..”The phonemic characters are indicated below key visemes.To emphasize the mouth articulation the cosine trajectory was computed.167PERFORMANCEFigure6:The InterfaceGeometry Total frames/sec(not computed)Wireframe+Physics15.673053792Texture mapped+Physics15.91 Table2:Performance table for DECface on the DEC3000/500workstation.17 8DiscussionWe expect the naturalness of the synthetic speech,rather than the intelligibility measure,will be the most compelling factor in the presentation of a synthetic character’s voice.One could easily imagine a character with an expressive voice and face closely interacting with the user.While this is easy to envision,significant technical issues have to be addressed both from the graphics and audio domains.For example,how do you create a character capable of telling a joke,or of expressing anger,distress,or other emotional states?Perhaps the most perplexing challenge from a graphics perspective is our critical examination of the face model as it approaches reality.Shortcomings in bland, exaggerated plastic faces are readily dismissed;perhaps our expectations are too low for such caricatures.Whatever the reason,we become much less tolerant when faces of real people are manipulated,since we know precisely what they look like and how they articulate.Therefore,if images of real people are to be used,the articulations have to mimic reality extremely closely.The alternative is to exploit the nuances of caricatures and create a new synthetic character,in much the same way that cartoon animators do today.While lip synchronization is crucial to building articulate personable characters, there are other important characteristics,such as body posture and facial expression, that convey bining facial expression with speech is beyond the scope of this paper and has been deliberately omitted.However,our research has demonstrated that basic expressions can dramatically enhance the realism of the talking face;even simple eye blinks can bring the face to life.It is important to remember that a phonetic transcription is an interpretation im-posed on a continuously varying acoustic signal.Therefore visemes are extensions to phonemes.The DECface algorithm can be extended to co-articulated words,but the viseme synchronization can ultimately be only as good as the text-to-speech synthesizer.While DECface has been designed to operate on2D images,extensions to3D are straightforward.In fact,physically based face systems can be simply modified to incorporate DECface.We plan to incorporate DECface into a3D facial model. 9ConclusionWe have demonstrated an algorithm for automatically synchronizing lip motion to a speech synthesizer in real-time.Theflexibility of the algorithm is derived from the ability to completely synthesize the face and speech explicitly from a stream。

425 BibliographyH.A KAIKE(1974).Markovian representation of stochastic processes and its application to the analysis of autoregressive moving average processes.Annals Institute Statistical Mathematics,vol.26,pp.363-387. B.D.O.A NDERSON and J.B.M OORE(1979).Optimal rmation and System Sciences Series, Prentice Hall,Englewood Cliffs,NJ.T.W.A NDERSON(1971).The Statistical Analysis of Time Series.Series in Probability and Mathematical Statistics,Wiley,New York.R.A NDRE-O BRECHT(1988).A new statistical approach for the automatic segmentation of continuous speech signals.IEEE Trans.Acoustics,Speech,Signal Processing,vol.ASSP-36,no1,pp.29-40.R.A NDRE-O BRECHT(1990).Reconnaissance automatique de parole`a partir de segments acoustiques et de mod`e les de Markov cach´e s.Proc.Journ´e es Etude de la Parole,Montr´e al,May1990(in French).R.A NDRE-O BRECHT and H.Y.S U(1988).Three acoustic labellings for phoneme based continuous speech recognition.Proc.Speech’88,Edinburgh,UK,pp.943-950.U.A PPEL and A.VON B RANDT(1983).Adaptive sequential segmentation of piecewise stationary time rmation Sciences,vol.29,no1,pp.27-56.L.A.A ROIAN and H.L EVENE(1950).The effectiveness of quality control procedures.Jal American Statis-tical Association,vol.45,pp.520-529.K.J.A STR¨OM and B.W ITTENMARK(1984).Computer Controlled Systems:Theory and rma-tion and System Sciences Series,Prentice Hall,Englewood Cliffs,NJ.M.B AGSHAW and R.A.J OHNSON(1975a).The effect of serial correlation on the performance of CUSUM tests-Part II.Technometrics,vol.17,no1,pp.73-80.M.B AGSHAW and R.A.J OHNSON(1975b).The influence of reference values and estimated variance on the ARL of CUSUM tests.Jal Royal Statistical Society,vol.37(B),no3,pp.413-420.M.B AGSHAW and R.A.J OHNSON(1977).Sequential procedures for detecting parameter changes in a time-series model.Jal American Statistical Association,vol.72,no359,pp.593-597.R.K.B ANSAL and P.P APANTONI-K AZAKOS(1986).An algorithm for detecting a change in a stochastic process.IEEE rmation Theory,vol.IT-32,no2,pp.227-235.G.A.B ARNARD(1959).Control charts and stochastic processes.Jal Royal Statistical Society,vol.B.21, pp.239-271.A.E.B ASHARINOV andB.S.F LEISHMAN(1962).Methods of the statistical sequential analysis and their radiotechnical applications.Sovetskoe Radio,Moscow(in Russian).M.B ASSEVILLE(1978).D´e viations par rapport au maximum:formules d’arrˆe t et martingales associ´e es. Compte-rendus du S´e minaire de Probabilit´e s,Universit´e de Rennes I.M.B ASSEVILLE(1981).Edge detection using sequential methods for change in level-Part II:Sequential detection of change in mean.IEEE Trans.Acoustics,Speech,Signal Processing,vol.ASSP-29,no1,pp.32-50.426B IBLIOGRAPHY M.B ASSEVILLE(1982).A survey of statistical failure detection techniques.In Contribution`a la D´e tectionS´e quentielle de Ruptures de Mod`e les Statistiques,Th`e se d’Etat,Universit´e de Rennes I,France(in English). M.B ASSEVILLE(1986).The two-models approach for the on-line detection of changes in AR processes. In Detection of Abrupt Changes in Signals and Dynamical Systems(M.Basseville,A.Benveniste,eds.). Lecture Notes in Control and Information Sciences,LNCIS77,Springer,New York,pp.169-215.M.B ASSEVILLE(1988).Detecting changes in signals and systems-A survey.Automatica,vol.24,pp.309-326.M.B ASSEVILLE(1989).Distance measures for signal processing and pattern recognition.Signal Process-ing,vol.18,pp.349-369.M.B ASSEVILLE and A.B ENVENISTE(1983a).Design and comparative study of some sequential jump detection algorithms for digital signals.IEEE Trans.Acoustics,Speech,Signal Processing,vol.ASSP-31, no3,pp.521-535.M.B ASSEVILLE and A.B ENVENISTE(1983b).Sequential detection of abrupt changes in spectral charac-teristics of digital signals.IEEE rmation Theory,vol.IT-29,no5,pp.709-724.M.B ASSEVILLE and A.B ENVENISTE,eds.(1986).Detection of Abrupt Changes in Signals and Dynamical Systems.Lecture Notes in Control and Information Sciences,LNCIS77,Springer,New York.M.B ASSEVILLE and I.N IKIFOROV(1991).A unified framework for statistical change detection.Proc.30th IEEE Conference on Decision and Control,Brighton,UK.M.B ASSEVILLE,B.E SPIAU and J.G ASNIER(1981).Edge detection using sequential methods for change in level-Part I:A sequential edge detection algorithm.IEEE Trans.Acoustics,Speech,Signal Processing, vol.ASSP-29,no1,pp.24-31.M.B ASSEVILLE, A.B ENVENISTE and G.M OUSTAKIDES(1986).Detection and diagnosis of abrupt changes in modal characteristics of nonstationary digital signals.IEEE rmation Theory,vol.IT-32,no3,pp.412-417.M.B ASSEVILLE,A.B ENVENISTE,G.M OUSTAKIDES and A.R OUG´E E(1987a).Detection and diagnosis of changes in the eigenstructure of nonstationary multivariable systems.Automatica,vol.23,no3,pp.479-489. M.B ASSEVILLE,A.B ENVENISTE,G.M OUSTAKIDES and A.R OUG´E E(1987b).Optimal sensor location for detecting changes in dynamical behavior.IEEE Trans.Automatic Control,vol.AC-32,no12,pp.1067-1075.M.B ASSEVILLE,A.B ENVENISTE,B.G ACH-D EVAUCHELLE,M.G OURSAT,D.B ONNECASE,P.D OREY, M.P REVOSTO and M.O LAGNON(1993).Damage monitoring in vibration mechanics:issues in diagnos-tics and predictive maintenance.Mechanical Systems and Signal Processing,vol.7,no5,pp.401-423.R.V.B EARD(1971).Failure Accommodation in Linear Systems through Self-reorganization.Ph.D.Thesis, Dept.Aeronautics and Astronautics,MIT,Cambridge,MA.A.B ENVENISTE and J.J.F UCHS(1985).Single sample modal identification of a nonstationary stochastic process.IEEE Trans.Automatic Control,vol.AC-30,no1,pp.66-74.A.B ENVENISTE,M.B ASSEVILLE and G.M OUSTAKIDES(1987).The asymptotic local approach to change detection and model validation.IEEE Trans.Automatic Control,vol.AC-32,no7,pp.583-592.A.B ENVENISTE,M.M ETIVIER and P.P RIOURET(1990).Adaptive Algorithms and Stochastic Approxima-tions.Series on Applications of Mathematics,(A.V.Balakrishnan,I.Karatzas,M.Yor,eds.).Springer,New York.A.B ENVENISTE,M.B ASSEVILLE,L.E L G HAOUI,R.N IKOUKHAH and A.S.W ILLSKY(1992).An optimum robust approach to statistical failure detection and identification.IFAC World Conference,Sydney, July1993.B IBLIOGRAPHY427 R.H.B ERK(1973).Some asymptotic aspects of sequential analysis.Annals Statistics,vol.1,no6,pp.1126-1138.R.H.B ERK(1975).Locally most powerful sequential test.Annals Statistics,vol.3,no2,pp.373-381.P.B ILLINGSLEY(1968).Convergence of Probability Measures.Wiley,New York.A.F.B ISSELL(1969).Cusum techniques for quality control.Applied Statistics,vol.18,pp.1-30.M.E.B IVAIKOV(1991).Control of the sample size for recursive estimation of parameters subject to abrupt changes.Automation and Remote Control,no9,pp.96-103.R.E.B LAHUT(1987).Principles and Practice of Information Theory.Addison-Wesley,Reading,MA.I.F.B LAKE and W.C.L INDSEY(1973).Level-crossing problems for random processes.IEEE r-mation Theory,vol.IT-19,no3,pp.295-315.G.B ODENSTEIN and H.M.P RAETORIUS(1977).Feature extraction from the encephalogram by adaptive segmentation.Proc.IEEE,vol.65,pp.642-652.T.B OHLIN(1977).Analysis of EEG signals with changing spectra using a short word Kalman estimator. Mathematical Biosciences,vol.35,pp.221-259.W.B¨OHM and P.H ACKL(1990).Improved bounds for the average run length of control charts based on finite weighted sums.Annals Statistics,vol.18,no4,pp.1895-1899.T.B OJDECKI and J.H OSZA(1984).On a generalized disorder problem.Stochastic Processes and their Applications,vol.18,pp.349-359.L.I.B ORODKIN and V.V.M OTTL’(1976).Algorithm forfinding the jump times of random process equation parameters.Automation and Remote Control,vol.37,no6,Part1,pp.23-32.A.A.B OROVKOV(1984).Theory of Mathematical Statistics-Estimation and Hypotheses Testing,Naouka, Moscow(in Russian).Translated in French under the title Statistique Math´e matique-Estimation et Tests d’Hypoth`e ses,Mir,Paris,1987.G.E.P.B OX and G.M.J ENKINS(1970).Time Series Analysis,Forecasting and Control.Series in Time Series Analysis,Holden-Day,San Francisco.A.VON B RANDT(1983).Detecting and estimating parameters jumps using ladder algorithms and likelihood ratio test.Proc.ICASSP,Boston,MA,pp.1017-1020.A.VON B RANDT(1984).Modellierung von Signalen mit Sprunghaft Ver¨a nderlichem Leistungsspektrum durch Adaptive Segmentierung.Doctor-Engineer Dissertation,M¨u nchen,RFA(in German).S.B RAUN,ed.(1986).Mechanical Signature Analysis-Theory and Applications.Academic Press,London. L.B REIMAN(1968).Probability.Series in Statistics,Addison-Wesley,Reading,MA.G.S.B RITOV and L.A.M IRONOVSKI(1972).Diagnostics of linear systems of automatic regulation.Tekh. Kibernetics,vol.1,pp.76-83.B.E.B RODSKIY and B.S.D ARKHOVSKIY(1992).Nonparametric Methods in Change-point Problems. Kluwer Academic,Boston.L.D.B ROEMELING(1982).Jal Econometrics,vol.19,Special issue on structural change in Econometrics. L.D.B ROEMELING and H.T SURUMI(1987).Econometrics and Structural Change.Dekker,New York. D.B ROOK and D.A.E VANS(1972).An approach to the probability distribution of Cusum run length. Biometrika,vol.59,pp.539-550.J.B RUNET,D.J AUME,M.L ABARR`E RE,A.R AULT and M.V ERG´E(1990).D´e tection et Diagnostic de Pannes.Trait´e des Nouvelles Technologies,S´e rie Diagnostic et Maintenance,Herm`e s,Paris(in French).428B IBLIOGRAPHY S.P.B RUZZONE and M.K AVEH(1984).Information tradeoffs in using the sample autocorrelation function in ARMA parameter estimation.IEEE Trans.Acoustics,Speech,Signal Processing,vol.ASSP-32,no4, pp.701-715.A.K.C AGLAYAN(1980).Necessary and sufficient conditions for detectability of jumps in linear systems. IEEE Trans.Automatic Control,vol.AC-25,no4,pp.833-834.A.K.C AGLAYAN and R.E.L ANCRAFT(1983).Reinitialization issues in fault tolerant systems.Proc.Amer-ican Control Conf.,pp.952-955.A.K.C AGLAYAN,S.M.A LLEN and K.W EHMULLER(1988).Evaluation of a second generation reconfigu-ration strategy for aircraftflight control systems subjected to actuator failure/surface damage.Proc.National Aerospace and Electronic Conference,Dayton,OH.P.E.C AINES(1988).Linear Stochastic Systems.Series in Probability and Mathematical Statistics,Wiley, New York.M.J.C HEN and J.P.N ORTON(1987).Estimation techniques for tracking rapid parameter changes.Intern. Jal Control,vol.45,no4,pp.1387-1398.W.K.C HIU(1974).The economic design of cusum charts for controlling normal mean.Applied Statistics, vol.23,no3,pp.420-433.E.Y.C HOW(1980).A Failure Detection System Design Methodology.Ph.D.Thesis,M.I.T.,L.I.D.S.,Cam-bridge,MA.E.Y.C HOW and A.S.W ILLSKY(1984).Analytical redundancy and the design of robust failure detection systems.IEEE Trans.Automatic Control,vol.AC-29,no3,pp.689-691.Y.S.C HOW,H.R OBBINS and D.S IEGMUND(1971).Great Expectations:The Theory of Optimal Stop-ping.Houghton-Mifflin,Boston.R.N.C LARK,D.C.F OSTH and V.M.W ALTON(1975).Detection of instrument malfunctions in control systems.IEEE Trans.Aerospace Electronic Systems,vol.AES-11,pp.465-473.A.C OHEN(1987).Biomedical Signal Processing-vol.1:Time and Frequency Domain Analysis;vol.2: Compression and Automatic Recognition.CRC Press,Boca Raton,FL.J.C ORGE and F.P UECH(1986).Analyse du rythme cardiaque foetal par des m´e thodes de d´e tection de ruptures.Proc.7th INRIA Int.Conf.Analysis and optimization of Systems.Antibes,FR(in French).D.R.C OX and D.V.H INKLEY(1986).Theoretical Statistics.Chapman and Hall,New York.D.R.C OX and H.D.M ILLER(1965).The Theory of Stochastic Processes.Wiley,New York.S.V.C ROWDER(1987).A simple method for studying run-length distributions of exponentially weighted moving average charts.Technometrics,vol.29,no4,pp.401-407.H.C S¨ORG¨O and L.H ORV´ATH(1988).Nonparametric methods for change point problems.In Handbook of Statistics(P.R.Krishnaiah,C.R.Rao,eds.),vol.7,Elsevier,New York,pp.403-425.R.B.D AVIES(1973).Asymptotic inference in stationary gaussian time series.Advances Applied Probability, vol.5,no3,pp.469-497.J.C.D ECKERT,M.N.D ESAI,J.J.D EYST and A.S.W ILLSKY(1977).F-8DFBW sensor failure identification using analytical redundancy.IEEE Trans.Automatic Control,vol.AC-22,no5,pp.795-803.M.H.D E G ROOT(1970).Optimal Statistical Decisions.Series in Probability and Statistics,McGraw-Hill, New York.J.D ESHAYES and D.P ICARD(1979).Tests de ruptures dans un mod`e pte-Rendus de l’Acad´e mie des Sciences,vol.288,Ser.A,pp.563-566(in French).B IBLIOGRAPHY429 J.D ESHAYES and D.P ICARD(1983).Ruptures de Mod`e les en Statistique.Th`e ses d’Etat,Universit´e deParis-Sud,Orsay,France(in French).J.D ESHAYES and D.P ICARD(1986).Off-line statistical analysis of change-point models using non para-metric and likelihood methods.In Detection of Abrupt Changes in Signals and Dynamical Systems(M. Basseville,A.Benveniste,eds.).Lecture Notes in Control and Information Sciences,LNCIS77,Springer, New York,pp.103-168.B.D EVAUCHELLE-G ACH(1991).Diagnostic M´e canique des Fatigues sur les Structures Soumises`a des Vibrations en Ambiance de Travail.Th`e se de l’Universit´e Paris IX Dauphine(in French).B.D EVAUCHELLE-G ACH,M.B ASSEVILLE and A.B ENVENISTE(1991).Diagnosing mechanical changes in vibrating systems.Proc.SAFEPROCESS’91,Baden-Baden,FRG,pp.85-89.R.D I F RANCESCO(1990).Real-time speech segmentation using pitch and convexity jump models:applica-tion to variable rate speech coding.IEEE Trans.Acoustics,Speech,Signal Processing,vol.ASSP-38,no5, pp.741-748.X.D ING and P.M.F RANK(1990).Fault detection via factorization approach.Systems and Control Letters, vol.14,pp.431-436.J.L.D OOB(1953).Stochastic Processes.Wiley,New York.V.D RAGALIN(1988).Asymptotic solutions in detecting a change in distribution under an unknown param-eter.Statistical Problems of Control,Issue83,Vilnius,pp.45-52.B.D UBUISSON(1990).Diagnostic et Reconnaissance des Formes.Trait´e des Nouvelles Technologies,S´e rie Diagnostic et Maintenance,Herm`e s,Paris(in French).A.J.D UNCAN(1986).Quality Control and Industrial Statistics,5th edition.Richard D.Irwin,Inc.,Home-wood,IL.J.D URBIN(1971).Boundary-crossing probabilities for the Brownian motion and Poisson processes and techniques for computing the power of the Kolmogorov-Smirnov test.Jal Applied Probability,vol.8,pp.431-453.J.D URBIN(1985).Thefirst passage density of the crossing of a continuous Gaussian process to a general boundary.Jal Applied Probability,vol.22,no1,pp.99-122.A.E MAMI-N AEINI,M.M.A KHTER and S.M.R OCK(1988).Effect of model uncertainty on failure detec-tion:the threshold selector.IEEE Trans.Automatic Control,vol.AC-33,no12,pp.1106-1115.J.D.E SARY,F.P ROSCHAN and D.W.W ALKUP(1967).Association of random variables with applications. Annals Mathematical Statistics,vol.38,pp.1466-1474.W.D.E WAN and K.W.K EMP(1960).Sampling inspection of continuous processes with no autocorrelation between successive results.Biometrika,vol.47,pp.263-280.G.F AVIER and A.S MOLDERS(1984).Adaptive smoother-predictors for tracking maneuvering targets.Proc. 23rd Conf.Decision and Control,Las Vegas,NV,pp.831-836.W.F ELLER(1966).An Introduction to Probability Theory and Its Applications,vol.2.Series in Probability and Mathematical Statistics,Wiley,New York.R.A.F ISHER(1925).Theory of statistical estimation.Proc.Cambridge Philosophical Society,vol.22, pp.700-725.M.F ISHMAN(1988).Optimization of the algorithm for the detection of a disorder,based on the statistic of exponential smoothing.In Statistical Problems of Control,Issue83,Vilnius,pp.146-151.R.F LETCHER(1980).Practical Methods of Optimization,2volumes.Wiley,New York.P.M.F RANK(1990).Fault diagnosis in dynamic systems using analytical and knowledge based redundancy -A survey and new results.Automatica,vol.26,pp.459-474.430B IBLIOGRAPHY P.M.F RANK(1991).Enhancement of robustness in observer-based fault detection.Proc.SAFEPRO-CESS’91,Baden-Baden,FRG,pp.275-287.P.M.F RANK and J.W¨UNNENBERG(1989).Robust fault diagnosis using unknown input observer schemes. In Fault Diagnosis in Dynamic Systems-Theory and Application(R.Patton,P.Frank,R.Clark,eds.). International Series in Systems and Control Engineering,Prentice Hall International,London,UK,pp.47-98.K.F UKUNAGA(1990).Introduction to Statistical Pattern Recognition,2d ed.Academic Press,New York. S.I.G ASS(1958).Linear Programming:Methods and Applications.McGraw Hill,New York.W.G E and C.Z.F ANG(1989).Extended robust observation approach for failure isolation.Int.Jal Control, vol.49,no5,pp.1537-1553.W.G ERSCH(1986).Two applications of parametric time series modeling methods.In Mechanical Signature Analysis-Theory and Applications(S.Braun,ed.),chap.10.Academic Press,London.J.J.G ERTLER(1988).Survey of model-based failure detection and isolation in complex plants.IEEE Control Systems Magazine,vol.8,no6,pp.3-11.J.J.G ERTLER(1991).Analytical redundancy methods in fault detection and isolation.Proc.SAFEPRO-CESS’91,Baden-Baden,FRG,pp.9-22.B.K.G HOSH(1970).Sequential Tests of Statistical Hypotheses.Addison-Wesley,Cambridge,MA.I.N.G IBRA(1975).Recent developments in control charts techniques.Jal Quality Technology,vol.7, pp.183-192.J.P.G ILMORE and R.A.M C K ERN(1972).A redundant strapdown inertial reference unit(SIRU).Jal Space-craft,vol.9,pp.39-47.M.A.G IRSHICK and H.R UBIN(1952).A Bayes approach to a quality control model.Annals Mathematical Statistics,vol.23,pp.114-125.A.L.G OEL and S.M.W U(1971).Determination of the ARL and a contour nomogram for CUSUM charts to control normal mean.Technometrics,vol.13,no2,pp.221-230.P.L.G OLDSMITH and H.W HITFIELD(1961).Average run lengths in cumulative chart quality control schemes.Technometrics,vol.3,pp.11-20.G.C.G OODWIN and K.S.S IN(1984).Adaptive Filtering,Prediction and rmation and System Sciences Series,Prentice Hall,Englewood Cliffs,NJ.R.M.G RAY and L.D.D AVISSON(1986).Random Processes:a Mathematical Approach for Engineers. Information and System Sciences Series,Prentice Hall,Englewood Cliffs,NJ.C.G UEGUEN and L.L.S CHARF(1980).Exact maximum likelihood identification for ARMA models:a signal processing perspective.Proc.1st EUSIPCO,Lausanne.D.E.G USTAFSON, A.S.W ILLSKY,J.Y.W ANG,M.C.L ANCASTER and J.H.T RIEBWASSER(1978). ECG/VCG rhythm diagnosis using statistical signal analysis.Part I:Identification of persistent rhythms. Part II:Identification of transient rhythms.IEEE Trans.Biomedical Engineering,vol.BME-25,pp.344-353 and353-361.F.G USTAFSSON(1991).Optimal segmentation of linear regression parameters.Proc.IFAC/IFORS Symp. Identification and System Parameter Estimation,Budapest,pp.225-229.T.H¨AGGLUND(1983).New Estimation Techniques for Adaptive Control.Ph.D.Thesis,Lund Institute of Technology,Lund,Sweden.T.H¨AGGLUND(1984).Adaptive control of systems subject to large parameter changes.Proc.IFAC9th World Congress,Budapest.B IBLIOGRAPHY431 P.H ALL and C.C.H EYDE(1980).Martingale Limit Theory and its Application.Probability and Mathemat-ical Statistics,a Series of Monographs and Textbooks,Academic Press,New York.W.J.H ALL,R.A.W IJSMAN and J.K.G HOSH(1965).The relationship between sufficiency and invariance with applications in sequential analysis.Ann.Math.Statist.,vol.36,pp.576-614.E.J.H ANNAN and M.D EISTLER(1988).The Statistical Theory of Linear Systems.Series in Probability and Mathematical Statistics,Wiley,New York.J.D.H EALY(1987).A note on multivariate CuSum procedures.Technometrics,vol.29,pp.402-412.D.M.H IMMELBLAU(1970).Process Analysis by Statistical Methods.Wiley,New York.D.M.H IMMELBLAU(1978).Fault Detection and Diagnosis in Chemical and Petrochemical Processes. Chemical Engineering Monographs,vol.8,Elsevier,Amsterdam.W.G.S.H INES(1976a).A simple monitor of a system with sudden parameter changes.IEEE r-mation Theory,vol.IT-22,no2,pp.210-216.W.G.S.H INES(1976b).Improving a simple monitor of a system with sudden parameter changes.IEEE rmation Theory,vol.IT-22,no4,pp.496-499.D.V.H INKLEY(1969).Inference about the intersection in two-phase regression.Biometrika,vol.56,no3, pp.495-504.D.V.H INKLEY(1970).Inference about the change point in a sequence of random variables.Biometrika, vol.57,no1,pp.1-17.D.V.H INKLEY(1971).Inference about the change point from cumulative sum-tests.Biometrika,vol.58, no3,pp.509-523.D.V.H INKLEY(1971).Inference in two-phase regression.Jal American Statistical Association,vol.66, no336,pp.736-743.J.R.H UDDLE(1983).Inertial navigation system error-model considerations in Kalmanfiltering applica-tions.In Control and Dynamic Systems(C.T.Leondes,ed.),Academic Press,New York,pp.293-339.J.S.H UNTER(1986).The exponentially weighted moving average.Jal Quality Technology,vol.18,pp.203-210.I.A.I BRAGIMOV and R.Z.K HASMINSKII(1981).Statistical Estimation-Asymptotic Theory.Applications of Mathematics Series,vol.16.Springer,New York.R.I SERMANN(1984).Process fault detection based on modeling and estimation methods-A survey.Auto-matica,vol.20,pp.387-404.N.I SHII,A.I WATA and N.S UZUMURA(1979).Segmentation of nonstationary time series.Int.Jal Systems Sciences,vol.10,pp.883-894.J.E.J ACKSON and R.A.B RADLEY(1961).Sequential and tests.Annals Mathematical Statistics, vol.32,pp.1063-1077.B.J AMES,K.L.J AMES and D.S IEGMUND(1988).Conditional boundary crossing probabilities with appli-cations to change-point problems.Annals Probability,vol.16,pp.825-839.M.K.J EERAGE(1990).Reliability analysis of fault-tolerant IMU architectures with redundant inertial sen-sors.IEEE Trans.Aerospace and Electronic Systems,vol.AES-5,no.7,pp.23-27.N.L.J OHNSON(1961).A simple theoretical approach to cumulative sum control charts.Jal American Sta-tistical Association,vol.56,pp.835-840.N.L.J OHNSON and F.C.L EONE(1962).Cumulative sum control charts:mathematical principles applied to their construction and use.Parts I,II,III.Industrial Quality Control,vol.18,pp.15-21;vol.19,pp.29-36; vol.20,pp.22-28.432B IBLIOGRAPHY R.A.J OHNSON and M.B AGSHAW(1974).The effect of serial correlation on the performance of CUSUM tests-Part I.Technometrics,vol.16,no.1,pp.103-112.H.L.J ONES(1973).Failure Detection in Linear Systems.Ph.D.Thesis,Dept.Aeronautics and Astronautics, MIT,Cambridge,MA.R.H.J ONES,D.H.C ROWELL and L.E.K APUNIAI(1970).Change detection model for serially correlated multivariate data.Biometrics,vol.26,no2,pp.269-280.M.J URGUTIS(1984).Comparison of the statistical properties of the estimates of the change times in an autoregressive process.In Statistical Problems of Control,Issue65,Vilnius,pp.234-243(in Russian).T.K AILATH(1980).Linear rmation and System Sciences Series,Prentice Hall,Englewood Cliffs,NJ.L.V.K ANTOROVICH and V.I.K RILOV(1958).Approximate Methods of Higher Analysis.Interscience,New York.S.K ARLIN and H.M.T AYLOR(1975).A First Course in Stochastic Processes,2d ed.Academic Press,New York.S.K ARLIN and H.M.T AYLOR(1981).A Second Course in Stochastic Processes.Academic Press,New York.D.K AZAKOS and P.P APANTONI-K AZAKOS(1980).Spectral distance measures between gaussian pro-cesses.IEEE Trans.Automatic Control,vol.AC-25,no5,pp.950-959.K.W.K EMP(1958).Formula for calculating the operating characteristic and average sample number of some sequential tests.Jal Royal Statistical Society,vol.B-20,no2,pp.379-386.K.W.K EMP(1961).The average run length of the cumulative sum chart when a V-mask is used.Jal Royal Statistical Society,vol.B-23,pp.149-153.K.W.K EMP(1967a).Formal expressions which can be used for the determination of operating character-istics and average sample number of a simple sequential test.Jal Royal Statistical Society,vol.B-29,no2, pp.248-262.K.W.K EMP(1967b).A simple procedure for determining upper and lower limits for the average sample run length of a cumulative sum scheme.Jal Royal Statistical Society,vol.B-29,no2,pp.263-265.D.P.K ENNEDY(1976).Some martingales related to cumulative sum tests and single server queues.Stochas-tic Processes and Appl.,vol.4,pp.261-269.T.H.K ERR(1980).Statistical analysis of two-ellipsoid overlap test for real time failure detection.IEEE Trans.Automatic Control,vol.AC-25,no4,pp.762-772.T.H.K ERR(1982).False alarm and correct detection probabilities over a time interval for restricted classes of failure detection algorithms.IEEE rmation Theory,vol.IT-24,pp.619-631.T.H.K ERR(1987).Decentralizedfiltering and redundancy management for multisensor navigation.IEEE Trans.Aerospace and Electronic systems,vol.AES-23,pp.83-119.Minor corrections on p.412and p.599 (May and July issues,respectively).R.A.K HAN(1978).Wald’s approximations to the average run length in cusum procedures.Jal Statistical Planning and Inference,vol.2,no1,pp.63-77.R.A.K HAN(1979).Somefirst passage problems related to cusum procedures.Stochastic Processes and Applications,vol.9,no2,pp.207-215.R.A.K HAN(1981).A note on Page’s two-sided cumulative sum procedures.Biometrika,vol.68,no3, pp.717-719.B IBLIOGRAPHY433 V.K IREICHIKOV,V.M ANGUSHEV and I.N IKIFOROV(1990).Investigation and application of CUSUM algorithms to monitoring of sensors.In Statistical Problems of Control,Issue89,Vilnius,pp.124-130(in Russian).G.K ITAGAWA and W.G ERSCH(1985).A smoothness prior time-varying AR coefficient modeling of non-stationary covariance time series.IEEE Trans.Automatic Control,vol.AC-30,no1,pp.48-56.N.K LIGIENE(1980).Probabilities of deviations of the change point estimate in statistical models.In Sta-tistical Problems of Control,Issue83,Vilnius,pp.80-86(in Russian).N.K LIGIENE and L.T ELKSNYS(1983).Methods of detecting instants of change of random process prop-erties.Automation and Remote Control,vol.44,no10,Part II,pp.1241-1283.J.K ORN,S.W.G ULLY and A.S.W ILLSKY(1982).Application of the generalized likelihood ratio algorithm to maneuver detection and estimation.Proc.American Control Conf.,Arlington,V A,pp.792-798.P.R.K RISHNAIAH and B.Q.M IAO(1988).Review about estimation of change points.In Handbook of Statistics(P.R.Krishnaiah,C.R.Rao,eds.),vol.7,Elsevier,New York,pp.375-402.P.K UDVA,N.V ISWANADHAM and A.R AMAKRISHNAN(1980).Observers for linear systems with unknown inputs.IEEE Trans.Automatic Control,vol.AC-25,no1,pp.113-115.S.K ULLBACK(1959).Information Theory and Statistics.Wiley,New York(also Dover,New York,1968). K.K UMAMARU,S.S AGARA and T.S¨ODERSTR¨OM(1989).Some statistical methods for fault diagnosis for dynamical systems.In Fault Diagnosis in Dynamic Systems-Theory and Application(R.Patton,P.Frank,R. Clark,eds.).International Series in Systems and Control Engineering,Prentice Hall International,London, UK,pp.439-476.A.K USHNIR,I.N IKIFOROV and I.S AVIN(1983).Statistical adaptive algorithms for automatic detection of seismic signals-Part I:One-dimensional case.In Earthquake Prediction and the Study of the Earth Structure,Naouka,Moscow(Computational Seismology,vol.15),pp.154-159(in Russian).L.L ADELLI(1990).Diffusion approximation for a pseudo-likelihood test process with application to de-tection of change in stochastic system.Stochastics and Stochastics Reports,vol.32,pp.1-25.T.L.L A¨I(1974).Control charts based on weighted sums.Annals Statistics,vol.2,no1,pp.134-147.T.L.L A¨I(1981).Asymptotic optimality of invariant sequential probability ratio tests.Annals Statistics, vol.9,no2,pp.318-333.D.G.L AINIOTIS(1971).Joint detection,estimation,and system identifirmation and Control, vol.19,pp.75-92.M.R.L EADBETTER,G.L INDGREN and H.R OOTZEN(1983).Extremes and Related Properties of Random Sequences and Processes.Series in Statistics,Springer,New York.L.L E C AM(1960).Locally asymptotically normal families of distributions.Univ.California Publications in Statistics,vol.3,pp.37-98.L.L E C AM(1986).Asymptotic Methods in Statistical Decision Theory.Series in Statistics,Springer,New York.E.L.L EHMANN(1986).Testing Statistical Hypotheses,2d ed.Wiley,New York.J.P.L EHOCZKY(1977).Formulas for stopped diffusion processes with stopping times based on the maxi-mum.Annals Probability,vol.5,no4,pp.601-607.H.R.L ERCHE(1980).Boundary Crossing of Brownian Motion.Lecture Notes in Statistics,vol.40,Springer, New York.L.L JUNG(1987).System Identification-Theory for the rmation and System Sciences Series, Prentice Hall,Englewood Cliffs,NJ.。

第38卷第1期2024年1月山东理工大学学报(自然科学版)Journal of Shandong University of Technology(Natural Science Edition)Vol.38No.1Jan.2024收稿日期:20221209基金项目:陕西省自然科学基金项目(2018JQ1043)第一作者:栗雪娟,女,lxj_zk@;通信作者:王丹,女,1611182118@文章编号:1672-6197(2024)01-0073-06Cahn-Hilliard 方程的一个超紧致有限差分格式栗雪娟,王丹(西安建筑科技大学理学院,陕西西安710055)摘要:研究四阶Cahn-Hilliard 方程的数值求解方法㊂给出组合型超紧致差分格式,将其用于四阶Cahn-Hilliard 方程的空间导数离散,采用四阶Runge-Kutta 格式离散时间导数,将二者结合得到四阶Cahn-Hilliard 方程的离散格式,并给出了该格式的误差估计㊂通过编程计算得到其数值解,并与精确解进行对比,结果表明本文的数值方法误差小,验证了所提方法的有效性和可行性㊂关键词:四阶Cahn-Hilliard 方程;组合型超紧致差分方法;四阶Runge-Kutta 方法;误差估计中图分类号:TB532.1;TB553文献标志码:AA supercompact finite difference scheme for Cahn-Hilliard equationsLI Xuejuan,WANG Dan(School of Science,Xiᶄan University of Architecture and Technology,Xiᶄan 710055,China)Abstract :A numerical method for solving the fourth order Cahn-Hilliard equation is studied.The combi-national ultra-compact difference scheme is given and applied to the spatial derivative discretization of the fourth order Cahn-Hilliard equation.The fourth-order Runge-Kutta scheme is used to discrete time deriv-atives.The discrete scheme of the fourth order Cahn-Hilliard equation is obtained by combining the two methods,and the error estimate of the scheme is given.Finally,the numerical solution is obtained by programming and compared with the exact solution.The results show that the numerical method in this paper has a small error,verifying the effectiveness and feasibility of the proposed method.Keywords :fourth order Cahn-Hilliard equation;combinational supercompact difference scheme;fourthorder Runge-Kutta;error estimation㊀㊀本文考虑的四阶Cahn-Hilliard 方程为u t -f u ()xx +ku xxxx =0,x ɪ0,2π[],t >0,u x ,0()=u 0x (),x ɪ0,2π[],u 0,t ()=0,u 2π,t ()=0,t >0,ìîíïïïï(1)式中:求解区域为0,2π[],且kn ȡ0;f u ()为光滑函数;u 0x ()表示t =0时刻的初值;u t 表示u 关于时间t 求偏导数,u t =∂u∂t;f u ()xx表示f u ()关于x求二阶偏导数,f u ()xx=∂2f u ()∂x 2;u xxxx 表示u 关于x 求四阶偏导数,u xxxx=∂4u∂x4;u 是混合物中某种物质的浓度,被称为相变量㊂1958年,Cahn 和Hilliard 提出Cahn-Hilliard 方程,该方程最早被用来描述在温度降低时两种均匀的混合物所发生的相分离现象㊂随着学者对该方程的研究越来越深入,该方程的应用也越来越广泛,特别是在材料科学和物理学等领域中有广泛的应用[1-3]㊂㊀Cahn-Hilliard 方程的数值解法目前已有很多研究,文献[4]使用了全离散有限元方法,文献[5]使用了一类二阶稳定的Crank-Nicolson /Adams-Bashforth 离散化的一致性有限元逼近方法,文献[6-7]使用了有限元方法,文献[8]使用了不连续伽辽金有限元方法,文献[9]使用了Cahn-Hilliard 方程的完全离散谱格式,文献[10]使用了高阶超紧致有限差分方法,文献[11]使用了高阶优化组合型紧致有限差分方法㊂综上所述,本文拟对Cahn-Hilliard 方程构造一种新的超紧致差分格式,将空间组合型超紧致差分方法和修正的时间四阶Runge-Kutta 方法相结合,求解Cahn-Hilliard 方程的数值解,得到相对于现有广义格式精度更高的数值求解格式,并对组合型超紧致差分格式进行误差估计,最后通过数值算例验证该方法的可行性㊂1㊀高阶精度数值求解方法1.1㊀空间组合型超紧致差分格式早期的紧致差分格式是在Hermite 多项式的基础上构造而来的,Hermite 多项式中连续三个节点的一阶导数㊁二阶导数和函数值的数值关系可以表示为ð1k =-1a k f i +k +b k fᶄi +k +c k fᵡi +k ()=0㊂(2)1998年,Krishnan 提出如下紧致差分格式:a 1fᶄi -1+a 0fᶄi +a 2fᶄi +1+hb 1fᵡi -1+b 0fᵡi +b 2fᵡi +1()=1h c 1f i -2+c 2f i -1+c 0f i +c 3f i +1+c 4f i +2(),(3)式中:h 为空间网格间距;a 1,a 0,a 2,b 1,b 0,b 2,c 1,c 2,c 0,c 3,c 4均表示差分格式系数;f i 表示i 节点的函数值;fᶄi 和fᵡi 分别表示i 节点的一阶导数值和二阶导数值;f i -1,f i -2,f i +1,f i +2分别表示i 节点依次向前两个节点和依次向后两个节点的函数值;fᶄi -1,fᶄi +1分别表示i 节点依次向前一个节点和依次向后一个节点的一阶导数值;fᵡi -1,fᵡi +1分别表示i 节点依次向前一个节点和依次向后一个节点的二阶导数值㊂式(2)对应f (x )展开以x i 为邻域的泰勒级数为f x ()=f x i ()+hfᶄx i ()+h 2fᵡx i ()2!+㊀㊀㊀㊀㊀h3f‴x i ()3!+h 4f 4()x i ()4!+h 5f 5()x i ()5!+h 6f 6()x i ()6!+h 7f 7()x i ()7!㊂㊀㊀(4)㊀㊀差分格式的各项系数由式(3)决定,可得到如下的三点六阶超紧致差分格式:716fᶄi +1+fᶄi -1()+fᶄi -h 16fᵡi +1-fᵡi -1()=㊀㊀1516h f i +1-f i -1(),98h fᶄi +1-fᶄi -1()+fᵡi -18fᵡi +1+fᵡi -1()=㊀㊀3h 2f i +1-2f i +f i -1()ìîíïïïïïïïïïï(5)为优化三点六阶紧致差分格式,并保持较好的数值频散,将迎风机制[12]引入式(5),构造出如下三点五阶迎风型超紧致差分格式:78fᶄi -1+fᶄi +h 19fᵡi -1-718fᵡi -172fᵡi +1()=㊀㊀1h -10148f i -1+73f i -1148f i +1(),25fᵡi -1+fᵡi +1h 1910fᶄi -1+165fᶄi +910fᶄi +1()=㊀㊀1h 2-135f i -1-45f i +175f i +1()㊂ìîíïïïïïïïïïï(6)左右边界可达到三阶精度紧致格式:fᶄ1-132fᶄ2+fᶄ3()+3h4fᵡ2+fᵡ3()=㊀㊀-12h f 3-f 2(),fᵡ1+3728h fᶄ3-fᶄ2()+3914h fᶄ1-3356fᵡ3-fᵡ2()=㊀㊀f 3-2f 1+f 2(),ìîíïïïïïïïï(7)fᶄN -132fᶄN -2+fᶄN -1()-3h 4fᵡN -2+fᵡN -1()=㊀㊀12h f N -2-f N -1(),fᵡN -3728h (fᶄN -2-fᶄN -1)-3914h fᶄN -3356(fᵡN -2-㊀㊀fᵡN -1)=1314h 2f N -2-2f N +f N -1()㊂ìîíïïïïïïïïïï(8)上述组合型超紧致差分格式只需要相邻的三个节点便可以同时求得一阶导数和二阶导数的五阶精度近似值,比普通差分格式的节点更少,降低了计算量㊂为便于编程计算,将上述构造的组合型超紧致差分格式重写为矩阵表达形式㊂假设U 为位移矩阵,其大小为m ˑn ,则求一阶导数和二阶导数的离47山东理工大学学报(自然科学版)2024年㊀散过程可以用矩阵运算表示为AF=BU,(9)结合内点的三点五阶迎风型超紧致差分格式和边界点的三点三阶差分格式,组成式(9)中等式左边的矩阵A和等式右边的矩阵B,大小分别为2mˑ2n 和2mˑn;F为奇数行为空间一阶导数和偶数行为空间二阶导数组成的矩阵,大小为2mˑn㊂以上矩阵分别为:A=10-13/23h/4-13/23h/439/14h1-37/28h33/5637/28h-33/567/8h/91-7h/180-h/7219/10h2/516/5h19/1007/8h/91-7h/180-h/7219/10h2/516/5h19/100⋱⋱⋱⋱⋱⋱7/8h/91-7h/180-h/7219/10h2/516/5h19/100-13/2-3h/4-13/2-3h/410-37/28h-33/5637/28h33/56-39/14h1éëêêêêêêêêêêêêêêêêùûúúúúúúúúúúúúúúúú,(10)F=∂u∂x()1,1∂u∂x()1,2∂u∂x()1,n-1∂u∂x()1,n∂2u∂x2()1,1∂2u∂x2()1,2 ∂2u∂x2()1,n-1∂2u∂x2()1,n︙︙︙︙∂u∂x()m,1∂u∂x()m,2∂u∂x()m,n-1∂u∂x()m,n∂2u∂x2()m,1∂2u∂x2()m,2 ∂2u∂x2()m,n-1∂2u∂x2()m,néëêêêêêêêêêêêêêùûúúúúúúúúúúúúú,(11) B=012/h-12/h-13/7h213/14h213/14h2-101/48h7/3h-11/48h-13/5h2-4/5h217/5h2-101/48h27/3h-11/48h-13/5h2-4/5h217/5h2⋱⋱⋱-101/48h7/3h-11/48h-13/5h2-4/5h217/5h2012/h-12/h-13/7h213/14h213/14h2éëêêêêêêêêêêêêêêêêùûúúúúúúúúúúúúúúúú,(12)U=u1,1u1,2 u1,n-1u1,nu2,1u2,2 u2,n-1u2,n︙︙︙︙u m-1,1u m-1,2 u m-1,n-1u m-1,nu m,1u m,2 u m,n-1u m,néëêêêêêêêùûúúúúúúú㊂(13)㊀㊀由式(9)可得F=A-1BU㊂(14)㊀㊀解线性代数方程组(9)可得Cahn-Hilliard方程的空间一阶导数和二阶导数㊂对于四阶导数,可将已求得的二阶导数替代式(14)中的U,再次使用式(14)进行求取㊂57第1期㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀栗雪娟,等:Cahn-Hilliard方程的一个超紧致有限差分格式1.2㊀时间离散格式在对很多偏微分方程的数值求解中不仅需要高精度的空间离散格式,同时还需要高精度的时间离散格式㊂普通的一阶精度时间离散格式显然满足不了高精度计算要求,因此本文选用时间四阶Runge-Kutta 格式进行时间离散㊂Runge-Kutta 方法是基于欧拉方法改进后的求解偏微分方程的常用方法,这种方法不仅计算效率高,而且稳定性好㊂格式的推算过程如下:假设求解方程为∂u∂t+F u ()=0,(15)式中F 是对空间变量的微分算子,则修正的四阶Runge-Kutta 格式为u 0i =u n i ,u 1i =u n i-Δt 4F u ()()0i,u 2i =u ni -Δt 3F u ()()1i,u 3i =u n i-Δt 2F u ()()2i,u n +1i =u n i -Δt F u ()()3i ㊂ìîíïïïïïïïïïïïï(16)1.3㊀误差估计以五阶精度将fᶄi -1,fᶄi +1,fᵡi -1,fᵡi +1泰勒级数展开:fᶄi -1=fᶄi -hfᵡi +h 22!f (3)i -h 33!f (4)i +㊀㊀h 44!f (5)i -h 55!f (6)i ,fᶄi +1=fᶄi +hfᵡi +h 22!f (3)i +h 33!f (4)i+㊀㊀h 44!f (5)i +h 55!f (6)i ,fᵡi -1=fᵡi -hf (3)i +h 22!f (4)i -h 33!f (5)i+㊀㊀h 44!f (6)i -h 55!f (7)i ,fᵡi +1=fᵡi +hf (3)i +h 22!f (4)i +h 33!f (5)i +㊀㊀h 44!f (6)i +h 55!f (7)i ㊂ìîíïïïïïïïïïïïïïïïïïïïïïïïï(17)将式(17)代入式(6),所求得组合型超紧致差分格式的一阶导数及二阶导数对应的截断误差为:78fᶄi -1+fᶄi +h19fᵡi -1-718fᵡi -172fᵡi +1()=㊀1h -10148f i -1+73f i -1148f i +1()+78640f 6()ih 5,25fᵡi -1+fᵡi +1h 1910fᶄi -1+165fᶄi +910fᶄi +1()=㊀-135f i -1-45f i +175f i +1()-5125200f 7()i h 5,ìîíïïïïïïïïïï(18)78640f 6()i h 5ʈ8.101ˑ10-4f 6()i h 5,5125200f 7()ih 5ʈ2.023ˑ10-3f 7()i h 5㊂ìîíïïïï(19)㊀㊀使用组合型超紧致差分格式的好处是在每一个网格点上存在一个一阶和二阶连续导数的多项式㊂本文比较了组合型超紧致差分格式和现有广义格式的一阶导数和二阶导数的截断误差:fᶄi +αfᶄi +1+fᶄi -1()+βfᶄi +2+fᶄi -2()=㊀㊀a f i +1-f i -12h +b f i +2-f i -24h +c f i +3-f i -36h ,fᵡi +αfᵡi +1+fᵡi -1()+βfᵡi +2+fᵡi -2()=㊀㊀a f i +1-2f i +f i -1h 2+b f i +2-2f i +f i -24h2+㊀㊀c f i +3-2f i +f i -39h 2,ìîíïïïïïïïïïïï(20)式中参数α,β,a ,b ,c 在各种格式中取不同的值(表1,表2)㊂本文发现在各种方案中,组合型超紧致差分格式的截断误差最小㊂表1㊀不同格式一阶导数的截断误差格式αβa b c 截断误差二阶中心010013!f 3()ih 2标准Padeᶄ格式1/403/20-15f 5()ih 4六阶中心03/2-3/51/1036ˑ17!f 7()ih 6五阶迎风143ˑ16!f 6()ih 5表2㊀不同格式二阶导数的截断误差格式αβa b c 截断误差二阶中心01002ˑ14!f 4()ih 2标准Padeᶄ格式1/1006/50185ˑ16!f 6()ih 4六阶中心03/2-3/51/1072ˑ18!f 8()ih 6五阶迎风165ˑ17!f 7()ih 567山东理工大学学报(自然科学版)2024年㊀2㊀数值算例误差范数L 1和L 2的定义为:L 1=1N ðNi =1u -U ,L 2=1N ðNi =1u -U ()2㊂对四阶Cahn-Hilliard 取f u ()=u 2,k =2,在边界条件u 0,t ()=u 2π,t ()=0下的计算区域为0,2π[],方程的精确解为u x ,t ()=e -tsin x2,数值解为U ㊂对给出的数值算例,计算误差范数L 1和L 2,并采用四种方法进行数值模拟,对其数值结果进行误差分析和对比,结果见表3,本文所使用方法效果最佳,由此证明所提方法的有效性和可行性㊂表3㊀0.5s 时刻精确度测试结果(N =10)方法L 1误差L 2误差间断有限元格式1.56235ˑ10-21.37823ˑ10-2普通中心差分格式1.66667ˑ10-18.33333ˑ10-2紧致差分格式7.14286ˑ10-31.78571ˑ10-3组合型超紧致差分格式6.48148ˑ10-36.34921ˑ10-4㊀㊀用本文提出的式(6) 式(8)和式(16)计算算例,图1 图3给出了不同时刻数值解与精确解的(a)精确解(b)数值解图1㊀0.1s 的精确解与数值解(a)精确解(b)数值解图2㊀0.5s 的精确解与数值解(a)精确解(b)数值解图3㊀1s 的精确解与数值解77第1期㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀㊀栗雪娟,等:Cahn-Hilliard 方程的一个超紧致有限差分格式对比图,可以看出,数值解与精确解吻合很好,表明本文给出的数值格式是可行的,并且精度较高㊂3　结论本文研究了组合型超紧致差分方法和四阶Runge-Kutta方法,并将其运用于四阶Cahn-Hilliard 方程的数值求解,通过研究与分析,得到如下结论: 1)使用泰勒级数展开锁定差分格式系数,得到本文的组合型超紧致差分格式精度更高,误差更小㊂2)在边界点处有效地达到了降阶,并提高了精度㊂3)通过数值算例验证了数值格式的有效性㊂4)预估该方法可应用于高阶偏微分方程的数值求解㊂参考文献:[1]HUANG Q M,YANG J X.Linear and energy-stable method with en-hanced consistency for the incompressible Cahn-Hilliard-Navier-Stokes two-phase flow model[J].Mathematics,2022,10 (24):4711.[2]AKRIVIS G,LI B Y,LI D F.Energy-decaying extrapolated RK-SAV methods for the allen-Cahn and Cahn-Hilliard equations[J].SIAM Journal on Scientific Computing,2019,41(6):3703-3727. [3]YOUNAS U,REZAZADEH H,REN J,et al.Propagation of diverse exact solitary wave solutions in separation phase of iron(Fe-Cr-X(X =Mo,Cu))for the ternary alloys[J].International Journal of Mod-ern Physics B,2022,36(4):2250039.[4]HE R J,CHEN Z X,FENG X L.Error estimates of fully discrete finite element solution for the2D Cahn-Hilliard equation with infinite time horizon[J].Numerical Methods for Partial Differential Equati-ions,2017,33(3):742-762.[5]HE Y N,FENG X L.Uniform H2-regularity of solution for the2D Navier-Stokes/Cahn-Hilliard phase field model[J].Journal of Math-ematical Analysis and Applications,2016,441(2):815-829. [6]WEN J,HE Y N,HE Y L.Semi-implicit,unconditionally energy sta-ble,stabilized finite element method based on multiscale enrichment for the Cahn-Hilliard-Navier-Stokes phase-field model[J]. Computers and Mathematics with Applications,2022,126:172 -181.[7]MESFORUSH A,LARSSON S.A posteriori error analysis for the Cahn-Hilliard equation[J].Journal of Mathematical Modeling, 2022,10(4):437-452.[8]XIA Y,XU Y,SHU C W.Local discontinuous Galerkin methods for the Cahn-Hilliard type equation[J].Journal of Computational Phys-ics,2007,227(1):472-491.[9]CHEN L,LüS J.A fully discrete spectral scheme for time fractional Cahn-Hilliard equation with initial singularity[J].Computers and Mathematics with Applications,2022,127:213-224. [10]周诚尧,汪勇,桂志先,等.二维黏弹介质五点八阶超紧致有限差分声波方程数值模拟[J].科学技术与工程,2020,20(1):54 -63.[11]汪勇,徐佑德,高刚,等.二维黏滞声波方程的优化组合型紧致有限差分数值模拟[J].石油地球物理勘探,2018,53(6):1152 -1164,1110.[12]程晓晗,封建湖,郑素佩.求解对流扩散方程的低耗散中心迎风格式[J].应用数学,2017,30(2):344-349.(编辑:杜清玲)87山东理工大学学报(自然科学版)2024年㊀。

Finding local community structure in networksAaron Clauset*Department of Computer Science,University of New Mexico,Albuquerque,New Mexico87131,USA͑Received22February2005;published29August2005͒Although the inference of global community structure in networks has recently become a topic of great interest in the physics community,all such algorithms require that the graph be completely known.Here,we define both a measure of local community structure and an algorithm that infers the hierarchy of communities that enclose a given vertex by exploring the graph one vertex at a time.This algorithm runs in time O͑k2d͒for general graphs when d is the mean degree and k is the number of vertices to be explored.For graphs where exploring a new vertex is time consuming,the running time is linear,O͑k͒.We show that on computer-generated graphs the average behavior of this technique approximates that of algorithms that require global knowledge.As an application,we use this algorithm to extract meaningful local clustering information in the large recommender network of an online retailer.DOI:10.1103/PhysRevE.72.026132PACS number͑s͒:89.75.Hc,05.10.Ϫa,87.23.Ge,89.20.HhI.INTRODUCTIONRecently,physicists have become increasingly interested in representing the patterns of interactions in complex sys-tems as networks͓1–4͔.Canonical examples include the In-ternet͓5͔,the World Wide Web͓6͔,social networks͓7͔,ci-tation networks͓8,9͔,and biological networks͓10͔.In each case,the system is modeled as a graph with n vertices and m edges,e.g.,physical connections between computers,friend-ships between people,and citations among academic papers.Within these networks,the global organization of vertices into communities has garnered broad interest both inside and beyond the physics community.Conventionally,a commu-nity is taken to be a group of vertices in which there are more edges between vertices within the group than to vertices out-side of it.Although the partitioning of a network into such groups is a well-studied problem,older algorithms tend to only work well in special cases͓11–15͔.Several algorithms have recently been proposed within the physics community, and have been shown to reliably extract the known commu-nity structure in real world networks͓16–21͔.Similarly,the computer science community has proposed algorithms based on the concept offlow͓22͔.However,each of these algorithms require knowledge of the entire structure of the graph.This constraint is problem-atic for networks like the World Wide Web,which for all practical purposes is too large and too dynamic to ever be known fully,or networks that are larger than can be accom-modated by the fastest algorithms͓21͔.In spite of these limi-tations,we would still like to make quantitative statements about community structure,albeit confined to some acces-sible and known region of the graph in question.For in-stance,we might like to quantify the local communities of either a person given their social network or a particular website given its local topology in the World Wide Web.Here,we propose a general measure of local community structure,which we call local modularity,for graphs in which we lack global knowledge and that must be exploredone vertex at a time.We then define a fast agglomerativealgorithm that maximizes the local modularity in a greedyfashion,and test the algorithm’s performance on a series ofcomputer-generated networks with known community struc-ture.Finally,we use this algorithm to analyze the local com-munity structure of the online retailer ’s recom-mender network,which is composed of more than400000vertices and2million edges.Through this analysis,we dem-onstrate the existence of a mesoscopic network structure thatis distinct from both the microstructure of vertex statisticsand the global community structure previously given in͓21͔.Interestingly,wefind a wide variety of local communitystructures,and that generally,the local modularity of the net-work surrounding a vertex is negatively correlated with itsdegree.II.LOCAL MODULARITYThe inference of community structure can generally bereduced to identifying a partitioning of the graph that maxi-mizes some quantitative notion of community structure.However,when we lack global knowledge of the graph’stopology,a measure of community structure must necessarilybe independent of those global properties.For instance,thisrequirement precludes the use of the modularity metric Q,due to Newman and Girvan͓17͔,as it depends on m,thetotal number of edges.Suppose that in the graph G,we have perfect knowledgeof the connectivity of some set of vertices,i.e.,the knownportion of the graph,which we denote C.This necessarilyimplies the existence of a set of vertices U about which weknow only their adjacencies to C.Further,let us assume thatthe only way we may gain additional knowledge about G isby visiting some neighboring vertex v i෈U,which yields a list of its adjacencies.As a result,v i becomes a member of C,and additional unknown vertices may be added to U.Thisvertex-at-a-time discovery process is directly analogous tothe manner in which“spider”or“crawler”programs harvestthe hyperlink structure of the World Wide Web.*Electronic address:aaron@PHYSICAL REVIEW E72,026132͑2005͒The adjacency matrix of such a partially explored graph is given byA ij =Ά1,if vertices i and j are connected,and either vertex is in C ;0,otherwise.͑1͒If we consider C to constitute a local community,the most simple measure of the quality of such a partitioning of G is simply the fraction of known adjacencies that are completely internal to C .This quantity is given by͚ijA ij ␰͑i ,j ͚͒ijA ij =12m *͚ijA ij ␰͑i ,j ͒,͑2͒where m *=12͚ij A ij ,the number of edges in the partial adja-cency matrix,and ␰͑i ,j ͒is 1if both v i and v j are in C ,and 0otherwise.This quantity will be large when C has many in-ternal connections and few connections to the unknown por-tion of the graph.This measure also has the property that when ͉C ͉ӷ͉U ͉,the partition will almost always appear to be good.If we restrict our consideration to those vertices in the subset of C that have at least one neighbor in U ,i.e.,the vertices that make up the boundary of C ͑Fig.1͒,we obtain a direct measure of the sharpness of that boundary.Addition-ally,this measure is independent of the size of the enclosed community.Intuitively,we expect that a community with a sharp boundary will have few connections from its boundary to the unknown portion of the graph,while having a greater proportion of connections from the boundary back into the local community.In the interest of keeping the notation con-cise,let us denote those vertices that comprise the boundary as B ,and the boundary-adjacency matrix asB ij =Ά1,if vertices i and j are connected,and either vertex is in B ;0,otherwise.͑3͒Thus,we define the local modularity R to beR =͚ijB ij ␦͑i ,j ͚͒ijB ij=IT,͑4͒where ␦͑i ,j ͒is 1when either v i ෈B and v j ෈C ,or vice versa ,and is 0otherwise.Here,T is the number of edges with one or more end points in B ,while I is the number of those edges with neither end point in U .This measure assumes an un-weighted graph,although the weighted generalization is straightforward ͓23͔.A few comments regarding this formulation are worth-while before proceeding.By considering the fraction of boundary edges that are internal to C ,we ensure that our measure of local modularity lies on the interval 0ϽR Ͻ1,where its value is directly proportional to sharpness of the boundary given byB .This is true except when the entire component has been discovered,at which point R is unde-fined.If we like,we may set R =1in that case in order to match the intuitive notion that an entire component consti-tutes the strongest kind of community.Finally,there are cer-tainly alternative measures that can be defined on B ,how-ever,in this paper we consider only the one given.III.THE ALGORITHMFor graphs like the World Wide Web,in which one must literally crawl the network in order to discover the adjacency matrix,any analysis of local community structure must nec-essarily begin at some source vertex v 0.In general,if the explored portion of the graph has k vertices,the number of ways to partition it into two sets,those vertices considered a part of the same local community as the source vertex and those considered outside of it,is given by 2k −2−1,which is exponential in the size of the explored portion of the net-work.In this section,we describe an algorithm that only takes time polynomial in k ,and that infers local community structure by using the vertex-at-a-time discovery process subject to maximizing our measure of local modularity.Initially,we place the source vertex in the community,v 0=C ,and place its neighbors in U .At each step,the algo-rithm adds to C ͑and to B ,if necessary ͒the neighboring vertex that results in the largest increase ͑or smallest de-crease ͒in R ,breaking ties randomly.That is,for each vertex v j ෈U ,we calculate the ⌬R j that corresponds to the change in local modularity as a result of joining v j to C .The vertex that results in the largest positive change is then joined.Fi-nally,we add to U any newly discovered vertices,and update our estimate of R .This process continues until it has agglom-erated either a given number of vertices k ,or it has discov-ered the entire enclosing component,whichever happens first.Generally,these operations can be done quite quickly by an appropriate use of arrays and linked lists.The pseudocode for this process is given in Table I.As we will see in the two subsequent sections,this algorithm performs well on both computer-generated graphs with some known community structure and on real world graphs in which we can manually verify the sensibility of the inferred communi-ties.FIG.1.An illustration of the division of an abstract graph into the local community C ,its boundary B ,and the edges that connect B to the largely unknown neighbors U .AARON CLAUSET PHYSICAL REVIEW E 72,026132͑2005͒The computation of the⌬R j associated with each v j෈Ucan be done quickly using an expression derived from Eq.͑4͒:⌬R j=x−Ry−z͑1−R͒T−z+y,͑5͒where x is the number of edges in T that terminated at v j,y is the number of edges that will be added to T by the agglom-eration of v j͑i.e.,the degree of v j is x+y͒,and z is the number of edges that will be removed from T by the agglom-eration.Because⌬R j depends on the current value of R,and on the y and z that correspond to v j,each step of the algo-rithm takes time proportional to the number of vertices in U. This is roughly kd,where d is the mean degree of the graph; we note that this will be a significant overestimate for graphs with nontrivial clustering coefficients,significant community structure,or when C is a large portion of the graph.Thus,in general,the running time for the algorithm is O͑k2d͒,or simply O͑k2͒for a sparse graph,i.e.,when mϳn.As it ag-glomerates vertices,the algorithm outputs a function R͑t͒, the local modularity of the community centered on v0after t steps,and a list of vertices paired with the time of their agglomeration.The above calculation of the running time is somewhat misleading,as it assumes that the algorithm is dominated by the time required to calculate the⌬R j for each vertex in U; however,for graphs like the World Wide Web,where adding a new vertex to U requires the algorithm to fetch a web page from a remote server,the running time will instead be domi-nated by the time-consuming retrieval.When this is true,the running time is linear in the size of the explored subgraph, O͑k͒.A few comments regarding this algorithm are due.Be-cause we greedily maximize the local modularity,a neigh-boring high degree vertex will only be agglomerated if a sufficient number of its neighbors have been explored.It is this behavior that prevents the algorithm from crossing a community boundary until absolutely necessary.Addition-ally,the algorithm is somewhat sensitive to the degree dis-tribution of the source vertex’s neighbors:when the source degree is high,the algorithm willfirst explore its low degree neighbors.This implicitly assumes that high degree vertices are likely to sit at the boundary of several local communities. While certainly not the case in general,this may be true for some real world networks.We shall return to this idea in a later section.Finally,although one could certainly stop the algorithm once thefirst enclosing community has been found,in prin-ciple,there is no reason that it cannot continue until some arbitrary number of vertices have been agglomerated.Doing so yields the hierarchy of communities that enclose the source vertex.In a sense,this process is akin to the follow-ing:given the dendrogram of the global community hierar-chy,walk upward toward the root from some leaf v0and observe the successive hierarchical relationships,as repre-sented by junctions in the dendrogram.Thus,the enclosing communities inferred by our algorithm for some source ver-tex is the community hierarchy from the perspective of that vertex.PUTER-GENERATED GRAPHS As has become standard with testing community infer-ence techniques,we apply our algorithm to a set of computer-generated random graphs that have known com-munity structure͓17͔.In these graphs,n=128vertices are divided into four equal-sized communities of32vertices. Each vertex has a total expected degree z that is divided between intra-and intercommunity edges such that z=z in+z out.These edges are placed independently and at ran-dom so that,in expectation,the values of z in and z out are respected.By holding the expected degree constant z=16, we may tune the sharpness of the community boundaries by varying z out.Note that for these graphs,when z out=12,edges between vertices in the same group are just as likely as edges between vertices that are not.Figure2shows the average local modularity R as a func-tion of the number of steps t,over500realizations of the graphs described above,with a randomly selected sourceTABLE I.The general algorithm for the greedy maximization of local modularity,as given in the text.addv0to Cadd all neighbors of v0to Uset B=v0while͉C͉Ͻk dofor each v j෈U docompute⌬R j from Eq.͑5͒end forfind v j such that its⌬R j is maximum add that v j to Cadd all new neighbors of that v j to Uupdate R and B end whileFIG.2.͑Color online͒Average local modularity R as a function of the number of agglomerations t͑details described in the text; error bars are omitted for clarity͒.By varying the expected number of intercommunity edges per node z out,the strength of the commu-nity boundaries are varied.FINDING LOCAL COMMUNITY STRUCTURE IN NETWORKS PHYSICAL REVIEW E72,026132͑2005͒vertex.For the sake of clarity,only data series for z out ഛ6.0are shown and error bars are omitted.Sharp community boundaries correspond to peaks in the curve.As z out grows,the sharpness of the boundaries and the height of the peaks decrease proportionally.When the first derivative is positive everywhere,e.g.,for z out Ͼ5,the inferred locations of the community boundaries may be extracted by finding local minima in the second derivative,possibly after some smoothing.From this information we may grade the performance of the algorithm on the computer-generated graphs in the fol-lowing manner:we label each vertex with its true community membership,run the local inference algorithm with a ran-domly selected source vertex,detect the community bound-aries as described above,and compute the fraction of cor-rectly classified vertices for each community.Averaging over 500instantiations provides an estimate,as a function of z out ,of the mean fraction of correctly classified vertices for each community ͑Fig.3͒;error bars depict one standard deviation.As a method for inferring the first enclosing community,our algorithm correctly classifies,on average,more than 50%of the vertices even when the boundaries are weak,i.e.,when z out =8.We note that for such computer-generated graphs,our algorithm’s average performance approximates that of more global methods ͓17–19͔,although the variance in the quality of classification is substantially larger for weak boundaries.Indeed,such an increase is to be expected as the algorithm we describe here may be misled by large fluctua-tions in local structure,while global algorithms may not.Recently,another approach to inferring community struc-ture using only local information appeared ͓24͔.This alter-native technique relies upon growing a breadth-first tree out-ward from the source vertex v 0,until the rate of expansion falls below an arbitrary threshold.The uniform exploration has the property that some level in the tree will correspond to a good partitioning only when v 0is equidistant from all parts of its enclosing community’s boundary.On the other hand,by exploring the surrounding graph one vertex at a time,our algorithm will avoid crossing boundaries until it has ex-plored the remainder of the enclosing community.V .LOCAL COPURCHASING HABITSIn this section,we apply our local inference algorithm to the recommender network of ,collected in Au-gust 2003,which has n =409687vertices,m =2464630edges,and thus a mean degree of 12.03.We note that the degree distribution is fairly right-skewed,having a standard deviation of 14.64.Here,vertices are items such as books and digital media sold on Amazon’s website,while edges connect pairs of items that are frequently purchased together by customers.It is this copurchasing data that yields recom-mendations for customers as they browse the online store.Although,in general,the algorithm we have described is intended for graphs like the World Wide Web,the Amazon recommender network has the advantage that,by virtue of being both very large and fully known,we may explore glo-bal regularities in local community structure without concern for sampling bias in the choice of source vertices.Addition-ally,we may check the inferred the community structures against our,admittedly heuristic,notions of correctness.As illustrative examples,we choose three qualitatively different items as source vertices:the compact disk Alegria by Cirque du Soleil,the book Small Worlds and the book Harry Potter and the Order of the Phoenix by J.K Rowling.These items have degree 15,19,and 3117,respectively.At the time the network data was collected,the Harry Potter book was the highest degree vertex in the network,its release date having been June 2003.For each of these items,we explore k =25000associated vertices.Figure 4illustrates the local modularity as a function of the number of steps t for each item;an analogous data series for a random graph with the same degree distribution ͓25͔has been plotted for com-parison.We mark the locations of the five principle enclosing communities with large open symbols,where their location was determined by the same second derivative test described above.These time series have several distinguishing features.First,Alegria has the smallest enclosing communities,com-posed of t =͕10,30,39,58,78͖vertices,and thesecommuni-FIG.3.͑Color online ͒Fraction of correctly classified nodes,by community,as a function of the number of inter-community edges z out .Although there the variance increases as the community bound-aries become less sharp,the average behavior degradesgracefully.FIG.4.͑Color online ͒Local modularity R for three items in the recommender network,shown on log-linear axes.For comparison,the time series for a random graph with the same de-gree distribution is shown.The large open symbols indicate the locations of the five strongest enclosing communities.AARON CLAUSET PHYSICAL REVIEW E 72,026132͑2005͒ties are associated with high values of local modularity.The first five enclosing communities all have R Ͼ0.62,while the third community corresponds to R =0.81,indicating that only about 20%of the boundary edges reach out to the rest of the network.In contrast,the communities of Small Worlds con-tain t =͕36,48,69,82,94͖vertices,while the Harry Potter book’s communities are extremely large,containing t =͕607,883,1270,1374,1438͖vertices.Both sets have only moderate values of local modularity,R ഛ0.43.It is notable that the local modularity functions for all three items follow relatively distinct trajectories until the algorithm has agglom-erated roughly 10000items.Beyond that point,the curves begin to converge,indicating that,from the perspectives of the source vertices,the local community structure has be-come relatively similar.To illustrate the inferred local structure,we show the par-tial subgraph that corresponds to the first three enclosing local communities for the compact disc Alegria in Fig.5.Here,communities are distinguished by shape according to the order of discovery ͑circle,diamond,and square,respec-tively ͒,and vertices beyond these communities are denoted by triangles.Items in the first enclosing community are uni-formly compact disks produced by Cirque du Soleil.Items in the second are slightly more diverse,including movies and books about the troupe,the Cirque du Soleil compact disk entitled Varekai ,and one compact disk by a band called Era;the third group contains both new and old Cirque du Soleil movies.Varekai appears to have been placed outside the first community because it has fewer connections to those items than to items in the subsequent enclosing communities.Briefly,we find that the enclosing local communities of Small Worlds are populated by texts in sociology and social network analysis,while the Harry Potter book’s communities have little topical similarity.In Fig.5,the labels 1and 4denote the items Alegria and Order of the Phoenix ,respectively.It is notable that these items are only three steps away in the graph,yet have ex-tremely different local community structures ͑Fig.4͒.If an item’s popularity is reflected by its degree,then it is reason-able to believe that the strength of the source vertex’s local community structure may be inversely related to its degree.That is,popular items like Order of the Phoenix may connect to many well-defined communities by virtue of being pur-chased by a large number of customers with diverse interests,while niche items like Cirque du Soleil’s Alegria exhibit stronger local community structure from more specific copurchasing habits.Such a structure appears to be distinct from both the macroscopic structure discovered using global community inference methods ͓21͔and the microscopic structure of simple vertex statistics such as the clustering coefficient or assortative mixing,i.e.,topological correla-tions based on degree or vertex-type similarity ͓26͔.We call this intermediate level of topological structure “mesoscopic.”With the exception of social networks,the degree of ad-jacent vertices in most complex networks appears to be nega-tively correlated in most networks.This property is often called “disassortative”mixing ͓26͔,and can be caused by a high clustering coefficient,global community structure,or a specific social mechanism ͓27͔.However,for the Amazon recommender network,we find that the assortativity coefficient is statistically no different from zero,r =−3.01ϫ10−19±1.49ϫ10−4,yet the network exhibits a nontrivial clustering coefficient,c =0.17and strong global community structure with a peak modularity of Q =0.745͓21͔.Returning to the suggested inverse relationship between the degree of the source vertex and the strength of its sur-rounding community structure,we sample for 100000ran-dom vertices the average local modularity over the first k =250steps.We find the average local modularity to be relatively high,͗R amzn ͘=0.49±0.08,while a random graph with the same degree distribution yields ͗R rand ͘=0.16±0.01.The variance for the Amazon graph is due to the contribu-tions of high degree vertices.In Fig.6,we plot from our random sample the average local modularity for all source vertices with degree of at least d .Notably,the averageisFIG.5.The first three enclosing communities for Cirque du Soleil’s Alegria in ’s recommender network;commu-nities are distinguished by shape ͑circles,diamonds,squares,re-spectively ͒.Connections to triangles represent connections to items in the remaining unknown portion of the graph.Alegria and Order of the Phoenix are denoted by 1and 4,respectively.FIG.6.͑Color online ͒The average local modularity over the first 250steps for source vertices with degree at least d .The “knee”in the upper data series is located at d =13;the mean degree for the network is 12.03.The logarithmic falloff illustrates the negative correlation between the source vertex degree and the strength of the surrounding local community.FINDING LOCAL COMMUNITY STRUCTURE IN NETWORKS PHYSICAL REVIEW E 72,026132͑2005͒relatively constant until d=13,after which it falls off loga-rithmically.This supports the hypothesis that,in the recom-mender network,there is a weak inverse relationship be-tween the degree of the source vertex and the strength of itssurrounding local community.Naturally,there are many ways to use the concept of localcommunity structure to understand the properties of realworld networks.Further characterizations of the Amazongraph are beyond the scope of this paper,but we propose arigorous exploration of the relationship between the sourcevertex degree and its surrounding local community structureas a topic for future work.VI.CONCLUSIONSAlthough many recent algorithms have appeared in thephysics literature for the inference of community structurewhen the entire graph structure is known,there has beenlittle consideration of graphs that are either too large for eventhe fastest known techniques,or that are,like the WorldWide Web,too large or too dynamic to ever be fully known.Here,we define a measure of community structure that de-pends only on the topology of some known portion of agraph.We then give a simple fast,agglomerative algorithmthat greedily maximizes our measure as it explores the graphone vertex at a time.When the time it takes to retrieve theadjacencies of a vertex is small,this algorithm runs in timeO͑k2d͒for general graphs when it explores k vertices,and the graph has mean degree d.For sparse graphs,i.e.,whenmϳn,this is simply O͑k2͒.On the other hand,when visiting a new vertex to retrieve its adjacencies dominates the run-ning time,e.g.,downloading a web page on the World WideWeb,the algorithm takes time linear in the size of the ex-plored subgraph,O͑k͒.Generally,if we are interested in making quantitative statements about local structure,that is, when kӶn,it is much more reasonable to use an algorithm which is linear or even quadratic in k,than an algorithm that is linear in the size of the graph n.Finally,we note that our algorithm’s simplicity will make it especially easy to inco-porate into web spider or crawler programs for the discovery of local community structures on the World Wide Web graph.Using computer-generated graphs with known community structure,we show that our algorithm extracts this structure and that its average performance approximates that of more global inference algorithms͓17–19͔that rely on global infor-mation.We then apply our algorithm to the large recom-mender network of the online retailer ,and ex-tract the local hierarchy of communities for several qualitatively distinct items.Additionally,we show that a ver-tex’s degree is inversely related to the strength of its sur-rounding local structure.This discovery points to the exis-tence of topological structure in real world networks that is above the single-vertex structure such as the clustering coef-ficient,but below that of the global community structure. Finally,this algorithm should allow researchers to character-ize the structure of a wide variety of other graphs,such as the World Wide Web,and we look forward to seeing such appli-cations.ACKNOWLEDGMENTSThe author is grateful to Cristopher Moore for persistent guidance,Mark Newman for discussions about community structure,and and Eric Promislow for provid-ing the recommender network data.This work was supported in part by the National Science Foundation under Grants No. PHY-0200909and No.ITR-0324845.͓1͔S.H.Strogatz,Nature͑London͒410,268͑2001͒.͓2͔R.Albert and A.-L.Barabási,Rev.Mod.Phys.74,47͑2002͒.͓3͔S.N.Dorogovtsev and J.F.F.Mendes,Adv.Phys.51,1079͑2002͒.͓4͔M.E.J.Newman,SIAM Rev.45,167͑2003͒.͓5͔M.Faloutsos,P.Faloutsos,and C.Faloutsos,-mun.Rev.29,251͑1999͒.͓6͔J.M.Kleinberg,S.R.Kumar,P.Raghavan,S.Rajagopalan, and A.Tomkins,“The Web as a graph:Measurements,models and methods,”in Proceedings of the International Conference on Combinatorics and Computing,No.1627in Lecture Notes in Computer Science͑Springer-Verlag,Berlin,1999͒,pp.1–18.͓7͔S.Wasserman and K.Faust,Social Network Analysis͑Cam-bridge University Press,Cambridge,1994͒.͓8͔D.J.de S.Price,Science149,510͑1965͒.͓9͔S.Redner,e-print physics/0407137.͓10͔T.Ito,T.Chiba,R.Ozawa,M.Yoshida,M.Hattori,and Y.Sakaki,Proc.Natl.Acad.Sci.U.S.A.98,4569͑2001͒.͓11͔B.W.Kernighan and S.Lin,Bell Syst.Tech.J.49,291͑1970͒.͓12͔M.Fiedler,Czech.Math.J.23,298͑1973͒.͓13͔A.Pothen,H.Simon,and K.-P.Liou,SIAM J.Matrix Anal.Appl.11,430͑1990͒.͓14͔J.Scott,Social Network Analysis:A Handbook,2nd ed.,Sage, London͑2000͒.͓15͔M.E.J.Newman,Eur.Phys.J.B38,321͑2004͒.͓16͔M.Girvan and M.E.J.Newman,Proc.Natl.Acad.Sci.U.S.A.99,7821͑2002͒.͓17͔M.E.J.Newman,Phys.Rev.E69,026113͑2004͒.͓18͔F.Radicchi,C.Castellano,F.Cecconi,V.Loreto,and D.Pa-risi,Proc.Natl.Acad.Sci.U.S.A.101,2658͑2004͒.͓19͔M.E.J.Newman,Phys.Rev.E69,066133͑2004͒.͓20͔F.Wu and B.A.Huberman,Eur.Phys.J.B38,331͑2004͒.͓21͔A.Clauset,M.E.J.Newman,and C.Moore,Phys.Rev.E70, 066111͑2004͒.͓22͔G.W.Flake,wrence,C.L.Giles,and F.M.Coetzee, IEEE Computer35͑3͒,66͑2002͒.͓23͔M.E.J.Newman,Phys.Rev.E70,056131͑2004͒.͓24͔J.P.Bagrow and E.M.Bollt,e-print cond-mat/0412482.͓25͔B.Bollobás,Random Graphs͑Academic,New York,1985͒.͓26͔M.E.J.Newman,Phys.Rev.E67,026126͑2003͒.͓27͔M.E.J.Newman and J.Park,Phys.Rev.E68,036122͑2003͒.AARON CLAUSET PHYSICAL REVIEW E72,026132͑2005͒。

英文回答：The Newton-Raphson method is an iterative optimization algorithm utilized for locating the local minimum or maximumof a given function. Within the realm of optimization, the Newton-Raphson method iteratively updates the current solution by leveraging the second derivative information of the objective function. This approach enables the method to converge towards the optimal solution at an accelerated pacepared to first-order optimization algorithms, such as the gradient descent method. Nonetheless, the Newton-Raphson method necessitates the solution of a system of linear equations involving the Hessian matrix, which denotes the second derivative of the objective function. Of particular note, when the Hessian matrix possesses a rank of one, it introduces a special case for the Newton-Raphson method.牛顿—拉弗森方法是一种迭代优化算法，用于定位特定函数的局部最小或最大值。

2011THOMAS J. SARGENT , CHRISTOPHER A. SIMS托马斯·J·萨金特克里斯托弗·A·西姆斯for their empirical research on cause and effect in the macroeconomy.两名获奖者的研究成果解答了许多有关经济政策与宏观经济变量之间的关系问题，例如提高利率或减税将对国内生产总值和通货膨胀产生何种影响，中央银行调整通货膨胀目标将产生何种后果等。

2010PETER A. DIAMOND DALE T. MORTENSEN CHRISTOPHER A. PISSARIDES彼得·戴蒙德戴尔·莫特森克里斯托弗·皮萨里德斯for their analysis of markets with search frictions.瑞典皇家科学院表示，他们对市场的分析使其可以得到这个奖项。

“市场大部分交易都是为贸易而进行的，当然会出现一些贸易摩擦，买者很难得到想要买的买品，而卖者很难找到消费者。

在劳动力市场上许多公司也发现会有许多工作空缺，而一些失业人员找不到适合的工作岗位。

”彼特-戴蒙德等人所开发的理论是解释了市场上这种冲突，他们的理论是基于微观经济学理论的，也就是市场合理产出，他们的工作也就是意味着雇佣工人要更加合理，在招聘人员和需求工作应该提供合理的机制。

2009ELINOR OSTROM OLIVER E. WILLIAMSON埃莉诺·奥斯特罗姆奥利弗·E·威廉姆森for her analysis of economic governance, especially the commons表彰她对经济治理的分析，尤其是对普通人经济治理活动的研究for his analysis of economic governance, especially the boundaries of the firm表彰他对经济治理的分析，特别是对公司的经济治理边界的分析2008PAUL KRUGMAN保罗·克鲁格曼for his analysis of trade patterns and location of economic activity.在贸易模式上所做的分析工作和对经济活动的定位2007LEONID HURWICZ ERIC S. MASKIN ROGER B. MYERSON莱昂尼德·赫维奇埃里克·马斯金罗杰·迈尔森for having laid the foundations of mechanism design theory.为机制设计理论（Mechanism Design Theory）奠定了基础2006EDMUND S. PHELPS埃德蒙德·菲尔普斯for his analysis of intertemporal tradeoffs in macroeconomic policy.加深人们对于通货膨胀和失业预期关系的理解方面所做的贡献2005ROBERT J. AUMANN THOMAS C. SCHELLING罗伯特·约翰·奥曼托马斯·克罗姆比·谢林for having enhanced our understanding of conflict and cooperation through game-theory analysis.通过博弈论分析促进了对冲突与合作的理解2004FINN E. KYDLAND EDWARD C. PRESCOTT芬恩·基德兰德爱德华·普雷斯科特for their contributions to dynamic macroeconomics: the time consistency of economic policy and the driving forces behind business cycles通过对宏观经济政策运用中“时间连贯性难题”的分析研究，为经济政策特别是货币政策的实际有效运用提供了思路；二是在对商业周期的研究中，通过对引起商业周期波动的各种因素和各因素间相互关系的分析，使人们对于这一现象的认识更加深入2003ROBERT F. ENGLE CLIVE W. J. GRANGER罗伯特·恩格尔克莱夫·格兰杰for methods of analyzing economic time series with time-varying volatility (ARCH)表彰罗伯特·恩格尔使用“随着时间变化的易变性”新方法分析经济时间序列，从而给经济学研究和经济发展带来巨大影响for methods of analyzing economic time series with common trends (cointegration)表彰克莱夫·格兰杰使用“共同趋势”新方法分析经济时间序列，从而给经济学研究和经济发展带来巨大影响2002DANIEL KAHNEMAN VERNON L. SMITH,丹尼尔·卡纳曼弗农·史密斯for having integrated insights from psychological research into economic science, especially concerning human judgment and decision-making under uncertainty表彰丹尼尔·卡纳曼用认知心理学分析方法研究人类的判断和决策行为的领域for having established laboratory experiments as a tool in empirical economic analysis, especially in the study of alternative market mechanisms表彰弗农·史密斯通过实验室实验来测试或检验根据经济学理论作出预测的未知或不确定性领域2001GEORGE A. AKERLOF A. MICHAEL SPENCE JOSEPH E. STIGLITZ乔治·阿克尔洛夫迈克尔·斯宾塞约瑟夫·斯蒂格利茨for their analyses of markets with asymmetric information.为不对称信息市场的一般理论奠定了基石2000JAMES J. HECKMAN DANIEL L. MCFADDEN詹姆斯·赫克曼丹尼尔·麦克法登for his development of theory and methods for analyzing selective samples表彰詹姆斯·赫克曼对分析选择性抽样的原理和方法所做出的发展和贡献for his development of theory and methods for analyzing discrete choice.表彰丹尼尔·麦克法登对分析离散选择的原理和方法所做出的发展和贡献1999ROBERT A. MUNDELL罗伯特·蒙代尔for his analysis of monetary and fiscal policy under different exchange rate regimes and his analysis of optimum currency areas.表彰他对不同汇率体制下货币与财政政策以及最适宜的货币流通区域所做的分析1998AMATYA SEN阿马蒂亚·森for his contributions to welfare economics.对福利经济学几个重大问题做出了贡献，包括社会选择理论、对福利和贫穷标准的定义、对匮乏的研究等1997ROBERT C. MERTON MYRON S. SCHOLES罗伯特·默顿迈伦·斯科尔斯for a new method to determine the value of derivatives.前者对布莱克－斯科尔斯公式所依赖的假设条件做了进一步减弱，在许多方面对其做了推广。

Network impacts of a road capacity reduction:Empirical analysisand model predictionsDavid Watling a ,⇑,David Milne a ,Stephen Clark baInstitute for Transport Studies,University of Leeds,Woodhouse Lane,Leeds LS29JT,UK b Leeds City Council,Leonardo Building,2Rossington Street,Leeds LS28HD,UKa r t i c l e i n f o Article history:Received 24May 2010Received in revised form 15July 2011Accepted 7September 2011Keywords:Traffic assignment Network models Equilibrium Route choice Day-to-day variabilitya b s t r a c tIn spite of their widespread use in policy design and evaluation,relatively little evidencehas been reported on how well traffic equilibrium models predict real network impacts.Here we present what we believe to be the first paper that together analyses the explicitimpacts on observed route choice of an actual network intervention and compares thiswith the before-and-after predictions of a network equilibrium model.The analysis isbased on the findings of an empirical study of the travel time and route choice impactsof a road capacity reduction.Time-stamped,partial licence plates were recorded across aseries of locations,over a period of days both with and without the capacity reduction,and the data were ‘matched’between locations using special-purpose statistical methods.Hypothesis tests were used to identify statistically significant changes in travel times androute choice,between the periods of days with and without the capacity reduction.A trafficnetwork equilibrium model was then independently applied to the same scenarios,and itspredictions compared with the empirical findings.From a comparison of route choice pat-terns,a particularly influential spatial effect was revealed of the parameter specifying therelative values of distance and travel time assumed in the generalised cost equations.When this parameter was ‘fitted’to the data without the capacity reduction,the networkmodel broadly predicted the route choice impacts of the capacity reduction,but with othervalues it was seen to perform poorly.The paper concludes by discussing the wider practicaland research implications of the study’s findings.Ó2011Elsevier Ltd.All rights reserved.1.IntroductionIt is well known that altering the localised characteristics of a road network,such as a planned change in road capacity,will tend to have both direct and indirect effects.The direct effects are imparted on the road itself,in terms of how it can deal with a given demand flow entering the link,with an impact on travel times to traverse the link at a given demand flow level.The indirect effects arise due to drivers changing their travel decisions,such as choice of route,in response to the altered travel times.There are many practical circumstances in which it is desirable to forecast these direct and indirect impacts in the context of a systematic change in road capacity.For example,in the case of proposed road widening or junction improvements,there is typically a need to justify econom-ically the required investment in terms of the benefits that will likely accrue.There are also several examples in which it is relevant to examine the impacts of road capacity reduction .For example,if one proposes to reallocate road space between alternative modes,such as increased bus and cycle lane provision or a pedestrianisation scheme,then typically a range of alternative designs exist which may differ in their ability to accommodate efficiently the new traffic and routing patterns.0965-8564/$-see front matter Ó2011Elsevier Ltd.All rights reserved.doi:10.1016/j.tra.2011.09.010⇑Corresponding author.Tel.:+441133436612;fax:+441133435334.E-mail address:d.p.watling@ (D.Watling).168 D.Watling et al./Transportation Research Part A46(2012)167–189Through mathematical modelling,the alternative designs may be tested in a simulated environment and the most efficient selected for implementation.Even after a particular design is selected,mathematical models may be used to adjust signal timings to optimise the use of the transport system.Road capacity may also be affected periodically by maintenance to essential services(e.g.water,electricity)or to the road itself,and often this can lead to restricted access over a period of days and weeks.In such cases,planning authorities may use modelling to devise suitable diversionary advice for drivers,and to plan any temporary changes to traffic signals or priorities.Berdica(2002)and Taylor et al.(2006)suggest more of a pro-ac-tive approach,proposing that models should be used to test networks for potential vulnerability,before any reduction mate-rialises,identifying links which if reduced in capacity over an extended period1would have a substantial impact on system performance.There are therefore practical requirements for a suitable network model of travel time and route choice impacts of capac-ity changes.The dominant method that has emerged for this purpose over the last decades is clearly the network equilibrium approach,as proposed by Beckmann et al.(1956)and developed in several directions since.The basis of using this approach is the proposition of what are believed to be‘rational’models of behaviour and other system components(e.g.link perfor-mance functions),with site-specific data used to tailor such models to particular case studies.Cross-sectional forecasts of network performance at specific road capacity states may then be made,such that at the time of any‘snapshot’forecast, drivers’route choices are in some kind of individually-optimum state.In this state,drivers cannot improve their route selec-tion by a unilateral change of route,at the snapshot travel time levels.The accepted practice is to‘validate’such models on a case-by-case basis,by ensuring that the model—when supplied with a particular set of parameters,input network data and input origin–destination demand data—reproduces current mea-sured mean link trafficflows and mean journey times,on a sample of links,to some degree of accuracy(see for example,the practical guidelines in TMIP(1997)and Highways Agency(2002)).This kind of aggregate level,cross-sectional validation to existing conditions persists across a range of network modelling paradigms,ranging from static and dynamic equilibrium (Florian and Nguyen,1976;Leonard and Tough,1979;Stephenson and Teply,1984;Matzoros et al.,1987;Janson et al., 1986;Janson,1991)to micro-simulation approaches(Laird et al.,1999;Ben-Akiva et al.,2000;Keenan,2005).While such an approach is plausible,it leaves many questions unanswered,and we would particularly highlight two: 1.The process of calibration and validation of a network equilibrium model may typically occur in a cycle.That is to say,having initially calibrated a model using the base data sources,if the subsequent validation reveals substantial discrep-ancies in some part of the network,it is then natural to adjust the model parameters(including perhaps even the OD matrix elements)until the model outputs better reflect the validation data.2In this process,then,we allow the adjustment of potentially a large number of network parameters and input data in order to replicate the validation data,yet these data themselves are highly aggregate,existing only at the link level.To be clear here,we are talking about a level of coarseness even greater than that in aggregate choice models,since we cannot even infer from link-level data the aggregate shares on alternative routes or OD movements.The question that arises is then:how many different combinations of parameters and input data values might lead to a similar link-level validation,and even if we knew the answer to this question,how might we choose between these alternative combinations?In practice,this issue is typically neglected,meaning that the‘valida-tion’is a rather weak test of the model.2.Since the data are cross-sectional in time(i.e.the aim is to reproduce current base conditions in equilibrium),then in spiteof the large efforts required in data collection,no empirical evidence is routinely collected regarding the model’s main purpose,namely its ability to predict changes in behaviour and network performance under changes to the network/ demand.This issue is exacerbated by the aggregation concerns in point1:the‘ambiguity’in choosing appropriate param-eter values to satisfy the aggregate,link-level,base validation strengthens the need to independently verify that,with the selected parameter values,the model responds reliably to changes.Although such problems–offitting equilibrium models to cross-sectional data–have long been recognised by practitioners and academics(see,e.g.,Goodwin,1998), the approach described above remains the state-of-practice.Having identified these two problems,how might we go about addressing them?One approach to thefirst problem would be to return to the underlying formulation of the network model,and instead require a model definition that permits analysis by statistical inference techniques(see for example,Nakayama et al.,2009).In this way,we may potentially exploit more information in the variability of the link-level data,with well-defined notions(such as maximum likelihood)allowing a systematic basis for selection between alternative parameter value combinations.However,this approach is still using rather limited data and it is natural not just to question the model but also the data that we use to calibrate and validate it.Yet this is not altogether straightforward to resolve.As Mahmassani and Jou(2000) remarked:‘A major difficulty...is obtaining observations of actual trip-maker behaviour,at the desired level of richness, simultaneously with measurements of prevailing conditions’.For this reason,several authors have turned to simulated gaming environments and/or stated preference techniques to elicit information on drivers’route choice behaviour(e.g. 1Clearly,more sporadic and less predictable reductions in capacity may also occur,such as in the case of breakdowns and accidents,and environmental factors such as severe weather,floods or landslides(see for example,Iida,1999),but the responses to such cases are outside the scope of the present paper. 2Some authors have suggested more systematic,bi-level type optimization processes for thisfitting process(e.g.Xu et al.,2004),but this has no material effect on the essential points above.D.Watling et al./Transportation Research Part A46(2012)167–189169 Mahmassani and Herman,1990;Iida et al.,1992;Khattak et al.,1993;Vaughn et al.,1995;Wardman et al.,1997;Jou,2001; Chen et al.,2001).This provides potentially rich information for calibrating complex behavioural models,but has the obvious limitation that it is based on imagined rather than real route choice situations.Aside from its common focus on hypothetical decision situations,this latter body of work also signifies a subtle change of emphasis in the treatment of the overall network calibration problem.Rather than viewing the network equilibrium calibra-tion process as a whole,the focus is on particular components of the model;in the cases above,the focus is on that compo-nent concerned with how drivers make route decisions.If we are prepared to make such a component-wise analysis,then certainly there exists abundant empirical evidence in the literature,with a history across a number of decades of research into issues such as the factors affecting drivers’route choice(e.g.Wachs,1967;Huchingson et al.,1977;Abu-Eisheh and Mannering,1987;Duffell and Kalombaris,1988;Antonisse et al.,1989;Bekhor et al.,2002;Liu et al.,2004),the nature of travel time variability(e.g.Smeed and Jeffcoate,1971;Montgomery and May,1987;May et al.,1989;McLeod et al., 1993),and the factors affecting trafficflow variability(Bonsall et al.,1984;Huff and Hanson,1986;Ribeiro,1994;Rakha and Van Aerde,1995;Fox et al.,1998).While these works provide useful evidence for the network equilibrium calibration problem,they do not provide a frame-work in which we can judge the overall‘fit’of a particular network model in the light of uncertainty,ambient variation and systematic changes in network attributes,be they related to the OD demand,the route choice process,travel times or the network data.Moreover,such data does nothing to address the second point made above,namely the question of how to validate the model forecasts under systematic changes to its inputs.The studies of Mannering et al.(1994)and Emmerink et al.(1996)are distinctive in this context in that they address some of the empirical concerns expressed in the context of travel information impacts,but their work stops at the stage of the empirical analysis,without a link being made to net-work prediction models.The focus of the present paper therefore is both to present thefindings of an empirical study and to link this empirical evidence to network forecasting models.More recently,Zhu et al.(2010)analysed several sources of data for evidence of the traffic and behavioural impacts of the I-35W bridge collapse in Minneapolis.Most pertinent to the present paper is their location-specific analysis of linkflows at 24locations;by computing the root mean square difference inflows between successive weeks,and comparing the trend for 2006with that for2007(the latter with the bridge collapse),they observed an apparent transient impact of the bridge col-lapse.They also showed there was no statistically-significant evidence of a difference in the pattern offlows in the period September–November2007(a period starting6weeks after the bridge collapse),when compared with the corresponding period in2006.They suggested that this was indicative of the length of a‘re-equilibration process’in a conceptual sense, though did not explicitly compare their empiricalfindings with those of a network equilibrium model.The structure of the remainder of the paper is as follows.In Section2we describe the process of selecting the real-life problem to analyse,together with the details and rationale behind the survey design.Following this,Section3describes the statistical techniques used to extract information on travel times and routing patterns from the survey data.Statistical inference is then considered in Section4,with the aim of detecting statistically significant explanatory factors.In Section5 comparisons are made between the observed network data and those predicted by a network equilibrium model.Finally,in Section6the conclusions of the study are highlighted,and recommendations made for both practice and future research.2.Experimental designThe ultimate objective of the study was to compare actual data with the output of a traffic network equilibrium model, specifically in terms of how well the equilibrium model was able to correctly forecast the impact of a systematic change ap-plied to the network.While a wealth of surveillance data on linkflows and travel times is routinely collected by many local and national agencies,we did not believe that such data would be sufficiently informative for our purposes.The reason is that while such data can often be disaggregated down to small time step resolutions,the data remains aggregate in terms of what it informs about driver response,since it does not provide the opportunity to explicitly trace vehicles(even in aggre-gate form)across more than one location.This has the effect that observed differences in linkflows might be attributed to many potential causes:it is especially difficult to separate out,say,ambient daily variation in the trip demand matrix from systematic changes in route choice,since both may give rise to similar impacts on observed linkflow patterns across re-corded sites.While methods do exist for reconstructing OD and network route patterns from observed link data(e.g.Yang et al.,1994),these are typically based on the premise of a valid network equilibrium model:in this case then,the data would not be able to give independent information on the validity of the network equilibrium approach.For these reasons it was decided to design and implement a purpose-built survey.However,it would not be efficient to extensively monitor a network in order to wait for something to happen,and therefore we required advance notification of some planned intervention.For this reason we chose to study the impact of urban maintenance work affecting the roads,which UK local government authorities organise on an annual basis as part of their‘Local Transport Plan’.The city council of York,a historic city in the north of England,agreed to inform us of their plans and to assist in the subsequent data collection exercise.Based on the interventions planned by York CC,the list of candidate studies was narrowed by considering factors such as its propensity to induce significant re-routing and its impact on the peak periods.Effectively the motivation here was to identify interventions that were likely to have a large impact on delays,since route choice impacts would then likely be more significant and more easily distinguished from ambient variability.This was notably at odds with the objectives of York CC,170 D.Watling et al./Transportation Research Part A46(2012)167–189in that they wished to minimise disruption,and so where possible York CC planned interventions to take place at times of day and of the year where impacts were minimised;therefore our own requirement greatly reduced the candidate set of studies to monitor.A further consideration in study selection was its timing in the year for scheduling before/after surveys so to avoid confounding effects of known significant‘seasonal’demand changes,e.g.the impact of the change between school semesters and holidays.A further consideration was York’s role as a major tourist attraction,which is also known to have a seasonal trend.However,the impact on car traffic is relatively small due to the strong promotion of public trans-port and restrictions on car travel and parking in the historic centre.We felt that we further mitigated such impacts by sub-sequently choosing to survey in the morning peak,at a time before most tourist attractions are open.Aside from the question of which intervention to survey was the issue of what data to collect.Within the resources of the project,we considered several options.We rejected stated preference survey methods as,although they provide a link to personal/socio-economic drivers,we wanted to compare actual behaviour with a network model;if the stated preference data conflicted with the network model,it would not be clear which we should question most.For revealed preference data, options considered included(i)self-completion diaries(Mahmassani and Jou,2000),(ii)automatic tracking through GPS(Jan et al.,2000;Quiroga et al.,2000;Taylor et al.,2000),and(iii)licence plate surveys(Schaefer,1988).Regarding self-comple-tion surveys,from our own interview experiments with self-completion questionnaires it was evident that travellersfind it relatively difficult to recall and describe complex choice options such as a route through an urban network,giving the po-tential for significant errors to be introduced.The automatic tracking option was believed to be the most attractive in this respect,in its potential to accurately map a given individual’s journey,but the negative side would be the potential sample size,as we would need to purchase/hire and distribute the devices;even with a large budget,it is not straightforward to identify in advance the target users,nor to guarantee their cooperation.Licence plate surveys,it was believed,offered the potential for compromise between sample size and data resolution: while we could not track routes to the same resolution as GPS,by judicious location of surveyors we had the opportunity to track vehicles across more than one location,thus providing route-like information.With time-stamped licence plates, the matched data would also provide journey time information.The negative side of this approach is the well-known poten-tial for significant recording errors if large sample rates are required.Our aim was to avoid this by recording only partial licence plates,and employing statistical methods to remove the impact of‘spurious matches’,i.e.where two different vehi-cles with the same partial licence plate occur at different locations.Moreover,extensive simulation experiments(Watling,1994)had previously shown that these latter statistical methods were effective in recovering the underlying movements and travel times,even if only a relatively small part of the licence plate were recorded,in spite of giving a large potential for spurious matching.We believed that such an approach reduced the opportunity for recorder error to such a level to suggest that a100%sample rate of vehicles passing may be feasible.This was tested in a pilot study conducted by the project team,with dictaphones used to record a100%sample of time-stamped, partial licence plates.Independent,duplicate observers were employed at the same location to compare error rates;the same study was also conducted with full licence plates.The study indicated that100%surveys with dictaphones would be feasible in moderate trafficflow,but only if partial licence plate data were used in order to control observation errors; for higherflow rates or to obtain full number plate data,video surveys should be considered.Other important practical les-sons learned from the pilot included the need for clarity in terms of vehicle types to survey(e.g.whether to include motor-cycles and taxis),and of the phonetic alphabet used by surveyors to avoid transcription ambiguities.Based on the twin considerations above of planned interventions and survey approach,several candidate studies were identified.For a candidate study,detailed design issues involved identifying:likely affected movements and alternative routes(using local knowledge of York CC,together with an existing network model of the city),in order to determine the number and location of survey sites;feasible viewpoints,based on site visits;the timing of surveys,e.g.visibility issues in the dark,winter evening peak period;the peak duration from automatic trafficflow data;and specific survey days,in view of public/school holidays.Our budget led us to survey the majority of licence plate sites manually(partial plates by audio-tape or,in lowflows,pen and paper),with video surveys limited to a small number of high-flow sites.From this combination of techniques,100%sampling rate was feasible at each site.Surveys took place in the morning peak due both to visibility considerations and to minimise conflicts with tourist/special event traffic.From automatic traffic count data it was decided to survey the period7:45–9:15as the main morning peak period.This design process led to the identification of two studies:2.1.Lendal Bridge study(Fig.1)Lendal Bridge,a critical part of York’s inner ring road,was scheduled to be closed for maintenance from September2000 for a duration of several weeks.To avoid school holidays,the‘before’surveys were scheduled for June and early September.It was decided to focus on investigating a significant southwest-to-northeast movement of traffic,the river providing a natural barrier which suggested surveying the six river crossing points(C,J,H,K,L,M in Fig.1).In total,13locations were identified for survey,in an attempt to capture traffic on both sides of the river as well as a crossing.2.2.Fishergate study(Fig.2)The partial closure(capacity reduction)of the street known as Fishergate,again part of York’s inner ring road,was scheduled for July2001to allow repairs to a collapsed sewer.Survey locations were chosen in order to intercept clockwiseFig.1.Intervention and survey locations for Lendal Bridge study.around the inner ring road,this being the direction of the partial closure.A particular aim wasFulford Road(site E in Fig.2),the main radial affected,with F and K monitoring local diversion I,J to capture wider-area diversion.studies,the plan was to survey the selected locations in the morning peak over a period of approximately covering the three periods before,during and after the intervention,with the days selected so holidays or special events.Fig.2.Intervention and survey locations for Fishergate study.In the Lendal Bridge study,while the‘before’surveys proceeded as planned,the bridge’s actualfirst day of closure on Sep-tember11th2000also marked the beginning of the UK fuel protests(BBC,2000a;Lyons and Chaterjee,2002).Trafficflows were considerably affected by the scarcity of fuel,with congestion extremely low in thefirst week of closure,to the extent that any changes could not be attributed to the bridge closure;neither had our design anticipated how to survey the impacts of the fuel shortages.We thus re-arranged our surveys to monitor more closely the planned re-opening of the bridge.Unfor-tunately these surveys were hampered by a second unanticipated event,namely the wettest autumn in the UK for270years and the highest level offlooding in York since records began(BBC,2000b).Theflooding closed much of the centre of York to road traffic,including our study area,as the roads were impassable,and therefore we abandoned the planned‘after’surveys. As a result of these events,the useable data we had(not affected by the fuel protests orflooding)consisted offive‘before’days and one‘during’day.In the Fishergate study,fortunately no extreme events occurred,allowing six‘before’and seven‘during’days to be sur-veyed,together with one additional day in the‘during’period when the works were temporarily removed.However,the works over-ran into the long summer school holidays,when it is well-known that there is a substantial seasonal effect of much lowerflows and congestion levels.We did not believe it possible to meaningfully isolate the impact of the link fully re-opening while controlling for such an effect,and so our plans for‘after re-opening’surveys were abandoned.3.Estimation of vehicle movements and travel timesThe data resulting from the surveys described in Section2is in the form of(for each day and each study)a set of time-stamped,partial licence plates,observed at a number of locations across the network.Since the data include only partial plates,they cannot simply be matched across observation points to yield reliable estimates of vehicle movements,since there is ambiguity in whether the same partial plate observed at different locations was truly caused by the same vehicle. Indeed,since the observed system is‘open’—in the sense that not all points of entry,exit,generation and attraction are mon-itored—the question is not just which of several potential matches to accept,but also whether there is any match at all.That is to say,an apparent match between data at two observation points could be caused by two separate vehicles that passed no other observation point.Thefirst stage of analysis therefore applied a series of specially-designed statistical techniques to reconstruct the vehicle movements and point-to-point travel time distributions from the observed data,allowing for all such ambiguities in the data.Although the detailed derivations of each method are not given here,since they may be found in the references provided,it is necessary to understand some of the characteristics of each method in order to interpret the results subsequently provided.Furthermore,since some of the basic techniques required modification relative to the published descriptions,then in order to explain these adaptations it is necessary to understand some of the theoretical basis.3.1.Graphical method for estimating point-to-point travel time distributionsThe preliminary technique applied to each data set was the graphical method described in Watling and Maher(1988).This method is derived for analysing partial registration plate data for unidirectional movement between a pair of observation stations(referred to as an‘origin’and a‘destination’).Thus in the data study here,it must be independently applied to given pairs of observation stations,without regard for the interdependencies between observation station pairs.On the other hand, it makes no assumption that the system is‘closed’;there may be vehicles that pass the origin that do not pass the destina-tion,and vice versa.While limited in considering only two-point surveys,the attraction of the graphical technique is that it is a non-parametric method,with no assumptions made about the arrival time distributions at the observation points(they may be non-uniform in particular),and no assumptions made about the journey time probability density.It is therefore very suitable as afirst means of investigative analysis for such data.The method begins by forming all pairs of possible matches in the data,of which some will be genuine matches(the pair of observations were due to a single vehicle)and the remainder spurious matches.Thus, for example,if there are three origin observations and two destination observations of a particular partial registration num-ber,then six possible matches may be formed,of which clearly no more than two can be genuine(and possibly only one or zero are genuine).A scatter plot may then be drawn for each possible match of the observation time at the origin versus that at the destination.The characteristic pattern of such a plot is as that shown in Fig.4a,with a dense‘line’of points(which will primarily be the genuine matches)superimposed upon a scatter of points over the whole region(which will primarily be the spurious matches).If we were to assume uniform arrival rates at the observation stations,then the spurious matches would be uniformly distributed over this plot;however,we shall avoid making such a restrictive assumption.The method begins by making a coarse estimate of the total number of genuine matches across the whole of this plot.As part of this analysis we then assume knowledge of,for any randomly selected vehicle,the probabilities:h k¼Prðvehicle is of the k th type of partial registration plateÞðk¼1;2;...;mÞwhereX m k¼1h k¼1172 D.Watling et al./Transportation Research Part A46(2012)167–189。

Location and Mapping for AUV Using aForward-Looking SonarTiedong Zhang, Yongjie Pang, Dapeng Jiang, Wenjing ZengThe Key Laboratory of Underwater V ehicle, Harbin Engineering University, Harbin 150001(E-mail: zhangtiedong@)Abstract—a method of location and mapping based on the forward-looking sonar is proposed to resolve the problem of the singleness of u nderwater navigation and the low precision of location. At first, the featu res of working environment models were represented based on the theory of Symmetries and Pertu rbation model, and the information of environment was obtained by a forward-looking sonar, and the sonar data were corrected according to DVL data. Secondly, Hou gh transform was u sed to extract geometric features in the corrected data of environment, and the current position of AUV was located by the matching of geometric features with features of working environment models. At last, the simulation results show the viability of the proposed approach.Keywords—Location and Mapping, AUV, forward-looking sonar෎Ѣࠡᢊໄ൷Ո∈ϟᴎ఼Ҏᅮԡᷛ೒ᓴ⎅ᷟᑲ∌ᵄྰ໻؃᳒᭛☝જᇨⒼᎹ࣏໻ᄺ∈ϟᴎ఼Ҏₑ⚍ᅲɀᅸજᇨⒼ 150001ᨬᡅ⍌ᇍ∈ϟᇐხ᠟↉ऩϔˈᅮԡϡޚܲՈ⒲Lˈᦤߎњϔ࢑෎Ѣࠡᢊໄ൷೒ڣՈ∈ϟᴎ఼Ҏᅮԡᷛ೒ᮍ⊩DŽŊܜ₋ϬѠමᇍࢴᨘࡼ 63 ῵ൟˈᓎএњ∈ϟᴎ఼Ҏ԰Ϯɳ๗῵ൟĽᕕDŽᔧ∈ϟᴎ఼ҎԡѢ԰Ϯɳ๗ᯊˈỞẋࠡᢊ೒ڣໄ൷ቻᕫɳ๗ֵᙃˈᑊḍ᥂໮᱂ࢦợᑺᩥ᭄᥂ẟᜐ᭄᥂᷵ℷDŽ݊⃵ˈỞẋ₋Ϭ+RXJKবᤶᮍ⊩ᦤপߎɳ๗᭄᥂ЁՈ޴ԩĽᕕDŽᇚ޴ԩĽᕕϢɳ๗ЁՈ޴ԩֵᙃּऍ‑ˈܲᅮᔧࠡᯊࠏ∈ϟᴎ఼ҎՈԡาDŽ᳔ৢˈ᭛Ёඝߎњӓף᪙ɀ඗ᵰDŽ᪙ɀ඗ᵰᜬᯢ᠔ᦤߎՈᮍ⊩ᰃ᳝ᬜՈDŽ݇⏲᪑ᅮԡᷛ೒ˈ∈ϟᴎ఼Ҏˈࠡᢊໄ൷1ˊᓩᣄᴎ఼Ҏᇍ਼ೈɳ๗Ոᛳکϔָᰃࢿࡼᴎ఼Ҏۘऺ:ඳЁՈϔϾۘऺᮍ৥ˈফࠄۘऺҎਬՈᑓ⊯݇⊼ˈᑊᦤߎњ᩼໮ᓔ߯ᗻՈਜ਼⊩ˈϔѯᏆඓ៤ࡳᑨϬ೼┊Ϟᴎ఼ҎЁDŽ໐ᇍѢ∈ϟᴎ఼Ҏ(AUV)ᴹ᪸ˈϵѢফࠄ᥶⌟Ӵᛳ఼└ࠊˈϔָ≵᳝ᕫࠄ⏅ܹՈۘऺDŽ಴ℸֲࠡAUV໮᭄₋ϬỞẋ໮᱂ࢦẟᜐხԡ᥼ਜ਼Ոᮍ⊩ᴹᛳکԡาֵᙃˈԚϵѢ᥼ਜ਼ẋ࣏Ёᄬ೼௳ࢳ᪳Ꮒˈ᠔ቻᕫՈԡาֵᙃᰃϡஂܲՈˈ◄ᡅAUVḍ᥂∈☦GPSᇐხֵᙃẟᜐԡา᷵ℷˈ᠔ҹẝ࢑ᛳکᮍ⊩ᕜ▂⒵ᱷAUVᅲ┉ՈᎹ԰◄ᡅDŽẕ޴ᑈˈ╓؄೒ڣໄ൷ᡔᴃՈ៤cˈՓᕫ∈ϟ᥶⌟᠟↉᳈ࡴЄᆠˈ߽Ϭໄ೒ڣẟᜐɳ๗᥶⌟г៤Ўৃ࿁DŽ⍌ᇍẝ࢑ᚙމˈᴀ᭛ᦤߎњ෎Ѣࠡᢊໄ൷೒ڣՈAUVɳ๗ᛳکਜ਼⊩DŽ᭛Ёҹऩ⊶ᴳࠡᢊໄ൷԰ЎЏᡅӴᛳ఼ˈ໮᱂ࢦ԰ЎṉࡽӴᛳ఼ˈ₋ϬѠමᇍࢴᨘࡼ(SP)ᮍ⊩ˈᓎএњ޴ԩܗ௤῵ൟ੠ɳ๗ഄ೒῵ൟˈᇚҢໄ൷᭄᥂ЁᦤপՈָඃĽᕕϢɳ๗᭄᥂ЁՈ޴ԩĽᕕẟᜐऍ‑ˈܲᅮߎAUV᠔໘ՈԡาDŽ2ˊ೒ڣ᷵ℷЎߚᵤ⒲Lᮍ֓ˈᓎএ໻ഄതᷛிOXYZ੠ẔࡼതᷛிZ YXO′′′ˈབ೒ ᠔߾DŽᔧAUV☝ℶᯊˈϸതᷛிതᷛḸₑড়ˈໄ൷෎⓹ԡѢYḸˈ⊩ඃ⊓XḸᮍ৥ˈᇍѢĭԧ⚍),,(βαρMˈρЎᲡࡿₓˈαЎ∈ᑇᮍԡᢖˈβЎؒᢖˈ݊Ҫখ᭄৪োՈọপখᢅ᭛Dz[3]DŽ᪂ໄ⊶೼∈ϟՈᑇഛӴ᪁ợᑺsmV/15002=ˈໄ൷ϵথễໄ⊶ࠄ᥹ফ⚍Mডᇘໄ⊶ՈᯊⒸЎ1tˈ᳝߭˖2/12tV=ρ (1)ϵѢࠡᢊໄ൷೼ൖָᮍ৥Ϟ≵᳝ߚṬɋˈ಴ℸҙ࿁ඝߎM⚍೼∈ᑇ☦ԡาֵᙃˈ໐ϡ࿁ඝߎM⚍ൖָᮍ৥ϞProceedings of the 7thWorld Congress on Intelligent Control and Automation June 25 - 27, 2008, Chongqing, China೒ +RXJKָඃᦤপD ໄ൷ᠿᦣܼ᱃E ᦤপָඃF ᡩࠬाⒸՈԡาֵᙃDŽḍ᥂ໄ⊶Ӵ᪁ॳˊˈM ⚍೼Ѡමᑇ☦ՈԡาֵᙃϢৠϔໄ⊶⊶☦੠∈ᑇ☦ϞՈѸ⚍M ′ᰃϔႸՈ> @> @ˈेM ⚍ᇍᑨ೼∈ᑇ☦ՈതᷛᑨЎ¯®­⋅=⋅=αραρsin cos y x (2) ϵℸৃᢅˈĭԧ⚍мবɴᬥҙϢαρ,ϸϾখ᭄᳝݇DŽAUV ೼∈ϟẔࡼᯊ݋᳝݁ϾႮϵᑺˈߚ߿ᰃ⊓؄Z Y X ,,ḸՈᑇࢿẔࡼ੠඙؄Z Y X ,,ḸՈᮟḰẔࡼˈԚϵѢAUV ೼∈ϟẔࡼᯊἹᕾᑇ☦Ẕࡼ؛᪂ˈ಴ℸৃҹᩨЎϞẴẔࡼᰃּѦưএՈˈৃߚ߿໇ᔕߚᵤDŽŊܜ໇ᔕᑇࢿẔࡼˈህֲࠡՈۘऺᴹ׏ˈAUV ᳔໻ხợs m V /4max ≤ˈϵ(1)ᓣˈAUV ᳔໻ԡࢿS Ϣ᥶⌟Სࡿ↨ؐЎ˖1500422max==V V l S (3) े˖ 005.015008≈=l S (4) ϵϞᓣৃᢅˈּᇍѢ᥶⌟Სࡿᴹ᪸ˈᑇࢿẔࡼ᠔ᓩ᰻ՈأᏂৃҹᗑЩϡᩥˈ໇ᔕࠄໄ൷ᤶ࿁఼ᬊথᯊⒸᕜڱˈ಴ℸৃҹ؛᪂AUV ܜᑇࢿẔࡼݡথ᥹ໄ൷᭄᥂DŽ෎ѢϞ☦؛᪂ˈ߭↣⃵থᇘ੠᥹ᬊẋ࣏῁াᄬ೼ᮟḰẔࡼᕅડˈेρ೼↣Ͼথᬊẋ࣏ЁᰃֱᣕϡবՈˈ໐া᳝αᰃব࣪ՈDŽᇍѢḰࡼẔࡼᚙމˈϵѢAUV ໮ЎಲḰԧˈ಴ℸ῾ᨛẔࡼৃҹᗑЩϡᩥDŽᇍѢ൹ؒẔࡼᴹ᪸ˈໄ൷෎⓹ԡา≵᳝থϣᬍবˈাᰃ׃ӄᢖᑺথϣব࣪ˈ಴ℸাᰃಲ⊶᭄᥂Ոᔎᑺ᳝᠔ব࣪ˈ໐ϡӮᕅડᲡࡿؐρ੠αˈᬙℸᚙމϡӮᕅડໄ൷᭄᥂៤ڣDŽᇍѢأხՈᚙމˈ෎⓹തᷛிӮথϣأḰˈे⚍M ′ԡาϵ),(αρবࠄ),(αρ′′ˈབ೒ ᠔߾DŽҸy θᜬ߾ხأᢖˈϵ޴ԩ݇ிˈৃ᥼ᇐߎϸϾԡาПⒸ݇ிЎ˖¯®­′=+′=ρρθααy(5) ϵϞẴߚᵤৃᢅˈࠡᢊໄ൷мবᚙމՈߎɴˈЏᡅᰃϵѢAUV ՈأხỤ៤ՈDŽḍ᥂݀ᓣ(5)ᇍмব∴ຕẟᜐ᷵ℷˈ᷵ℷ඗ᵰབ೒33ˊĽᕕᦤপḍ᥂Hough বᤶՈ෎ᴀᗱᛇ> @ˈᇚָᢖതᷛிЁ೒ڣाⒸЁՈ↣ϔ⚍),(i i y x ḍ᥂ᓣ(6))sin()cos(θθρy x += (6) ᯴ᇘࠄHough ाⒸЁՈϔඈ௳ࡴ఼),(θρH ˈ⒵ᱷ(6)ᓣՈ↣ϔ⚍ˈᇚՓᇍᑨՈ᠔᳝௳ࡴ఼ЁՈؐࡴ DŽབᵰ೒ڣЁࣙ৿ϔᴵָඃˈ᳝߭ϔϾᇍᑨՈ௳ࡴ఼Ӯߎɴሔᾬ᳔໻ؐˈỞẋẔ⌟Hough ाⒸЁՈሔᾬ᳔໻ؐˈ ेৃܲᅮϢ᪩ᴵָඃᇍᑨՈϔᇍখ᭄),(θρˈᅲɴָඃՈẔ⌟DŽϵѢໄ൷᭄᥂ҹᵕതᷛ),(θρᔶᓣẟᜐᄬټˈ಴ℸ◄ᡅᇚ(1)ᓣḰ࣪ࠄᵕതᷛிϟẟᜐ໘ˊDŽḍ᥂ָᢖതᷛி੠ᵕതᷛிՈḰᤶ݀ᓣ˖¯®­==θρθρsin cos y x ੠(6)ᓣৃᕫᵕതᷛிЁ↣ϔ⚍),(i i θρᇍيবᤶЎ˖)sin()sin()cos()cos(θθρθθρρi i i i += (7) े˖)cos(i i θθρρ−= (8)ЎޣᇥHough খ᭄Ͼ᭄ˈ඗ড়ऩ⊶ᴳໄ൷ՈᠿᦣĽ⚍ˈࠊᅮབϟᢈ߭˖ᡩࠬⓌؐDŽᔧ᥹ᬊࠄᶤᣛ৥ϞՈ₋ḋ᭄᥂ˈ߭Ϣᡩ ࠬⓌؐẟᜐ↨ṇˈҹܲᅮᰃ৺݋᳝ᡩࠬᰈḐDŽᆩ᭄᥂໻ѢᡩࠬⓌؐᯊˈ߭݋᳝ᡩࠬᰈḐˈ᭄ֱؐᣕϡবDŽᆩᇣѢᡩࠬⓌؐˈ߭প⍜ᡩࠬᰈḐˈ₋ḋ᭄ؐ⏙0DŽỞẋ᪂এᡩࠬⓌؐˈ┑Ԣᡩࠬখ᭄ˈޣᇣᩥਜ਼ₓDŽߚࡆⓌؐDŽߚࡆⓌؐϬᴹᇍড়៤Ոܼ᱃೒ڣẟᜐߚࡆˈএ┨ᥝాໄ⚍ˈᦤপߎऎඳDŽᔧ᭄᥂໻ѢߚࡆⓌؐᯊˈ߭᪩⚍Ўऎඳݙ⚍៪ṽො⚍ˈϵ᪩⚍ᓔྟẟᜐऎඳϣ⑃ˈָࠄϣ⑃⚍᭄ؐᇣѢ0ᯊذℶˈ᠔ϣ៤Ո⚍Ўऎඳ⚍៪ṽА⚍ˈ݊ᅗ⚍ᩨЎాໄ⚍ˈ᭄ؐ⏙0DŽ᪙ɀ඗ᵰབ೒4᠔߾DŽ೒ мব᷵ℷ඗ᵰD мবՈ∴ຕໄ೒ڣE ᷵ℷৢՈ∴ຕໄ೒ڣ ೒ ໻ഄതᷛிϢẔࡼതᷛி߾೒ мবॳˊ߾ᛣ೒ 63῵ൟЁതᷛ݇ி4ˊᅮԡᷛ೒4.1SP ޴ԩܗ௤῵ൟབ೒ ೼SP ῵ൟЁˈ↣Ͼ޴ԩܗ௤F ೎༘ϔϾതᷛிE തᷛிՈX ḸϢ޴ԩܗ௤Ո޴ԩḸৠ৥> @DŽ޴ԩܗ௤F ೼ܼሔതᷛிW ЁՈԡาϬ৥ₓWF X ᜬ߾DŽ೼ѠමՈᚙމϟˈϬϸϾതᷛ੠ϔϾḰᢖᴹẟᜐᦣẴˈे˖T WF y x X ),,(φ= (9)ᢆ⌟ࠄՈĽᕕ೼Ꮖکഄ೒ЁՓϬখ໇ிM ᜬ߾ˈ೼ܼሔതᷛிՈԡาЎWM X DŽ೼῵ൟЁˈӴᛳ఼ֵᙃՈϡܲᅮᗻᰃỞẋϔϾඝᅮഛؐ੠णᮍᏂՈʌᮃՁాໄᜬ߾Ոˈे᠔᳝Ոϡܲᅮᗻ῁ᰃ ೼ሔᾬതᷛிϟᜬ߾Ոˈ಴ℸᇚԄᩥՈ᪳ᏂϬϔϾּᇍѢ޴ԩܗ௤খ໇ிՈᖂߚԡา৥ₓF d ᜬ߾DŽᡞF ೼ܼሔതᷛிЁՈԡาԄᩥᩴ԰WFXˆˈ߭F ೼ܼሔതᷛிЁՈԡาৃᜬ߾Ў FWF WF d X X ⊕=ˆ (10) ᓣЁ˖T F yF xF F d d d d ),,(ϕ=DŽ⊕ᰃᇍ܄⃵തᷛẔਜ਼Ո੄ݭᔶᓣ খᢅ᭛Dz[7] ϟৠDŽᅮНᨘࡼ৥ₓF p ҹ⍜┨ᜬ߾ԡาႮϵᑺՈₓˈे˖F FT F p Bd = (11)F B ᰃ޴ԩܗ௤ՈԈ╓ڭ⓹ˈ޴࢑݌ൟՈ޴ԩܗ௤Ոᨘࡼ৥ₓপؐᢅ᭛Dz> @DŽϵᓣ(9)-(11)ৃᢅˈ೼SP ῵ൟЁˈϔϾ޴ԩܗ௤Ոԡาֵᙃৃᜬ߾Ў˖),,ˆ,ˆ(F F F WFWF B C p X L = (12) WF L ࢴЎ޴ԩܗ௤F ೼ܼሔതᷛிЁՈϡܲᅮԡาˈᓣЁ৘৥ₓՈ݇ிЎ˖])ˆ)(ˆ[(][ˆˆT F F F F F F F FT F WF WFp p pp E C p E pp B X X −−==⊕= (13)݊ЁF C ᜬ߾णᮍᏂڭ⓹ˈϬᴹᜬ߾źᗕԄᩥՈϡܲᅮᗻDŽ4.2 SP ഄ೒῵ൟ೼ഄ೒῵ൟЁˈЏᡅࣙᣀ޴ԩܗ௤῵ൟ੠ᴎ఼Ҏ῵ൟˈᇍѢ޴ԩܗ௤῵ൟৃᢅᓣ(12)ˈᇍѢᴎ఼ҎR ῵ൟˈӓ+ᓣ(12) ৃݭЎ˖()RR R WR WRB C d xL ,,ˆ,ˆ= (14) ϵѢR ೼ɳ๗Ёϡᄬ೼ᇍࢴˈ߭˖3I B R =໇ᔕᴎ఼ҎՈϡܲᅮԡา੠N Ͼഄ೒ĽᕕՈϡܲᅮԡาҹঞᅗӀПⒸՈּ݇ᗻˈϬ৥ₓW xˆᜬ߾ᴎ఼Ҏ੠ഄ೒ՈԡาԄᩥˈ߭ϵᓣ(12) (14)ৃکSP ഄ೒ՈԡาԄᩥ੠ᨘࡼ৥ₓЎ˖»»»»»»¼º««««««¬ª=»»»»»»¼º««««««¬ª=N N F F F R W WF WF WF WR W p p p d p x x x x x##2121ˆˆˆˆˆ (15)ेSP ഄ೒῵ൟՈᷛޚᔶᓣৃݭЎ˖()W W W W map B C p xSP ,,ˆ,ˆ= (16) ᓣЁ৘খ᭄ĭˊ৿НৠϞDŽЎߚᵤᮍ֓ˈᇚSP ഄ೒Ոᨘࡼ৥ₓW p ੠णᮍᏂڭ⓹W C ᩴЎ˖¸¸¹·¨¨©§=»¼º«¬ª=F RF RF RW F R W C CC C C p d p ;4.3 ĽᕕԄᩥ4.3.1źᗕᮍ࣏ᓎএ᪂k ᯊࠏᴎ఼Ҏԡา1|−k k R W x Ў˖k k k k k k R R R W R W x x x 11|11|−−−−⊕= (17)಴Ўᴎ఼Ҏ೼k ᯊࠏּᇍѢ1−k ᯊࠏՈԡࢿЎ˖()3,,ˆ,ˆ1111I C d x L kk k k k k kk R R R R R R R R −−−−= (18)᠔ҹּᇍԡา৥ₓk k R R x 1−Ў˖k k k k k k R R R R R R d x x 111ˆˆ−−−⊕= (19)ᓣЁ˖k k R R x 1ˆ−ᜬ߾њᴎ఼ҎՈԡࢿԄᩥˈ᭛ЁỞẋჽԡ᥼ਜ਼ቻᕫDŽ ()k k kk R R R R C N d11,0~ˆ−−ˈᜬ߾њჽԡ᥼ਜ਼ՈஂᑺDŽᇚᓣ(19)ҷܹ(17)ᓣˈᭈˊᕫ˖k k k k k k k k R R R R WR WR d x x x 111|11|ˆ−−−−−⊕⊕= (20) ᭈˊᕫ˖1|1|1|ˆ−−−⊕=k k k k k k R WR WR d xx (21)ϵ(21)ᓣৃᢅˈᔧᴎ఼ҎࢿࡼৢˈSP ഄ೒೼k ᯊࠏՈԡา৥ₓৃᜬ߾Ў˖»»»»»»»»¼º««««««««¬ª⊕=»»»»»»»»¼º««««««««¬ª=−−−−−−−−−−−−−11,11,211,1|1|111,1,21,11ˆˆˆˆ1k k N k k k k k k k k k k N k k k k k k WF WF WF R WR WF WF WF WR W k k x x x d x xx x x x ##(22) ϵѢᴎ఼ҎخԢợẔࡼˈ߭ৃҹᩨЎ᪳Ꮒṇᇣˈ߭ 1|−k k R d ੄࣪Ў˖k k k k k k k k R R R R R R d d J d 11|111|−−−−−+≈ (23)ᓣЁ:1−k k R R J ᜬ߾▉ৃ↨ᜐ߫ᓣ[9]ிඣՈźᗕᮍ࣏k k k k w x F x +=+1 (24) ৃᬍݭЎ˖k k R R k W k k W kd G pF p 11−+=− (25)ᓣЁ˖¸¸¹·¨¨©§=¸¸¹·¨¨©§=×××××−113330 , 001N k N N N N R R k I G I J F k k SP ഄ೒Ոᨘࡼ৥ₓ೼k ᯊࠏՈϔℹԄᩥW k k p 1|ˆ−੠णᮍᏂڭ⓹Wk k C 1|−ߚ߿Ў˖»»¼º««¬ª=»»¼º««¬ª=−−−−−−−−1|11|111|1|ˆˆˆˆˆ1|k k k k k k k k k k M R R R M R Wk k p d J p d pϵवᇨ᳐Ⓒ⊶݀ᓣˈৃݭߎSP ഄ೒೼k ᯊࠏՈᨘࡼ৥ₓ੠णᮍᏂڭ⓹ҹঞवᇨ᳐๲֎བϟ˖11111)( ,)( , ˆˆ−−−−−+=−=−=T kk k T k W k k k T kWk k k Wk k k k W k k W k k W k G S G H CH H CK C H K I C h K p pᇚ⌟ₓֵᙃϢᏆکഄ೒ֵᙃּᙑড়ˈ࿁໳ᇍSP ഄ೒Ոᨘࡼ৥ₓ੠णᮍᏂڭ⓹ẟᜐₑᮄԄᩥˈᦤʌஂᑺDŽSP ഄ೒ЁՈ৘Ͼܗ௤ˈࣙᣀᴎ఼Ҏԡา੠ഄ೒ĽᕕՈԡาԄᩥˈ῁ᝯᙑড়ՈᮄֵᙃₑᮄԄᩥˈ಴ЎᅗӀ῁ᝯࣙ৿೼ৠϔϾźᗕ৥ₓЁˈ໐ϨĽᕕⒸՈּ݇ᗻгϵѢிඣणᮍᏂڭ⓹Ոₑᮄᩥਜ਼໐᳈ᮄˈ໐ᔧAUV ẟܹϔϾᄺдẋՈɳ๗Ёˈɳ๗ЁՈ᠔᳝ՈĽᕕՈϡܲᅮᗻ῁ᇚ┑ԢˈࣙᣀὧѯAUV ≵᳝ᢆ⌟ࠄԚᰃϢᴎ఼Ҏᢆ⌟ࠄՈĽᕕּ݇ՈĽᕕDŽ⌟ₓᮍ࣏᪂),,1(,ˆ,n m y m k ""∈ᰃᇍźᗕ৥ₓՈϔඈ⌟ₓˈm k y ,ᰃᢆ⌟Ոˊ᩾ؐˈ߭˖ m k m k m k u y y,,.,ˆ+= ),0(~,.m k m k S N u (26)ᇍᓣ(26)ՓϬϔϾ╔ᓣߑ᭄ᴹᜬ߾> @ˈे˖ 0),(,,=m k k m k y x f (27)ϵѢ(27)ᓣЁՈźᗕ৥ₓ৿᳝ᮍԡ-ˈ಴ℸᰃϔϾ☢ඃᗻՈᮍ࣏DŽᇍ݊ඃᗻ࣪ˈ೼)ˆ,ˆ(,m k k y x໘⋄ࢦሩᓔˈᗑЩʌ⓺ˈাপϔ⓺ˈৃᕫ˖0)ˆ()ˆ(),(,,,,,,,=−+−+≅m k m k m k k k m k m k m k k m k yy G xx H h y x f (28)ᓣЁ˖)ˆ,ˆ(,,,)ˆ,ˆ(,,,,,,, ,)ˆ,ˆ(m k k m k k y xmk m k m k y xkm k m k m k k m k m k y f G x f H y x f h ∂∂=∂∂==প˖Tmk m k m k m k m k m k m k k m k m k m k G S G R u G v xH h z ,,,,,,,,,, ,ˆ=−=+−=߭(28)ᓣৃᭈˊЎ˖),0(~ ,,,,,,m k m k m k k m k m k R N v v x H z += (29)Ľᕕऍ‑೼Ľᕕऍ‑Ёˈ₋ϬϹḐ൪ᴳᴵӊˈᇏᡒᢆ⌟Ľᕕ੠ഄ೒ĽᕕՈऍ‑DŽ✊ৢৠᯊ໇ᔕ೼ഄ೒ĽᕕЁ᳝ऍ‑Ո᠔᳝Ոᢆ⌟Ľᕕˈᇚ᠔᳝ᢆ⌟Ľᕕᔧ԰ϔϾ⏋ড়Ľᕕ໘ˊ> @DŽ ᇍѢ↣ϔϾ݋᳝ऍ‑ՈĽᕕˈඃᗻ࣪ৢՈ⌟ₓᮍ࣏Ոி᭄Ո᪪ඊ᥼ᇐৃᢅ᭛Dz> @ˈᇚ᠔᳝Ոᢆ⌟ݭ៤ڭ⓹Ոᔶᓣ᳝߭˖;21»»»»»¼º«««««¬ª=N E E E p p p y # ¸¸¸¸¸¹·¨¨¨¨¨©§=NN N N N E T E ET EE E E E T E E E E E E E C C C C C C C C C S "#%##""2121211211߭(28)ᓣবЎ˖()()()0ˆˆ,=−+−+≈y y G pp H h y p f W W W (30) ඃᗻ࣪ৢˈᕫࠄՈ৘Ͼখ᭄Ў˖»»»»¼º««««¬ª=N h h h h #21 ¸¸¸¸¸¹·¨¨¨¨¨©§=N E E E G G G G "#%##""00000021¸¸¸¸¸¹·¨¨¨¨¨©§=0000000000000002211""#########""""NNE R E RE R H H H H H H HᮄĽᕕ⏏ࡴᴎ఼Ҏ೼k ᯊࠏՈᢆ⌟Ёˈৃ࿁ࣙ৿᳝Ꮖکഄ೒Ё≵᳝ऍ‑ՈĽᕕˈỞẋᴎ఼ҎՈԡาԄᩥ੠ᢆ⌟Ľᕕּᇍᴎ఼ҎՈԡาԄᩥˈָ᥹ᇚ᳾ऍ‑Ոᢆ⌟ĽᕕࡴܹࠄSP ഄ೒DŽՓϬখ໇ிE ′ᜬ߾ϔϾ೼ഄ೒Ё≵᳝ऍ‑Ոᢆ⌟ĽᕕˈᅗּᇍѢᴎ఼ҎՈԡาЎ˖()E E E E R E R B C p x L ′′′′′=,,ˆ,ˆ (31) ϵѢℸᯊᏆکᴎ఼Ҏ೼ܼሔതᷛிՈԡาԄᩥˈ᠔ҹᢆ⌟Ľᕕ೼ܼሔതᷛிՈԡาЎ˖W E E W E R WR E W d x x x x ′′′′⊕=⊕=ˆ (32)ϵѢ᪳Ꮒṇᇣˈᬙ˖E T E R ER E T E R ER W E p B d J p B d J d +≈⊕= (33) ߭ᮄĽᕕՈᨘࡼ৥ₓৃݭЎ˖ E R ER E W E p d J B p += (34) ϵᓣ(33)(34)ৃکˈᇚᮄՈĽᕕࡴܹഄ೒ˈּᔧѢ೼SP ഄ೒Ոźᗕ৥ₓЁࡴܹϔ-ˈे˖()»»»¼º«««¬ª⊕=′RE WR WM WR W x x x xx ˆˆˆˆˆ (35) ᨘࡼ৥ₓবЎ˖()»»»¼º«««¬ª+=»»»¼º«««¬ª=′E R ER E M R W E M R W p d J B d d p d d p (36) ◄ᡅᣛߎՈᰃˈӏԩᮄࡴܹՈĽᕕ῁া࿁┑ԢிඣՈஂᑺˈফࠄᕜ໻ాໄᑆᡄՈֵᙃˈབᵰࡴܹՈ῵ൟЁˈϞႷ᳝ৃ࿁ᓩ᰻ிඣՈϡࣷᅮˈ಴ℸˈᇍѢফࠄ໻ాໄᑆᡄՈ᭄᥂ˈᡅঞᮽഄࠨ┨DŽ5. ӓף᪙ɀ඗ᵰϵѢ∈∴ɳ๗੄ऩˈᑆᡄᇣˈ಴ℸᓎњњּᇍṇ໡ᴖՈӓףഎ᱃ˈᑊࡴܹњ⍋⌕Ոᑆᡄˈབ೒6᠔߾DŽ೒Ёതᷛऩԡ₋Ϭỿṕऩԡˈ೼ഄ೒Ёọ߭ Ͼ޴ԩĽᕕˈᷛোߚ߿Ў ǃ ǃ ǃ ǃ ǃ ǃ ˈབ೒ ᠔߾DŽ೼೒Ё᠔߾ɳ๗Ёˈ᪂าϡৠՈܹ∈⚍ˈ$89ߚ߿ḍ᥂ Ͼ޴ԩĽᕕˈऍ‑ߎᔧࠡ᠔໘Ոԡาˈ໐$89ףᅲԡาỞẋ⌟ᲡӴᛳ఼ቻᕫˈϸ້↨ṇ඗ᵰᢅ೒ DŽҢ೒8(a)-8(e)Ёৃҹ׏ߎˈऍ‑ᅮԡৢˈᴎ఼ҎՈԡาഛᮍᏂϢӴᛳ఼Ո⌟ᲡഛᮍᏂ೼ৠϔϾ᭄ₓ൫ϞˈẂࠄњᳳᳯՈஂᑺDŽ໐೼೒8(f)Ёˈᴎ఼ҎՈԡา೼X ᮍ৥Ϟߎɴњṇ໻Ո᪳ᏂˈẝЏᡅᰃϵѢᴎ఼Ҏ೼ẝϾԡาϞᢆ⌟ࠄՈĽᕕЁˈা᳝ϔᴵָඃᦤկњ೼X ᮍ৥ϞՈ൪ᴳỤ៤ˈབ೒7᠔߾DŽᔧẔࡼᯊˈAUV Ởẋձ᥂Ľᕕऍ‑඗ᵰᴹׂℷხԡ᥼ਜ਼᭄᥂ᴹᅲɴᅮԡˈ᪙ɀ඗ᵰབ೒9੠೒10DŽ೒Ё᠔߾AUV ೼ߚ߿ᓔਃϢ݇⒱ऍ‑ׂℷՈᚙމϟˈ໮⃵ඓẋϔϾɳᔶᲳᕘˈᕫࠄՈϸඈᅲɀ඗ᵰDŽҢᅲɀ඗ᵰᴹ׏ˈỞẋऍ‑ׂℷৢՈხԡ᥼ਜ਼Ո඗ᵰ⊶ࡼṇᇣˈAUV ხẽࣷᅮˈপᕫᕜདՈᅮԡᬜᵰDŽ᳔ৢˈ೼ӓףЁᇍඝᅮϡᅠᭈഄ೒ՈᚙމϟˈAUV ႮЏᷛ೒࿁࡯ẟᜐɀ᪅ˈ᪙ɀ඗ᵰབ೒13᠔߾DŽ೒11ᜬ߾8ܜᓎএՈAUV ԰Ϯɳ๗ഄ೒ˈ໐೒12߭ᰃAUV ףᅲ԰Ϯɳ๗ഄ೒ˈϢ8کഄ೒ּ↨ˈɳ๗ֵᙃথϣњব࣪DŽ೒ӓףഄ೒ ೒ ܹ∈⚍ ໘ᢆ⌟Ľᕕ೒ ᥼ਜ਼Ϣऍ‑೒ D ᮍԡᢖ᪳Ꮒ ೒ E Სࡿ᪳Ꮒ೒ $89ẔࡼḬẽ೒(a) ხԡ᥼ਜ਼Ḭẽ೒ E ׂℷḬẽ ೒ ϡৠܹ∈⚍Ոᅮԡ඗ᔧAUV ೼ɳ๗Ё⓿␌ϔ਼ৢˈ֕⌟ࠄᮄՈɳ๗޴ԩĽᕕˈᑊᇚ݊ࡴܹࠄഄ೒῵ൟЁˈᕫࠄׂℷৢՈഄ೒བ೒13DŽ೼᪙ɀЁˈAUV ᇚϸᴵᢆ⌟ࠄ໐೼ഄ೒Ё≵᳝ᷛᩴՈָඃᷛ⊼њߎᴹˈẝₐᇚṇڱՈָඃᷛᩴԡLine1ˈṇ⑃ՈָඃᷛᩴЎLine2ˈᇚᷛᩴߎՈָඃϢףᅲɳ๗ЁՈָඃ↨ṇˈৃҹᕫࠄϟ☦Ո඗ᵰ˖Line126.032.391.138=−==−==−=mapping real mapping real mapping real y y y x x x ααδαδδ Line287.268.691.32=−==−==−=mapping real mapping real mapping real y y y x x x ααδαδδ Ң඗ᵰЁৃᢅˈṇЎޚܲഄݡɴњɳ๗ЁՈ޴ԩĽᕕDŽ6.඗᩾᪙ɀᜬᯢˈᴀ᭛ᦤߎՈAUV ₋Ϭࠡᢊໄ൷Ոᅮԡऍ‑ᮍ⊩ˈ೼∈ϟԢాໄɳ๗Ёˈ࿁໳ቻᕫṇʌՈᅮԡஂᑺˈ೼᳝ాໄᑆᡄՈᚙމϟˈৃṇདഄᇍხԡ᥼ਜ਼Ո඗ᵰẟᜐׂℷˈᇍѢAUV ∈ϟᇐხᰃৃᜐՈDŽখ໇᭛Dz[1] Sonar Log Data & Decode Manual. Tritec. International Limited.[2] Zhang Tiedong. “The post-processing of the image for theforward looking sonar.” M. S. dissertatio n, Dept. Ship. Eng., Harbin Eng. Univ., Harbin, China, 2004 [3] ᮑϣẂ ┰ᄋ᪡൹ᗻ ࣫Ҁ ೑⓶ᎹϮߎČࠂ 1995.pp.10-11 [4] Lei Wanming, Hu Xuecheng. “High Resolution Airborne SARMotion Compensation B ased on Echo Data.” The Jo urnal o f Electro nics & Info rmatio n Techno lo gy . vol. 26, no. 12, pp. 1908-1914, 2004[5] Li Yanping, Xing Meng dao, Bao Zheng. “Motioncompensation method based on radar returns.” The Journal of Data acquisition & Processing . vol. 22, no. 1, pp.1-7, 2007. [6] Hough ˈP V C .A method and means for recognizing complexpatterns ˈU. S Pattern , 3069654,1962[7] Juan D. Tarods. “Representing partial and uncertain sensorialinformation using the theory of symmetries,” in IEEE Int . R o b o tics and Aut o mati o n C onf ., Nice, France, pp.1799-1084,1992.[8] H.F Durrant-Whyte. Integration Coordination and Controlof Multi-Sensor Robot Systems . Kluwer Academic Pub., Massachusetts, 1988.[9] R. Smith, M. Self, and P. Cheeseman, “Estimating uncertainspatial relationships in robotics,” in Autonomous robot vehicles . Springer Verlag New York, Inc., pp. 167–193, 1990. [10] S. Roumeliotis and J. Burdick, “Stochastic cloning: ageneralized framework for processing relative statemeasurements,” in Proceedings .ICRA ’02. IEEE International Co nference o n Ro bo tics and Auto matio n , vol. 2, pp. 1788–1795, 2002.[11] S. Thrun, Y. Liu, D. Koller, A. Y. Ng, Z. Ghahramani, and H. Durrant-Whyte, “Simultaneous Localization and Mapping with Sparse Extended Information Filters,” The International Jo urnal o f Ro bo tics Research , vol. 23, no. 7-8, pp. 693–716, 2004.೒ᷛ೒඗ᵰ೒ ॳ᳝ഄ೒೒ ᅲ┉ഄ೒。

一种基于区域局部搜索的NSGA II 算法栗三一 1王延峰 1乔俊飞 2黄金花3摘 要 针对局部搜索类非支配排序遗传算法 (Nondominated sorting genetic algorithms, NSGA II)计算量大的问题,提出一种基于区域局部搜索的NSGA II 算法(NSGA II based on regional local search, NSGA II-RLS). 首先对当前所有种群进行非支配排序, 根据排序结果获得交界点和稀疏点, 将其定义为交界区域和稀疏区域中心; 其次, 围绕交界点和稀疏点进行局部搜索. 在局部搜索过程中, 同时采用极限优化策略和随机搜索策略以提高解的质量和收敛速度, 并设计自适应参数动态调节局部搜索范围. 通过ZDT 和DTLZ 系列基准函数对NSGA II-RLS 算法进行验证, 并将结果与其他局部搜索类算法进行对比, 实验结果表明NSGA II-RLS 算法在较短时间内收敛速度和解的质量方面均优于所对比算法.关键词 非支配排序遗传算法, 分区搜索, 局部搜索, 多目标优化引用格式 栗三一, 王延峰, 乔俊飞, 黄金花. 一种基于区域局部搜索的NSGA II 算法. 自动化学报, 2020, 46(12): 2617−2627DOI 10.16383/j.aas.c180583A Regional Local Search Strategy for NSGA II AlgorithmLI San-Yi 1 WANG Yan-Feng 1 QIAO Jun-Fei 2 HUANG Jin-Hua 3Abstract In order to reduce the amount of calculation and keep the advantage of local search strategy simultan-eously, this paper proposed a kind of nondominated sorting genetic algorithms (NSGA II) algorithm based on re-gional local search (NSGA II-RLS). Firstly, get corner points and sparse point according to the results of non-dom-inated sorting of current populations, define those points as the centers of border areas and sparse area respectively;secondly, search around the corner points and sparse point locally during every genetic process; NSGA II-RLS ad-opts extreme optimization strategy and random search strategy simultaneously to improve the quality of solutions and convergence rate; adaptive parameter is designed to adjust local search scope dynamically. ZDT and DTLZ functions are used to test the effectiveness of NSGA II-RLS, the performance is compared with four other reported local search algorithms. Results show that: the solution quality of NSGA II-RLS is better than the other methods within limited time; the time complexity of NSGA II-RLS needed to achieve the set IGD value is less than the oth-er methods.Key words Nondominated sorting genetic algorithms (NSGA II), regional search, local search, multi-objective op-timizationCitation Li San-Yi, Wang Yan-Feng, Qiao Jun-Fei, Huang Jin-Hua. A regional local search strategy for NSGA II algorithm. Acta Automatica Sinica , 2020, 46(12): 2617−2627带精英策略的非支配排序遗传算法(Nondom-inated sorting genetic algorithms, NSGA II)作为一种启发式算法, 通过模拟进化论的自然选择和遗传学机理, 可以在不考虑实际工程内部工作方式的情况下求解高度复杂的非线性最优值问题, 被广泛应用于经济结构优化[1]、路径规划[2]、生产调度[3]等实际工程中. 然而作为一种类随机搜索算法, NSGA II 存在收敛速度慢的问题[4].针对NSGA II 收敛速度慢的问题, 已有的研究表明局部搜索策略可以有效提高种群收敛速度, 并且在靠近Pareto 前沿时避免陷入局部极优[5]. 目前已经提出的局部搜索算法可以分为两类: 随机搜索算法和定向搜索算法.X )Y )Y X Y X,Y 随机搜索算法将指定解(设为 周围区域作为搜索区间, 对该解增加一较小的随机值形成新的解(设为 , 若 支配 则 取代 之后以 为中心继续进行搜索. 一些研究者认为初始种群对局部搜索算法的效果有重要影响[6−7], 初始种群分布范收稿日期 2018-09-01 录用日期 2019-01-18Manuscript received September 1, 2018; accepted January 18,2019全国教育科学规划一般课题(BJA170096), 湖北省教育科学规划课题 (2018GB148), 教育部新一代信息技术创新项目(2019ITA04002),河南省科技攻关项目基金 (202102310284)资助Supported by General Project of National Education Science (BJA170096), Education Science Project of Hubei Province (2018GB148), Innovation Project in Information Technology of Education Ministry (2019ITA04002), and Key Projects of Sci-ence and Technology of Henan Province (202102310284)本文责任编委 魏庆来Recommended by Associate Editor WEI Qing-Lai1. 郑州轻工业学院 郑州 4500022. 北京工业大学信息学部 北京 1001243. 武汉船舶职业技术学院 武汉 4300001. Zhengzhou University of Light Industry, Zhengzhou 4500022. Faculty of Information Technology, Beijing University of Technology, Beijing 1001243. Wuhan Institute of Shipbuild-ing Technology, Wuhan 430000第 46 卷 第 12 期自 动 化 学 报Vol. 46, No. 122020 年 12 月ACTA AUTOMATICA SINICADecember, 2020围越大、分布越均匀, 随机搜索的效果越好, 从这一方面出发设计了基于任务分解的种群初始化方法,产生初始解之后使用随机搜索尝试从不同的初始解逼近Pareto前沿. 部分研究者尝试将单目标局部搜索算法扩展应用于解决多目标优化问题[8]. 还有一些专家尝试调整搜索区域和搜索范围以提高局部搜索的效率[9−10]. 在专家学者的努力下, 已有的随机搜索类算法可以有效提高种群的收敛速度, 避免种群陷入局部极优, 然而在每一次迭代过程中都需要对每个解进行局部搜索, 普遍存在计算复杂度高的问题.定向搜索算法通过梯度或任务分解等方法指定搜索方向, 使初始种群朝着指定方向收敛. 一些专家使用梯度求导等方法获得搜索方向[11−12], 可以有效指导种群向Pareto前沿逼近, 但是求导计算量太大. 为了避免求导, 研究者利用目标空间几何信息[13]、解的邻域信息[14]和父代与子代之间差别信息[15]等获得搜索方向. 定向搜索算法通过指定搜索方向,搜索效率高, 但由于搜索方向固定, 对初始种群的分布特性要求很高, 且方向函数的计算也增加了计算复杂度, 与随机算法相同的是随机算法在每一次迭代过程中对每个解进行局部搜索, 计算成本高.随机搜索算法和定向搜索算法在搜索过程中对每个解均进行局部搜索, 计算复杂度很高[9, 16], 限制了局部搜索算法在对优化速度要求较高场合的应用. 针对这一问题本文提出基于区域局部搜索的NSGA II算法(NSGA II based on regional local search, NSGA II-RLS). NSGA II-RLS以NSGA II 为框架, 在交叉变异操作过程后进行局部搜索, 首先根据非支配排序结果获得交界点(目标空间中单个目标向量方向上值最大的解)和稀疏点(除了交界点以外拥挤距离最大的点), 将其定义为交界区域和稀疏区域中心, 然后围绕交界点和稀疏点进行局部搜索. 局部解由极限优化变异策略[17]和随机搜索策略产生, 在局部搜索过程中设计自适应参数动态调整局部搜索范围, 提高了局部搜索的效率. NSGA II-RLS主要有以下优势:1)交界点和稀疏点可以直接根据非支配排序结果获得, 不需要计算密度和梯度, 计算量小.2)在交界区域和稀疏区域进行局部搜索可以同时提高收敛速度和增加种群分布的均匀性.3)只在交界区域和稀疏区域进行局部搜索可以避免计算资源浪费, 有效降低计算复杂度.4)自适应参数的设定可以使算法在初期具有较大的搜索范围, 靠近Pareto前沿时具有较小的搜索范围, 提高了局部搜索的搜索效率.通过基准多目标优化实验验证算法的有效性.实验结果证明NSGA II-RLS可以有效提高NSGA II 算法的收敛速度. NSGA II-RLS在有限时间内得到的解的质量明显优于其他算法; 评价指标达到设定值所消耗的计算量明显少于其他算法; 优化效果优于固定搜索范围的局部搜索方法.1 非支配、拥挤距离排序P,N,X i iS i C i,S i X i C iX i C i=0P1,C i C i=0P2,NSGA II算法结合非支配关系与拥挤距离对非支配解进行排序. 假设当前种群为 种群规模为 对每个个体(第个个个体), 设两个参数和 为被支配的个体的集合, 为支配的个体的数量. 将所有的个体组成集合 其为第1级非支配解, 也是当前的非支配解集.然后将剩下的所有个体的减1, 此时的个体组成集合 其为第2级非支配解, 重复上述过程直到所有解完成分级.iD i,i然后对每一级的解进行拥挤距离排序, 设第个解的拥挤距离为 定义为在个体周围包含个体本身但不包含其他个体的最小的长方形(周长最小的长方形), 如图1所示. 将同一级别的解按照拥挤距离从大到小排列, 最边界的解的拥挤度为无穷大, 同一级别的解拥挤距离越大越好.i图 1 个体的拥挤距离iFig. 1 Crowded distance of individual2 NSGA II-RLS算法为了解决目前局部搜索多目标优化算法计算复杂度高的问题, 本文提出基于一种基于分区局部搜索的多目标优化算法NSGA II-RLS. NSGA II-RLS根据NSGA II算法的非支配排序结果直接获2618自 动 化 学 报46 卷得交界点和稀疏点, 将其设为交界区域和稀疏区域中心, 不需要额外的计算量; 在交界点和稀疏点周围进行局部搜索, 在提高种群收敛速度的同时保证了进化过程中种群的多样性和均匀性; 分区搜索较以往的局部搜索算法计算复杂度明显降低; 为了提高局部搜索的效率, 对局部搜索范围进行自适应动态调整.2.1 获取区域中心局部搜索可以有效提高收敛速度[8−12], 然而对每一个解都进行局部搜索计算量很大, 而且对远离Pareto 前沿的被支配解进行局部搜索对种群逼近Pareto 前沿贡献不大. 因此人们尝试分区域进行搜索, 区域搜索的中心思想是在重点区域进行局部搜索, 有侧重点的搜索以增加搜索效率, 如在目标空间进行聚类得到聚类中心、通过拐点确定中心等,然后围绕中心点进行局部搜索.然而通过聚类、密度计算等确定中心, 任意两个解之间的欧氏距离都需要计算, 增加了算法的计算量. 因此本文直接利用非支配排序和拥挤距离排序获得搜索区域中心, 非支配排序和拥挤距离排序是NSGA II 算法的固有步骤, 因此利用其获得搜索区域中心不会增加计算复杂度.t P,m,P P 1P 1m m +1m +1P 1m m +1m 获取区域中心具体方法如下: 设第 代种群为 目标函数个数为 根据第1节对种群 进行非支配排序和拥挤距离排序, 则第1级非支配解 为当前种群的非支配解. 由拥挤距离排序机理可知,交界点的个数等于目标函数个数, 在种群 中前 个个体的拥挤距离为无穷大, 第 个个体的拥挤距离为除了边界点之外最大的点, 意味着第 个个体周围解的密度最低. 由此可知非支配解集 中前 个个体为交界点, 第 个个体为稀疏点, 对应着 个交界区域的中心和稀疏区域的中心.本文以交界点和稀疏点为中心进行局部搜索,交界点是至少在一个目标向量方向上的极大值或极小值, 以交界点为中心进行搜索是确保种群分布范围的一个措施. 拥挤距离最低的点周围是解最稀疏的区域, 围绕稀疏点进行局部搜索可以增加种群分布的均匀性. 因此本文以交界点和拥挤距离最小的点为中心进行局部搜索是有意义的.以两目标优化问题为例, 图2显示了NSGA II-RLS 算法的种群进化过程. 从图中可以看到以下三点: 局部搜索可以提高种群收敛速度; 围绕稀疏点进行局部搜索可以增加种群的分布均匀性; 围绕交界点进行局部搜索可以保持种群的分布范围.2.2 局部搜索搜索区域中心确定之后需要围绕区域中心进行局部搜索, 局部搜索的方法对局部搜索效果有重要影响. 本文采用文献[18]提出的局部搜索方法, 同时采用极限优化策略和随机优化策略产生局部解,可以平衡局部搜索算法的全局搜索能力和局部搜索能力.X =(x 1,x 2,···,x n ),n n,极限优化的具体方法如下[17]: 设被搜索区域的中心解为 为决策变量个数,种群规模为N , 则产生局部解个数为 变异公式为X i =(x 1,···,x ′i ,···,x n ),0<i ≤n(1)x ′i =x i +α×βmax (x i ),0<i ≤n(2)图 2 NSGA II-RLS 算法的种群进化过程Fig. 2 The evolution process of NSGA II-RLS algorithm12 期栗三一等: 一种基于区域局部搜索的NSGA II 算法2619βmax (x i )=max [x i −l i ,u i −x i ],0<i ≤n(4)x i h q q l i u i i βmax (x i )x i 其中, 为决策变量; 为0到1之间的随机数; 为正实数, 称为形状参数, 根据文献[17]将 设为11; 和 分别为第 个决策变量的下界和上界; 为当前决策变量 可变动的最大值.X =(x 1,x 2,···,x n ),n,n N,随机局部搜索策略为: 设当前搜索中心解为 为决策变量个数, 种群总数为 随机搜索产生局部解个数为种群总数的20%[19] (如不能整除则取整), 随机搜索变异公式为γrand (γ)(−γ,γ)0.1N (m +1)(n +⌈0.3N ⌉)m 其中, 为搜索范围参数, 用于确定搜索范围的大小; 表示取值在 之间的随机数. 同时还产生 个随机解以保证种群的多样性[17]. 以上变异策略共产生 个局部解, 为目标函数个数.2.3 自适应参数设定T max 对于随机局部搜索策略, 在算法运行初期较大的搜索范围可以提高全局搜索能力, 使种群快速靠近Pareto 前沿. 在靠近Pareto 前沿后, 较小的搜索范围可以提高种群逼近Pareto 前沿的能力. 本文根据最大优化时间( )设计参数动态调整方法.T max ,根据经验与已有的研究成果[9−10, 19]设定的取值范围为(0.05,0.2), 设最大优化时间为 搜索范围调整公式为t T max T max t T max T max γT max 其中, 为从优化开始到当前时刻的时间, 使用式(9)时需 的值大于5, 若小于5则单位量级可以降低一个等级, 如当 为3 s, 则可设其为3 000 ms, 与 统一单位. 图2和图3分别显示了 为5 s 和1 000 s 时 的变化曲线, 从图中可以看出, 其变化轨迹完全相同, 且其变化规律符合预期的搜索范围变化规律, 因此本文提出的自适应公式针对不同的 设定值具有很好的自适应性.γ注1. 的取值范围本文根据经验与相关文献[9–10, 19]进行设定, 并没有通过实验验证, 后期研究可以通过实验找到最佳的搜索范围.注2. 多数实际优化问题的Pareto 前沿是未知的, 本文提出的搜索范围调整方法并不需要知道Pareto 前沿, 而且需要设定的参数很少, 适用于解决实际问题.2.4 NSGA II -RLS 算法流程N ;T max ;n ;u =(u 1,u 2,···,u n )l =(l 1,l 2,···,l n );m ;q ;p c pm ;ηc ηm .根据实际MOP 问题设定算法参数: 初始种群数量为 最大优化时间 决策向量维数 决策变量取值上界 和下界 目标函数个数 形状参数 交叉概率和变异概率 交叉参数 和变异参数 NSGA II-RLS 的具体算法流程如下:P I ={X 1,X 2,···,X N }X i =(x ′1,x ′2,···,x ′n )i =1,2,···,n 步骤1. 在取值范围内随机初始化种群 , 其中, , .P I P C ,步骤2. 对 进行非支配、拥挤距离排序, 当前种群中所有非支配解记为 根据排序结果确定交界区域中心和稀疏区域中心.P I P M .步骤3. 按照标准NSGA II 算法对 中的种群进行交叉变异, 形成子代T max γ图 3 为5 s 时 的变化曲线γT max Fig. 3 the change curve of when is set to be 5 sT max γ图 4 为1 000 s 时 的变化曲线γT max Fig. 4 the change curve of when isset to be 1 000 s2620自 动 化 学 报46 卷(m +1)(n +⌈0.3N ⌉)P N .步骤4. 围绕交界中心和稀疏区域中心进行局部搜索, 共产生 个局部解, 将其集合设为种群 P I P M P N N P O ,P I =P O .步骤5. 将 、 和 合并, 并对所有解进行非支配、拥挤距离排序, 从中选取最优的 个解形成下一代种群 并设 T max 步骤6. 重复步骤2~5, 当达到最大优化时间 或目标精度时进行步骤7.P I 步骤7. 当前 中的非支配解即为得到的最优解.与其他局部搜索NSGA II 算法相比, NSGA II-RLS 算法只在边界和稀疏区域进行局部搜索, 计算量明显降低, 并且能够保证算法的收敛速度和种群的分布性; 采用两种局部搜索策略, 使算法在搜索前期和后期都有较好的搜索效率; 搜索范围的自适应调整平衡了全局搜索与局部搜索的比重; 本文提出的搜索范围调整方法不需要提前获知Pareto 前沿, 需要设定的参数也较少, 符合实际工程应用要求.下面分析NSGA II-RLS 算法运行一代的时间复杂度(函数计算次数, 每一次公式的计算都增加计算复杂度1). NSGA II-RLS 的时间开销主要集中在子代目标函数值求解部分.0.5N +(m +1)(n +⌈0.3N ⌉)n m N,O (mN )O (0.8N +0.3mN )子代目标函数值求解: 交叉变异和局部搜索共产生 个解, 因决策变量个数 和目标个数 一般远小于 因此该步骤计算复杂度为 , 一代函数计算次数为 .O (mN )O (N 2)O (N )N 综合以上分析, NSGA II-RLS 的计算复杂度为 . 当单个解的局部搜索解数量与本文相同时, 全部解都进行局部搜索的随机搜索算法的计算复杂度为 . 对于定向搜索算法, 由于需要通过梯度等方法计算方向, 计算复杂度一般高于随机搜索算法. 标准NSGA II 的一代计算复杂度为 ,由此可以看出, 本文提出算法的计算复杂度与标准NSGA II 都为 的一次方级别, 且其系数也不大,因此NSGA II-RLS 的计算复杂度远远低于其他局部搜索算法, 而且其搜索范围也可以自适应调整.注3. NSGA II-RLS 算法采用分区搜索机制,重点在交界区和解稀疏区域局部搜索, 因此, 对于Pareto 非连续的问题, 在优化过程中可能出现某些片段无解的情况, 所以NSGA II-RLS 适用于解决Pareto 前沿连续的问题.3 仿真实验本文通过双目标ZDT 系列函数和三目标DTLZ系列函数对算法NSGA II-RLS 进行验证, 测试函数的特征及参数如表1所示. 实验结果与基于密度O (0.8N +0.1N 2)O (0.5N +C 2nN )O (0.5N +15(n −3)N )O (2.5N +C 2n N )O (N 2)O (n 2N )O (nN )O (n 2N )O (0.8N +0.3mN )O (mN )的局部搜索算法NSGA II-DLS [18]、随机局部搜索算法ED-LS [9]、变深度随机局部搜索算法MOILS [10]和定向搜索算法DirLS [15]进行对比, 一代函数计算次数分别为 、 、 和 , 计算复杂度分别为 、 、 和 . 本文提出的NS-GA II-RLS 的一代函数计算次数为 ,计算复杂度为 . 采用MATLAB10.0b 软件进行仿真实验, 处理器为3.60 GHz, 8.00 GB 内存,Microsoft 实验环境.N n ηc ηm 1/n ;q 双目标实验初始种群规模 为100, 三目标实验初始种群规模为200; 决策变量个数 如表1所示; 所有实验均采用模拟二进制交叉变异方法, 交叉参数 和变异参数 均设为20, 交叉和变异概率分别为0.9和 形状参数 设为11.I 采用综合评价指标 IGD (Inverted generation-al distance) 和种群多样性指标 进行评价. IGD 的计算公式为P ∗P d (v,P )v P |P ∗|其中, 为帕累托解集; 为近似解; 为向量 与解集 中的向量的欧氏距离最小值; 为帕累托解集中解的个数.I 指标 的计算公式为d f d l d i ,i =1,2,···,N −1d m d i 其中, 和 为帕累托末端解和所得解集边界解之间的欧氏距离, 为所获得的连续非支配解之间的欧氏距离, 为所有 的平均值.实验结果使用秩和检验(Wilcoxon ranksum test), 在0.05显著性水平上说明不同结果之间的差异显著性.注4. 在使用优化算法做基准实验时, 大部分文表 1 测试函数参数Table 1 Paramter setting of the test functions函数Pareto 前沿特征决策变量维度目标维度种群规模ZDT1凸3021 000ZDT2凹3021 000ZDT4凸1021 000ZDT6凹1021 000DTLZ1非凸非凹732 500DTLZ2凹734 096DTLZ3凹734 096DTLZ4凹1234 09612 期栗三一等: 一种基于区域局部搜索的NSGA II 算法2621章将种群规模设定为100或200[3−7], ZDT 系列问题是双目标问题, 相对比较简单, 因此我们对ZDT 实验设定种群规模为100, DTLZ 系列为三目标问题,比较复杂, 我们设定种群规模为200. Deb 等[20]在提出NSGA II 算法时设定交叉变异参数为20, 交叉概率为0.9, 变异概率为1/n , 之后人们在使用和改进NSGA II 算法时, 有一些保持了该参数设定[21−23],而这一部分不是本文研究的重点, 因此本文也按照文献[20]对交叉变异参数进行设定. 文献[17]提出了极限优化算法, 本文根据文献[17]将形状参数q 设为11.3.1 实验1T max 本实验设定当ZDT 系列函数IGD 值达到0.01、三目标函数IGD 值达到0.1 (所对比算法双目标函数IGD 最优值在0.01以下但接近0.01, 三目标函数IGD 最优值在0.1以下并接近0.1)时停止优化, 对双目标和三目标函数最大优化时间 分别设为50 s 和200 s (各对比算法达到目标精度的时间多数情况下小于此设定值, 为了防止未达到目标精度而因达到最大优化时间停止优化, 将最大优化时间设定为较大的值比较合理), 计算停止时的总函数计算次数. 本实验用于验证算法在达到目标精度时函数计算总次数, 可以反映算法的收敛速度,实验结果如表2和表3所示.从表2和表3可以看出, 本文提出的NSGA II-RLS 在所有实验中达到目标精度时总函数计算次数远远低于其他对比算法. 在双目标ZDT 系列实验中, NSGA II-RLS 的总进化代数是最低的. 对于ZDT1和ZDT2, MOILS 的函数计算次数均低于ED-LS, 说明当n 较大时对搜索范围的调整可以有效提高搜索效率. NSGA II-RLS 的实验结果优于ED-LS, 说明极限优化策略的引入可以有效提高种群的收敛效果. DirLS 的函数计算次数低于ED-LS,且该两种算法的局部解生成机制同为2-opt (2-op-timization)方法, 说明在双目标实验中指向性参数的加入对收敛性有较大提升. 同时从表2中可以看到, NSGA II-RLS 结果的波动范围最小, 说明NSGA II-RLS 有较好的稳定性.在三目标DTLZ 系列实验中, MOILS 在DTLZ1、DTLZ3和DTLZ4中的进化代数最低, 除NSGA II-DLS 之外的算法总进化代数差异不超过30%, 因此总函数计算次数的差异主要由单步函数计算次数产生. 从总函数计算次数可以反映出NSGA II-RLS 单步计算复杂度低的优势. MOILS 在三目标实验中总函数计算次数是最高的, 结合双目标实验结果表 2 ZDT 系列函数IGD 值达到0.01时对比算法的总时间复杂度与进化代数 (连续10次实验求平均)Table 2 For ZDT series function the comparison of total time complexity and the evolution algebra of differentalgorithms when IGD value reaches 0.01 (mean value of ten consecutive experimental results)算法ZDT1ZDT2ZDT3ZDT4复杂度代数复杂度代数复杂度代数复杂度代数NSGA II-RLS 2 10015 2 38017 1 96014 1 40010NSGA II-DLS [17]46 440 (+)4334 560 (+)3231 320 (+)2932 400 (+)30ED-LS [9]2 569 450 (+)591 698 450 (+)39168 350 (+)37163 800 (+)36MOILS [10]1 419 250 (+)351 135 400 (+)28327 050 (+)31232 100 (+)22DirLS [15]1 312 500 (+)301 137 500 (+)26156 750 (+)33114 000 (+)24注: (+) 表示NSGA II-RLS 的结果明显优于相应的算法表 3 DTLZ 系列函数IGD 值达到0.1时对比算法的总时间复杂度与进化代数 (连续10次实验求平均)Table 3 For DTLZ series function the comparison of total time complexity and the evolution algebra of differentalgorithms when IGD value reaches 0.1 (mean value of ten consecutive experimental results)算法DTLZ1DTLZ2DTLZ3DTLZ4复杂度代数复杂度代数复杂度代数复杂度代数NSGA II-RLS 29 9208817 3401933 6609927 54041NSGA II-DLS [17]378 000 (+)17599 360 (+)46416 880 (+)193192 240 (+)89ED-LS [9]369 800 (+)86116 100 (+)27460 100(+)107545 300 (+)41MOILS [10]883 300 (+)73254 100 (+)211 028 500 (+)851 084 000 (+)40DirLS [15]385 400 (+)82159 800 (+)34451 200 (+)96671 300 (+)49注: (+) 表示NSGA II-RLS 的结果明显优于相应的算法2622自 动 化 学 报46 卷n n 可知, 当决策变量维数较高时, 变深度局部搜索的总函数计算次数优于2-opt 方法, 当决策变量维数低时2-opt 方法的总函数计算次数优于变深度搜索. NSGA II-DLS 在所有实验中消耗的进化代数最大, 原因是其只在稀疏解周围进行局部搜索, 导致算法先收敛到前沿附近的某一区域, 之后还需要继续对周边区域进行探索. 由以上分析可知, NSGA II-RLS 产生局部解时使用极限优化策略和随机搜索策略可以有效提高种群收敛速度, 其单步函数计算次数远远低于所对比算法. 从局部解生成公式(1)~(6)可知, 其生成解数量仅与种群规模N 相关,而目前大部分局部搜索方法产生局部解时采用2-opt 及其衍生方法, 产生局部解的个数不仅与N 相关, 而且与决策变量维数 成正相关, 实际问题中, 的维数一般较高, 这也是导致一般局部搜索计算复杂度高的原因之一.从以上实验结果分析可知, 达到目标精度时NSGA II-RLS 的总函数计算次数最低, 说明NSGA II-RLS 具有快速收敛的优点, 并且优化效果比较稳定.注5. 文献[14, 18, 24]以优化函数计算次数作为评价标准, 但对于局部搜索类算法而言, 遗传过程中梯度计算、密度计算等消耗的计算资源也不容忽视. 为了更公平地比较算法效果, 本文以总时间复杂度为评价标准.3.2 实验2T max T max T max T max T max 实际工程中经常以时间作为停止标准, 本实验设定当达到最大优化时间 时停止实验. 本文提出的NSGA II-RLS 优势在于计算复杂度低、收敛速度快, 因此为了体现算法的优势, 本实验对ZDT 系列实验设 为20 s, DTLZ1和DTLZ3实验设 为200 s, DTLZ2实验设 为40 s, DTLZ4实验设 为80 s, 比较实验停止时的IGD 值、GD 值和I 值. 所有实验除局部搜索步骤不同, 其他部分完全相同, 以保证实验对比的公平性. 本实验的目的是检验在有限的时间内各算法的优化效果.实验结果如表4和表5所示.T max 从表4可以看出, 除DTLZ2实验外NSGA II-RLS 的IGD 均值与标准差均为最小. DTLZ2实验各算法实验结果相差不超过30%, 说明对于DT-LZ2函数将 设为40 s 各算法有相对充足时间收敛到Pareto 前沿, 也说明当时间足够长时NSGA II-RLS 的效果反而不如已有的局部搜索算法, 但其优化效果与对比算法中最优的结果差距在5% 以内,从实验1的结果可以看出, 在达到目标精度时NSGA II-RLS 的总计算量不超过对比算法总计算量的10%,也就是说用10%的计算资源达到了95%的优化效果. 在其他实验中NSGA II-RLS 的IGD 均值与方差均最小, 对应实验1的结果可知, 在优化时间有限且较短时, 计算复杂度低、收敛速度快的NSGA II-RLS 算法具有明显优势.ED-LS 与DirLS 算法类似, DirLS 比ED-LS 多加了方向性指标指导种群进化, 在ZDT 函数实验中DirLS 的结果优于ED-LS, 而在三目标实验中ED-LS 的结果反而优于DirLS, 这是由于DirLS 随着种群进化, 根据进化前后两代种群更新方向指标, 在目标个数较少时方向指标可以加快收敛到Pareto 前沿的速度, 但当目标空间维数较高时, 虽然其也能加快收敛速度, 但解的分布性反而不如ED-LS,从而导致其综合评价指标IGD 值不理想. NSGA II-RLS 的实验结果均优于NSGA II-DLS, 说明本文提出的分区域搜索比只在稀疏区域搜索有更好的效果, 计算复杂度更低.从表5可以看出, 对于DTLZ1、DTLZ3和DT-LZ4, NSGA II-RLS 的I 值均值与标准差最低, 说明在较短时间内NSGA II-RLS 获得的种群有较好的分布特性. 对于优化时间充足的DTLZ2, NSGA II-RLS 的I 值较最优的ED-LS 的I 值降低了2.4%,说明当优化时间充足时NSGA II-RLS 的解的分布特性虽不如最优的算法, 但差距不大. 这与从表4得到的结论相一致.从表4和表5的实验结果可以看出, 对于ZDT 系列实验、DTLZ1、DTLZ3和DTLZ4, NSGA II-RLS 的实验结果优于对比算法且具有显著差异性,说明其结果具有统计学意义.由以上分析可知, 对于多目标优化问题, 在有限的较短时间内NSGA II-RLS 获得种群的IGD 值与I 值优于所对比算法, 说明在较短时间内NSGA II-RLS 具有较好的逼近性和分布特性. 从DTLZ2的实验结果可知, 当优化时间充足时NSGA II-RLS 算法效果与最优算法相比虽有所降低, 但差距在5%以内.注6. NSGA II-RLS 的提出主要为了解决局部搜索算法计算量大的问题, 因此将停止时间设定较低可以反映算法的优势. 当算法运行时间足够长时NSGA II-RLS 失去其优势, 但从DTLZ2的实验结果可以看出, 优化时间充足时NSGA II-RLS 的优化效果与其他优秀算法的优化结果差距在5%以内,效果相差不大. 本文提出的算法更适用于对优化快速性要求较高的场合, 时间充足的场合也可以应用.12 期栗三一等: 一种基于区域局部搜索的NSGA II 算法2623。

Discovering Similar Multidimensional TrajectoriesMichail VlachosUC Riverside mvlachos@George KolliosBoston Universitygkollios@Dimitrios GunopulosUC Riversidedg@AbstractWe investigate techniques for analysis and retrieval of object trajectories in a two or three dimensional space. Such kind of data usually contain a great amount of noise, that makes all previously used metrics fail.Therefore,here we formalize non-metric similarity functions based on the Longest Common Subsequence(LCSS),which are very ro-bust to noise and furthermore provide an intuitive notion of similarity between trajectories by giving more weight to the similar portions of the sequences.Stretching of sequences in time is allowed,as well as global translating of the se-quences in space.Efficient approximate algorithms that compute these similarity measures are also provided.We compare these new methods to the widely used Euclidean and Time Warping distance functions(for real and synthetic data)and show the superiority of our approach,especially under the strong presence of noise.We prove a weaker ver-sion of the triangle inequality and employ it in an indexing structure to answer nearest neighbor queries.Finally,we present experimental results that validate the accuracy and efficiency of our approach.1IntroductionIn this paper we investigate the problem of discovering similar trajectories of moving objects.The trajectory of a moving object is typically modeled as a sequence of con-secutive locations in a multidimensional(generally two or three dimensional)Euclidean space.Such data types arise in many applications where the location of a given object is measured repeatedly over time.Examples include features extracted from video clips,animal mobility experiments, sign language recognition,mobile phone usage,multiple at-tribute response curves in drug therapy,and so on.Moreover,the recent advances in mobile computing, sensor and GPS technology have made it possible to collect large amounts of spatiotemporal data and there is increas-ing interest to perform data analysis tasks over this data [4].For example,in mobile computing,users equipped with mobile devices move in space and register their location at different time instants via wireless links to spatiotemporal databases.In environmental information systems,tracking animals and weather conditions is very common and large datasets can be created by storing locations of observed ob-jects over time.Data analysis in such data include deter-mining andfinding objects that moved in a similar way or followed a certain motion pattern.An appropriate and ef-ficient model for defining the similarity for trajectory data will be very important for the quality of the data analysis tasks.1.1Robust distance metrics for trajectoriesIn general these trajectories will be obtained during a tracking procedure,with the aid of various sensors.Here also lies the main obstacle of such data;they may contain a significant amount of outliers or in other words incorrect data measurements(unlike for example,stock data which contain no errors whatsoever).Figure1.Examples of2D trajectories.Two in-stances of video-tracked time-series data representingthe word’athens’.Start&ending contain many out-liers.Athens 1Athens 2Boston 1Boston 2DTWLCSSFigure 2.Hierarchical clustering of 2D series (displayed as 1D for clariry).Left :The presence of many outliers in the beginning and the end of the sequences leads to incorrect clustering.DTW is not robust under noisy conditions.Right :The focusing on the common parts achieves the correct clustering.Our objective is the automatic classification of trajec-tories using Nearest Neighbor Classification.It has been shown that the one nearest neighbor rule has asymptotic er-ror rate that is at most twice the Bayes error rate[12].So,the problem is:given a database of trajectories and a query (not already in the database),we want to find the trajectory that is closest to .We need to define the following:1.A realistic distance function,2.An efficient indexing scheme.Previous approaches to model the similarity between time-series include the use of the Euclidean and the Dy-namic Time Warping (DTW)distance,which however are relatively sensitive to noise.Distance functions that are ro-bust to extremely noisy data will typically violate the trian-gular inequality.These functions achieve this by not consid-ering the most dissimilar parts of the objects.However,they are useful,because they represent an accurate model of the human perception,since when comparing any kind of data (images,trajectories etc),we mostly focus on the portions that are similar and we are willing to pay less attention to regions of great dissimilarity.For this kind of data we need distance functions that can address the following issues:Different Sampling Rates or different speeds.The time-series that we obtain,are not guaranteed to be the outcome of sampling at fixed time intervals.The sensors collecting the data may fail for some period of time,leading to inconsistent sampling rates.Moreover,two time series moving at exactly the similar way,but one moving at twice the speed of the other will result (most probably)to a very large Euclidean distance.Similar motions in different space regions .Objectscan move similarly,but differ in the space they move.This can easily be observed in sign language recogni-tion,if the camera is centered at different positions.If we work in Euclidean space,usually subtracting the average value of the time-series,will move the similar series closer.Outliers.Might be introduced due to anomaly in the sensor collecting the data or can be attributed to hu-man ’failure’(e.g.jerky movement during a track-ing process).In this case the Euclidean distance will completely fail and result to very large distance,even though this difference may be found in only a few points.Different lengths.Euclidean distance deals with time-series of equal length.In the case of different lengths we have to decide whether to truncate the longer series,or pad with zeros the shorter etc.In general its use gets complicated and the distance notion more vague.Efficiency.It has to be adequately expressive but suf-ficiently simple,so as to allow efficient computation of the similarity.To cope with these challenges we use the Longest Com-mon Subsequence (LCSS)model.The LCSS is a varia-tion of the edit distance.The basic idea is to match two sequences by allowing them to stretch,without rearranging the sequence of the elements but allowing some elements to be unmatched .The advantages of the LCSS method are twofold:1)Some elements may be unmatched,where in Eu-clidean and DTW all elements must be matched,even the outliers.2)The LCSS model allows a more efficient approximatecomputation,as will be shown later(whereas in DTW you need to compute some costly Norm).Infigure2we can see the clustering produced by the distance.The sequences represent data collected through a video tracking process.Originally they represent 2d series,but only one dimension is depicted here for clar-ity.The fails to distinguish the two classes of words, due to the great amount of outliers,especially in the begin-ning and in the end of the ing the Euclidean distance we obtain even worse results.The produces the most intuitive clustering as shown in the samefigure. Generally,the Euclidean distance is very sensitive to small variations in the time axis,while the major drawback of the is that it has to pair all elements of the sequences.Therefore,we use the model to define similarity measures for trajectories.Nevertheless,a simple extension of this model into2or more dimensions is not sufficient, because(for example)this model cannot deal with paral-lel movements.Therefore,we extend it in order to address similar problems.So,in our similarity model we consider a set of translations in2or more dimensions and wefind the translation that yields the optimal solution to the problem.The rest of the paper is organized as follows.In section2 we formalize the new similarity functions by extending the model.Section3demonstrates efficient algorithms to compute these functions and section4elaborates on the indexing structure.Section5provides the experimental validation of the accuracy and efficiency of the proposed approach and section6presents the related work.Finally, section7concludes the paper.2Similarity MeasuresIn this section we define similarity models that match the user perception of similar trajectories.First we give some useful definitions and then we proceed by presenting the similarity functions based on the appropriate models.We assume that objects are points that move on the-plane and time is discrete.Let and be two trajectories of moving objects with size and respectively,whereand.For a trajectory,let be the sequence.Definition1Given an integer and a real number,we define the as follows:A Ba and andThe constant controls how far in time we can go in order to match a given point from one trajectory to a point in another trajectory.The constant is the matching threshold(see figure3).Thefirst similarity function is based on the and the idea is to allow time stretching.Then,objects that are close in space at different time instants can be matched if the time instants are also close.Definition2We define the similarity function between two trajectories and,given and,as follows:Definition3Given,and the family of translations,we define the similarity function between two trajectories and,as follows:So the similarity functions and range from to. Therefore we can define the distance function between two trajectories as follows:Definition4Given,and two trajectories and we define the following distance functions:andNote that and are symmetric.is equal to and the transformation that we use in is translation which preserves the symmetric prop-erty.By allowing translations,we can detect similarities be-tween movements that are parallel in space,but not iden-tical.In addition,the model allows stretching and displacement in time,so we can detect similarities in move-ments that happen with different speeds,or at different times.Infigure4we show an example where a trajectory matches another trajectory after a translation is applied. Note that the value of parameters and are also important since they give the distance of the trajectories in space.This can be useful information when we analyze trajectory data.XFigure4.Translation of trajectory.The similarity function is a significant improvement over the,because:i)now we can detect parallel move-ments,ii)the use of normalization does not guarantee that we will get the best match between two u-ally,because of the significant amount of noise,the average value and/or the standard deviation of the time-series,that are being used in the normalization process,can be distorted leading to improper translations.3Efficient Algorithms to Compute the Simi-larity3.1Computing the similarity functionTo compute the similarity functions we have to run a computation for the two sequences.Thecan be computed by a dynamic programming algorithm in time.However we only allow matchings when the difference in the indices is at most,and this allows the use of a faster algorithm.The following lemma has been shown in[5],[11].Lemma1Given two trajectories and,with and,we canfind the in time.If is small,the dynamic programming algorithm is very efficient.However,for some applications may need to be large.For that case,we can speed-up the above computa-tion using random sampling.Given two trajectories and ,we compute two subsets and by sampling each trajectory.Then we use the dynamic programming algo-rithm to compute the on and.We can show that,with high probability,the result of the algorithm over the samples,is a good approximation of the actual value. We describe this technique in detail in[35].3.2Computing the similarity functionWe now consider the more complex similarity function .Here,given two sequences,and constants, we have tofind the translation that maximizes the length of the longest common subsequence of()over all possible translations.Let the length of trajectories and be and re-spectively.Let us also assume that the translationis the translation that,when applied to,gives a longest common subsequence,and it is also the translation that maximizes the length of the longest common subsequence.The key observation is that,although there is an infinite number of translations that we can apply to,each transla-tion results to a longest common subsequence between and,and there is afinite set of possible longest common subsequences.In this section we show that we can efficiently enumerate afinite set of translations,such that this set provably includes a translation that maximizes the length of the longest common subsequence of and .To give a bound on the number of transformations that we have to consider,we look at the projections of the two trajectories on the two axes separately.We define theprojection of a trajectoryto be the sequence of the valueson the -coordinate:.A one di-mensional translation is a function that adds a con-stant to all the elements of a 1-dimensional sequence:.Take the projections of and ,and respec-tively.We can show the following lemma:Lemma 2Given trajectories ,if ,then the length of the longest common subsequence of the one dimensional sequences and,is at least :.Also,.Now,consider and .A translation by ,applied to can be thought of as a linear transformation of the form .Such a transformation will allowto be matched to all for which ,and.It is instructive to view this as a stabbing problem:Con-sider the vertical line segments,where (Figure 5).Bx,i By,i+2Bx,i+3Bx,i+4Bx,i+5Ax,iAx,i+1Ax,i+2fc1(x) = x + c1fc2(x) = x + c2Ax axisBx axisFigure 5.An example of two translations.These line segments are on a two dimensional plane,where on the axis we put elements ofand on the axis we put elements of .For every pair of elementsin and that are within positions from eachother (and therefore can be matched by the algo-rithm if their values are within ),we create a vertical line segment that is centered at the point and extends above and below this point.Since each element in can be matched with at most elements in ,the total number of such line segments is .A translation in one dimension is a function of the form .Therefore,in the plane we de-scribed above,is a line of slope 1.After translatingby ,an element of can be matched to an el-ement of if and only if the line intersects the line segment .Therefore each line of slope 1defines a set of possi-ble matchings between the elements of sequences and.The number of intersected line segments is actually an upper bound on the length of the longest common sub-sequence because the ordering of the elements is ignored.However,two different translations can result to different longest common subsequences only if the respective lines intersect a different set of line segments.For example,the translations and in figure 5intersect different sets of line segments and result to longest common subsequences of different length.The following lemma gives a bound on the number of possible different longest common subsequences by bound-ing the number of possible different sets of line segments that are intersected by lines of slope 1.Lemma 3Given two one dimensional sequences ,,there are lines of slope 1that intersect different sets of line segments.Proof:Let be a line of slope 1.If we move this line slightly to the left or to the right,it still in-tersects the same number of line segments,unless we cross an endpoint of a line segment.In this case,the set of inter-sected line segments increases or decreases by one.There are endpoints.A line of slope 1that sweeps all the endpoints will therefore intersect at most different sets of line segments during the sweep.In addition,we can enumerate the trans-lations that produce different sets of potential matchings byfinding the lines of slope 1that pass through the endpoints.Each such translation corresponds to a line .This set of translations gives all possible matchings for a longest common subsequence of .By applying the same process on we can also find a set of translations that give all matchings of.To find the longest common subsequence of the se-quences we have to consider only thetwo dimensional translations that are created by taking the Cartesian product of the translations on and the trans-lations on .Since running the LCSS algorithm takeswe have shown the following theorem:Theorem 1Given two trajectories and,withand ,we can compute theintime.3.3An Efficient Approximate AlgorithmTheorem 1gives an exact algorithm for computing ,but this algorithm runs in cubic time.In this section we present a much more efficient approximate algorithm.The key in our technique is that we can bound the difference be-tween the sets of line segments that different lines of slope 1intersect,based on how far apart the lines are.Consider again the one dimensional projections. Lets us consider the translations that result to different sets of intersected line segments.Each translation is a line of the form.Let us sort these trans-lations by.For a given translation,let be the set of line segments it intersects.The following lemma shows that neighbor translations in this order intersect similar sets of line segments.Lemma4Let be the different translations for sequences and,where.Then the symmetric difference.We can now prove our main theorem:Theorem2Given two trajectories and,with and,and a constant,we canfind an ap-proximation of the similaritysuch that intime.Proof:Let.We consider the projections of and into the and axes.There exists a translation on only such that is a superset of the matches in the optimal of and.In addition,by the previous lemma,there are translations()that have at most different matchings from the optimal. Therefore,if we use the translations,fortime if we runpairs of translations in the plane.Since there is one that is away from the optimal in each dimension,there is one that is away from the optimal in2dimensions.Setting-th quantiles for each set,pairs of translations.4.Return the highest result.4Indexing Trajectories for Similarity Re-trievalIn this section we show how to use the hierarchical tree of a clustering algorithm in order to efficiently answer near-est neighbor queries in a dataset of trajectories.The distance function is not a metric because it does not obey the triangle inequality.Indeed,it is easy to con-struct examples where we have trajectories and, where.This makes the use of traditional indexing techniques diffi-cult.We can however prove a weaker version of the triangle inequality,which can help us avoid examining a large por-tion of the database objects.First we define:Clearly,or in terms of distance:In order to provide a lower bound we have to maximize the expression.Therefore,for every node of the tree along with the medoid we have to keep the trajectory that maximizes this expression.If the length of the query is smaller than the shortest length of the trajec-tories we are currently considering we use that,otherwise we use the minimum and maximum lengths to obtain an approximate result.4.2Searching the Index tree for Nearest Trajec-toriesWe assume that we search an index tree that contains tra-jectories with minimum length and maximum length .For simplicity we discuss the algorithm for the1-Nearest Neighbor query,where given a query trajectory we try tofind the trajectory in the set that is the most sim-ilar to.The search procedure takes as input a nodein the tree,the query and the distance to the closest tra-jectory found so far.For each of the children,we check if the child is a trajectory or a cluster.In case that it is a trajectory,we just compare its distance to with the current nearest trajectory.If it is a cluster,we check the length of the query and we choose the appropriate value for .Then we compute a lower bound to the distance of the query with any trajectory in the cluster and we compare the result with the distance of the current near-est neighbor.We need to examine this cluster only if is smaller than.In our scheme we use an approximate algorithm to compute the.Consequently,the value offrom the bound we compute for.Note that we don’t need to worry about the other terms since they have a negative sign and the approximation algorithm always underestimates the .5Experimental EvaluationWe implemented the proposed approximation and index-ing techniques as they are described in the previous sec-tions and here we present experimental results evaluating our techniques.We describe the datasets and then we con-tinue by presenting the results.The purpose of our experi-ments is twofold:first,to evaluate the efficiency and accu-racy of the approximation algorithm presented in section3 and second to evaluate the indexing technique that we dis-cussed in the previous section.Our experiments were run on a PC AMD Athlon at1GHz with1GB RAM and60 GB hard disk.5.1Time and Accuracy ExperimentsHere we present the results of some experiments using the approximation algorithm to compute the similarity func-tion.Our dataset here comes from marine mammals’satellite tracking data.It consists of sequences of geo-graphic locations of various marine animals(dolphins,sea lions,whales,etc)tracked over different periods of time, that range from one to three months(SEALS dataset).The length of the trajectories is close to.Examples have been shown infigure1.In table1we show the computed similarity between a pair of sequences in the SEALS dataset.We run the exact and the approximate algorithm for different values of and and we report here some indicative results.is the num-ber of times the approximate algorithm invokes the procedure(that is,the number of translations that we try).As we can see,for and we get very good results.We got similar results for synthetic datasets.Also, in table1we report the running times to compute the simi-larity measure between two trajectories of the same dataset. The approximation algorithm uses again from to differ-ent runs.The running time of the approximation algorithm is much faster even for.As can be observed from the experimental results,the running times of the approximation algorithm is not pro-portional to the number of runs().This is achieved by reusing the results of previous translations and terminat-ing early the execution of the current translation,if it is not going to yield a better result.The main conclusion of the above experiments is that the approximation algorithm can provide a very tractable time vs accuracy trade-off for computing the similarity between two trajectories,when the similarity is defined using the model.5.2Classification using the Approximation Algo-rithmWe compare the clustering performance of our method to the widely used Euclidean and DTW distance functions. Specifically:cover.htmlSimilarityApproximate for K tries Exact9494250.250.3160.18460.2530.00140.0022 0.50.5710.4100.5100.00140.0022 0.250.3870.1960.3060.00180.00281 0.50.6120.4880.5630.00180.00280 0.250.4080.2500.3570.001910.0031 0.50.6530.4400.5840.001930.0031 Table1.Similarity values and running times between two sequences from our SEALS dataset.1.The Euclidean distance is only defined for sequencesof the same length(and the length of our sequences vary considerably).We tried to offer the best possible comparison between every pair of sequences,by slid-ing the shorter of the two trajectories across the longer one and recording their minimum distance.2.For DTW we modified the original algorithm in orderto match both x and y coordinates.In both DTW and Euclidean we normalized the data before computing the distances.Our method does not need any normal-ization,since it computes the necessary translations.3.For LCSS we used a randomized version with andwithout sampling,and for various values of.The time and the correct clusterings represent the average values of15runs of the experiment.This is necessary due to the randomized nature of our approach.5.2.1Determining the values for&The values we used for and are clearly dependent on the application and the dataset.For most datasets we had at our disposal we discovered that setting to more than of the trajectories length did not yield significant improvement.Furthermore,after some point the similarity stabilizes to a certain value.The determination of is appli-cation dependent.In our experiments we used a value equal to the smallest standard deviation between the two trajec-tories that were examined at any time,which yielded good and intuitive results.Nevertheless,when we use the index the value of has to be the same for all pairs of trajectories.5.2.2Experiment1-Video tracking data.The2D time series obtained represent the X and Y position of a human tracking feature(e.g.tip offinger).In conjuc-tion with a”spelling program”the user can”write”various words[19].We used3recordings of5different words.The data correspond to the following words:’athens’,’berlin’,’london’,’boston’,’paris’.The average length of the series is around1100points.The shortest one is834points and the longest one1719points.To determine the efficiency of each method we per-formed hierarchical clustering after computing thepairwise distances for all three distance functions.We eval-uate the total time required by each method,as well as the quality of the clustering,based on our knowledge of whichword each trajectory actually represents.We take all possi-ble pairs of words(in this case pairs)and use the clustering algorithm to partition them into two classes.While at the lower levels of the dendrogram the clustering is subjective,the top level should provide an accurate divi-sion into two classes.We clustered using single,complete and average linkage.Since the best results for every dis-tance function are produced using the complete linkage,we report only the results for this approach(table2).The same experiment is conducted with the rest of the datasets.Exper-iments have been conducted for different sample sizes and values of(as a percentage of the original series length).The results with the Euclidean distance have many clas-sification errors and the DTW has some errors,too.For the LCSS the only real variations in the clustering are for sam-ple sizes.Still the average incorrect clusterings for these cases were constantly less than one().For 15%sampling or more,there were no errors.5.2.3Experiment2-Australian Sign LanguageDataset(ASL).The dataset consists of various parameters(such as the X,Y, Z hand position,azimuth etc)tracked while different writ-ers sign one the95words of the ASL.These series are rel-atively short(50-100points).We used only the X and Y parameters and collected5recordings of the following10 words:’Norway’,’cold’,’crazy’,’eat’,’forget’,’happy’,’innocent’,’later’,’lose’,’spend’.This is the experiment conducted also in[25](but there only one dimension was used).Examples of this dataset can be seen infigure6.Correct Clusterings(out of10)Complete Linkage Euclidean34.96DTW237.6412.7338.04116.17328.85145.06565.203113.583266.753728.277Distance Time(sec)CorrectClusterings(out of45)ASL with noiseEuclidean 2.271520Figure 7.ASL data :Time required to compute the pairwise distances of the 45combinations(same for ASL and ASL withnoise)Figure 8.Noisy ASL data :The correct clusterings of the LCSS method using complete linkage.Figure 9.Performance for increasing number of Near-est Neighbors.Figure 10.The pruning power increases along with the database size.jectories.We executed a set of -Nearest Neighbor (K-NN)queries for ,,,and and we plot the fraction of the dataset that has to be examined in order to guarantee that we have found the best match for the K-NN query.Note that in this fraction we included the medoids that we check during the search since they are also part of the dataset.In figure 9we show some results for -Nearest Neigh-bor queries.We used datasets with ,and clusters.As we can see the results indicate that the algorithm has good performance even for queries with large K.We also per-formed similar experiments where we varied the number of clusters in the datasets.As the number of clusters increased the performance of the algorithm improved considerably.This behavior is expected and it is similar to the behavior of recent proposed index structures for high dimensional data [9,6,21].On the other hand if the dataset has no clusters,the performance of the algorithm degrades,since the major-ity of the trajectories have almost the same distance to the query.This behavior follows again the same pattern of high dimensional indexing methods [6,36].The last experiment evaluates the index performance,over sets of trajectories with increasing cardinality.We in-dexed from to trajectories.The pruning power of the inequality is evident in figure 10.As the size of the database increases,we can avoid examining a larger frac-tion of the database.6Related WorkThe simplest approach to define the similarity between two sequences is to map each sequence into a vector and then use a p-norm distance to define the similarity measure.The p-norm distance between two n-dimensional vectors and is defined as。

&RS\ULJKW E\ 6SDWLDO $XWRPDWLRQ /DERUDWRU\6$/$ 0HVKIUHH 0HWKRG IRU ,QFRPSUHVVLEOH )OXLG '\QDPLFV 3UREOHPV, 7VXNDQRY 9 6KDSLUR 6 =KDQJA Meshfree Method for Incompressible Fluid Dynamics ProblemsI.Tsukanov a∗,V.Shapiro a,S.Zhang ba Spatial Automation LaboratoryUniversity of Wisconsin-Madison1513University AvenueMadison,WI53706,U.S.A.b General Motors R&D CenterWarren,MI48090,U.S.A.Accepted for publication in Int.Journal forNumerical Methods in EngineeringAbstractWe show that meshfree variational methods may be utilized for solution of incompressiblefluid dynamics prob-lems using the R-function method(RFM).The proposed approach constructs an approximate solution that satisfiesall prescribed boundary conditions exactly using approximate distancefields for portions of the boundary,transfiniteinterpolation,and computations on a non-conforming spatial grid.We give detailed implementation of the methodfor two common formulations of the incompressiblefluid dynamics problem:first using scalar stream function for-mulation and then using vector formulation of the Navier-Stokes problem with artificial compressibility approach.Extensive numerical comparisons with commercial solvers and experimental data for the benchmark back-facing stepchannel problem reveal strengths and weaknesses of the proposed meshfree method.Keywords:meshfree method,distancefield,solution structure,Navier-Stokes problem,stream function,artificialcompressibility approach1Introduction1.1Towards meshfree solution of computationalfluid dynamics problemsModeling of the incompressiblefluidflow involves solution of the Navier-Stokes equations inside a geometric domain. The interaction between thefluid and the boundary of the geometric domain,in terms of the mathematical model is described by boundary conditions,formulated for viscousfluid as known velocity profile at the inlet and zero velocity at the walls.The nature of this problem makes its treatment difficult:the solution algorithm needs to incorporate two distinct types of information—(1)analytical information that describes the Navier-Stokes equations and func-tions given as boundary conditions;and(2)geometric information about boundaries where the boundary conditions are prescribed.Conventional methods of engineering analysis solve this problem by employingfirst,the spatial dis-cretization of the geometric domain(a mesh that conforms to the boundary of the geometric domain),and second,the discretization of the Navier-Stokes equations and the boundary conditions over the discretized geometry domain.The resulting approximation,therefore,unifies both functional and geometric information.Such approach,despite its wide usage,has some drawbacks.For example,it is well known that the construction of a“good”mesh is a difficult and time consuming task.In engineering practice design iterations require efficient feedback from the analysis results to the geometric model.However,employment of conforming meshes for solution of engineering problems is not quite suitable for design purposes,because the spatial grid restricts changes of the parameters of the geometric model such that it is difficult or even impossible to change the shape of the model without remeshing.∗Corresponding author.E-mail:igor@These difficulties in the conventional approaches led to the development of methods which use non-conforming1 meshes or no meshes at all.These new meshfree(sometimes they are also called meshless)methods represent a solu-tion of the problem by linear combination of basis functions which may be constructed over meshes not conforming to the shape of the geometric model[3,23,4,17,8,21,25,18,7,9].However,the employment of non-conforming spatial discretization makes the treatment of boundary conditions more difficult.Proposed remedies include the combination of Element Free Galerkin Method(EFG)[4]withfinite element shape functions near the boundary[17],the use of modified variational principle[20],window or correction functions that vanish on the boundary[9],and Lagrange multipliers.Although these techniques appear promising,they often contradict the meshfree nature of the approxi-mation near the boundary,introduce additional constraints on solutions,or lead to systems with an increased number of unknowns[13].Several promising transformation-based approaches to satisfying essential boundary conditions at desired nodal locations have been recently proposed and compared by J.-S.Chen[8].The meshing problem can be substantially simplified by employment of the Cartesian grid methods[42,1,11,5]. These methods represent the geometric model by a hierarchical set of cubical/rectangular cells that simplify computa-tion of the partial derivatives using afinite difference scheme.Instead of requiring that cells conform to the boundaries of the domain,the geometric model of the domain is approximated by quad/octtree spatial decompositions to any prescribed accuracy.This approach is accompanied by introduction of additional sources of errors and potentially exponential(in the subdivision depth)increase in computational cost.In contrast to Cartesian grid methods,immersed boundary methods[24,14,15]solve the problem using a uniform non-conforming grid of points that cover the geometric model.Influence of the boundaries and boundary conditions is accounted for by modification of the differential equation of the problem,based on special case analysis.In this paper,we describe a method that also relies on a non-conforming uniform rectangular grid,but goes sub-stantially further.All prescribed boundary conditions are satisfied exactly by transfinitely interpolating individual boundary conditions inversely proportional to the approximate distance to each boundary portion.The technique can be applied systematically to any and all boundary value problems using the theory of R-functions[30],and the result-ing interpolant can be combined with just about any standard numerical solution method.The method is demonstrated with variational methods applied to the solution of incompressiblefluid dynamics problems:first using scalar stream function formulation,and then using vector formulation of the Navier-Stokes problem with artificial compressibility approach.1.2Brief History of the MethodKantorovich showed that Dirichlet boundary conditions could be satisfied exactly using functions vanishing on the boundary of a geometric object[16].He proposed to represent a solution satisfying Dirichlet boundary conditionu|∂Ω=ϕin the following form:u=ωNi=1C iχi+ϕ,(1)whereωis a function taking on zero value on the boundary of the domain;{χi}N i=1is a system of linearly independent basis functions;{C i}N i=1is a vector of unknown coefficients andϕis a function given as a boundary condition. Different sets of the coefficients{C i}N i=1give different functions u but all of them satisfy the prescribed boundary condition.Numerical values of the unknown coefficients can be obtained via variational or projectional methods. Application of Kantorovich’s method was limited to very simple geometric domains,because at that time it was unclear how to construct functionωfor arbitrary geometric domains.Several years later,Rvachev proposed that functions taking on zero value on the boundary of a geometric domain can be constructed for virtually any geometric object using R-functions[27,28,34].Informally,R-functions serve as a construction toolkit transforming a set-theoretic description of the boundary of a geometric object into a real valued function whose zero set coincides with the boundary.Detailed discussion on R-functions and construction techniques is outside of the scope of this paper,but it can be found in numerous references[28,34,35,26,29]and will be illustrated in section2.2.Functions constructed using R-functions are differentiable everywhere except a finite number of points[28,35]and behave as distances to the boundaries near the boundary points.We will refer 1This should not be confused with the another commonly used terminology of“conforming/non-conformingfinite element”.In this paper the non-conformance of the spatial grid to the shape of the geometric domain means that the grid is extended beyond,and unconstrained by the boundary of the geometric domain.to such functions as approximate distancefields.Besides techniques based on the theory of R-functions[28],other methods may also be applied for construction of approximate distancefields.For example,the level set method [33,32]results in a distance-like functions,albeit defined at a discrete set of points and usually implicitly.In contrast, the approximate distancefields constructed using R-functions are explicitly defined at all points of the space.The successful employment of the level set method to model holes and inclusions was discussed in[38].Similar technique was used to model crack development and propagation in[37].Approximate distancefields can be used for interpolation of the functions and their derivatives prescribed on the boundary pieces of a geometric object[31].Representation of boundaries of a geometric object by approximate distancefields made possible the extension of the Kantorovich’s method into the R-function method(RFM).The RFM allows the satisfaction of many types of boundary conditions exactly by employing solution structures that incorporate boundary conditions,approximate distancefields,and basis functions with unknown coefficients[30]. RFM is essentially a meshfree method because it places no restriction on the choice of the basis functions:they can be constructed over conforming or non-conforming mesh.For example,finite elements can be used as basis functions;in this case,RFM can be viewed as an enhancedfinite element method that treats all given boundary conditions exactly. But in this paper,all computations were performed over a uniform rectangular grid of B-splines and performed within the SAGE system developed by authors[41].In[36],we showed that the method is particularly effective in dealing with moving and deforming boundaryconditions.Figure1:Parametrization of the geometry of a back facing step channel1.3Scope and outlineThis paper serves two purposes.First,the application of the RFM to Computational Fluid Dynamics(CFD)is illus-trated through two different formulations;second,the numerical properties(accuracy and convergence)of the RFM, applied to these formulations,are investigated.The paper will show that the RFM approach to CFD provides a unique andflexible method to link both the functional and geometric information in a unified manner.It will also be shown that the artificial compressibility approach gives good results for low Reynolds numberflow(Re<400).For higher Reynolds numbersflows the RFM needs to be implemented either with Finite V olume(FV)/Finite Difference(FD) numerical schemes or using the regularization approach described in[10].Throughout the paper,we assume thatflow is laminar,fluid is Newtonian,and we focus on two-dimensional problems.In this case incompressible two dimensional viscousflow is described by Navier-Stokes equations and the continuity equation[19]:u ∂u∂x+v∂u∂y−1Re∇2u=−Eu∂p∂x;u ∂v∂x+v∂v∂y−1Re∇2v=−Eu∂p∂y∂u ∂x +∂v∂y=0,(2)where variables u and v are the velocity components in the x and y coordinate directions respectively,p is thepressure variable,Re=23u max2h inletνand Eu=Pρu2maxare Reynolds and Euler numbers respectively.We explorethe accuracy of the RFM and its convergence properties,solving a standard textbook benchmark problem:an incom-pressible viscousfluidflow in a two-dimensional back-facing step channel,whose parametrization is shown in Figure 1.For this problem experimental data[2],as well as the computer simulation results given by the conventionalfluiddynamics systems,are available.For concreteness we let the geometric parameters take on the following numericalvalues:L inlet=5,L channel=12,h inlet=0.5and s=0.471.In order to simplify the comparison of the RFM modeling results with experimental data,we use the same ratio between h inlet and s as in[2].Boundary conditions areformulated as a parabolic velocity profile with u max=1.5at the inlet and zero velocity at the walls of the channel.Most meshfree methods employ some variational principle in order to solve the problem,and the RFM is no ex-ception.Since RFM treats the given boundary conditions exactly,the variational principle is applied to the differentialequation(s)of the problem only.Because viscousfluidflows do not conserve energy,we employ a least squares method.In the paper we discuss the application of the RFM to two different formulation of the incompressiblefluid dynamics problem:stream function and artificial compressibility formulations.The stream function formulation discussed in Section2substantially simplifies the initial problem reducing the system of the Navier-Stokes and continuity equations to a single equation.We use the stream function formulation as an introductory example in order to explain the concept of the RFM solution structure and the RFM solution procedure. The velocity profiles given by the RFM are in good agreement with the experimental data for Reynolds number100. For higher Reynolds numbers the RFM overestimates the position of the reattachment point.The accuracy of the modeling results can be improved by applying the RFM to the primitive variables of the Navier-Stokes equations via an artificial compressibility approach,detailed in Section3.In contrast to the stream function formulation,the artificial compressibility formulation allows modeling offluidflows in channels with arbitrary geo-metric shape including multiple connected channels.Further,it can be easily extended to model three-dimensional and turbulentflows.Since the artificial compressibility formulation leads to solution of a vector problem,application of the RFM to this formulation of incompressiblefluid dynamics problem results in vector solution structures whose construction we explain in Section3.2.Section2.5and Section3.4contain the analysis of the RFM modeling results and their comparison with exper-imental data and numerical results obtained using the commercialfluid analysis system Fluent.Distributions of the velocity components and the pressurefield given by the RFM are in good agreement with experimental data,however the employment of variational methods appears to raise some issues.These observations and possible ways to improve effectiveness of the method are discussed in Section4.2RFM with stream function formulation2.1Stream function formulation and solutionIntroduction of a stream functionψsuch that u=∂ψ∂y and v=−∂ψ∂xallows us to satisfy the continuity equation andto exclude the pressure p from the momentum equations.Substitution of the stream functionψinstead of derivatives of velocity components gives a differential equation for the stream function[19]:1 Re ∇4ψ−∂ψ∂y∂∂x∇2ψ+∂ψ∂x∂∂y∇2ψ=0.(3)Boundary conditions for the stream function at the inlet can be derived from the velocity profile at the inlet which is usually known:ψ|inlet=sV(x,y)cos(n,V)dS,(4)where V=u i+v j is the velocity vector and n is the normal vector to the inlet section.Assuming a parabolic profile for the u component of the velocity vector,as shown in Figure1,with u max=1.5and v=0at the inlet,we obtain the boundary condition for stream function as:ψ|inlet=2h2inlety3+3h inlety2;∂ψ∂n|inlet =0.(5)On the walls of the channel the stream function should satisfy the following boundary conditions:ψ|lower wall =0;ψ|upper wall =h inlet .(6)Since we are dealing with viscous flow,velocity is assumed to be zero on all walls.This condition is expressed in terms of homogeneous Neumann boundary condition for the stream function:∂ψ∂n |walls =0.(7)Boundary conditions for the stream function at the outlet can be derived from two conditions:(a)the total discharge at the outlet should be the same as at the inlet,and (b)the velocity components at the outlet should respectively be:v =0,and u should possess a parabolic profile as is shown in Figure 1.After simplification we get:ψ|outlet =m −y 33+h inlet −s 2y 2+h inlet sy +k ;m =6h inlet (h inlet +s )3;k =m 4 3s 2h inlet +s 3 ;∂ψ∂n |outlet=0.(8)The geometry of the channel,equation (3)together with boundary conditions (4–8)constitutes a complete mathe-matical formulation of the problem.The first steps of solving this problem with RFM are construction of approximate distance fields for boundary pieces and the RFM solution structure for the problem.The solution structure interpolates all given boundary conditions over the specified geometry,but also includes a set of basis functions with undetermined coefficients.The subsequent solution procedure will determine the coefficients that best approximate the govern-ing differential equation in some sense.The solution structure,as well as approximate distance fields,are usually constructed automatically without user’s intervention,but below we show all the construction details manually and explicitly.2.2Theory of R -functions and approximate distance fieldsThe theory of R -functions was originally developed in Ukraine by V .L.Rvachev and his students [28,26,29].A complete list of references through 2001can be found in [22].A brief English summary of the theory of R -functions written by Shapiro in 1988[34]is available as a technical report.An R -function is real-valued function whose sign is completely determined by the signs of its arguments.For example,the function xyz can be negative only when the number of its negative arguments is odd.A similar property is possessed by functions x +y + xy +x 2+y 2and xy +z +|z −yx |,and so on.Such functions ‘encode’Boolean logic functions and are called R -functions .Every Boolean function is a companion to infinitely many R -functions,which form a branch of the set of R -functions.For example,it is well known that min(x 1,x 2)is an R -function whose companion Boolean function is logical “and”(∧),and max(x 1,x 2)is an R -function whose companion Boolean function is logical “or”(∨).But the same branches of R -functions contain many other functions,e.g.x 1∧αx 2≡11+α x 1+x 2− x 21+x 22−2αx 1x 2 ;x 1∨αx 2≡11+α x 1+x 2+x 21+x 22−2αx 1x 2 ,(9)where α(x 1,x 2)is an arbitrary symmetric function such that −1<α(x 1,x 2)≤1.The precise value of αmay or may not matter,and often it can be set to a constant.For example,setting α=1yields the min and max respectively,but setting α=0results in much nicer functions ∨0and ∧0that are analytic everywhere except when x 1=x 2=0.Similarly,R -functionsx 1∧m αx 2≡(x 1∧αx 2)(x 21+x 22)m 2;x 1∨m αx 2≡(x 1∨αx 2)(x 21+x 22)m 2(10)(a)(b)Figure2:(a)Halfspaces that constitute a CSG representation of the channel;(b)the corresponding approximate distancefield(a)(b)Figure3:Approximate distancefields for the portions of the boundary where non-slip boundary conditions are pre-scribed:(a)for the car presented in Figure23(b);(b)for the car presented in Figure23(c)are analytic everywhere except the origin(x1=x2=0),where they are m times differentiable.Many other systems of R-functions are studied in[28].The choice of an appropriate system of R-functions is dictated by many considerations, including simplicity,continuity,differential properties,and computational convenience.Just as Boolean functions,R-functions are closed under ing R-functions,any object defined by a predicate on“primitive”geometric regions(e.g.regions defined by a system of inequalities)can now also be represented by a single inequality,or equation.The latter can be evaluated,differentiated,and possesses many other useful properties.In particular:•the functions are constructed in a‘logical’fashion and can be controlled through intuitive user-defined parame-ters;•functions can be normalized,in which case they behave as distance functions near the boundary of the object and can be differentiated everywhere[28,35];•functions can also be constructed for individual cells and cells complexes,given prescribed values for the func-tions and their gradients;•the functions can be used to define time-varying geometry and used for modeling various complex physical phenomena.Theory of R-functions provides the connection between logical and set operations on geometric primitives and analytic constructions.For every logical or set-theoretic construction,there is a corresponding approximate distance function with the above properties.Furthermore,the translation from logical and set-theoretic description is a matter of simple syntactic substitution that does not require expensive symbolic computations.For example,the geometric domain of the channel in Figure2(a)can be defined as a Boolean(Constructive Solid Geometry)combination of six primitives:Ω=(f1∪f2)∩f3∩f4∩f5∩f6,where x denotes the regularized complement of x,and individual primitives f1through f6are defined by the following inequalities:f1=y≥0;f2=x−L inlet≥0;f3=y−s≥0;f4=L inlet+L channel−x≥0;f5=h inlet−y≥0;f6=x≥0Naturally,all numeric constants can be viewed as specific values for some parameters(size,position,etc.).The constructed Boolean representation can be translated into the approximate distancefield shown in Figure2(b)using R-functions:ω=(f1∨0f2)∧0f3∧0f4∧0f5∧0f6,(11) which is also parameterized by h inlet,s,L inlet and L channel.This example clearly shows that any Boolean representation may be translated into the corresponding approximate distancefield.Similarly,boundary representation of a solid is a union of solid’s faces,each face is a subset of some surface bounded by edges,and so on.This logical description can also be directly translated into a function such that is zero for every point on the boundary and positive elsewhere. Our recent results[35,41]indicate that such functions can be constructed directly from the commercially available solid modeling representations,as well as from a variety of other geometric data structures,such as cell complexes. For example,Figures3(a)and(b)show approximate distancefields for the portions of the boundary where non-slip boundary conditions are prescribed(compare to the car shapes in Figures23(b)and(c)respectively).In the next Section we explain the usage of approximate distancefields for construction of the RFM solution structures and transfinite interpolation of the prescribed boundary conditions.2.3RFM solution structure for stream functionA solution structure is a function that satisfies exactly all prescribed boundary conditions.In general,any RFM solution structure can be represented as a sum of two functions:ψ=ψ0+ψ1(12) whereψ0satisfies homogeneous boundary conditions and contains necessary degrees of freedom in order to approxi-mate the differential equation of the problem;functionψ1interpolates the functions given as boundary conditions(5),Figure4:The RFM solution structure that satisfies boundary conditions(4-8)exactly(6)and(8).The interpolation term is constructed using the transfinite interpolation method[31]which is a generaliza-tion of the inverse distance weighting technique;it matches all specified boundary conditions and extends them inside the domain by some arbitrary but well behaved function.For the problem considered here,it takes:ψ1=ψoutletω2outlet+ψinletω2inlet+ψupper wallω2upper wall+ψlower wallω2lower wall1ω2outlet+1ω2inlet+1ω2upper wall+1ω2lower wall,(13)whereωoutlet,ωinlet,ωupper wall andωlower wall are approximate distancefields that describe outlet,inlet and walls of the channel as it is shown in Figure4.Rasing these functions to the second power assures that boundary condition∂ψ∂n|whole boundary =0is satisfied.Functionψ0in the solution structure(12)serves for approximation of the differential equation of the problem.In our case it can be represented as a product of the second power of an approximate distance to the boundary of the channelωand unknown functionΦwhose sole purpose is the approximation of the differential equation of theproblem:ψ0=ω2Φ.Sinceωtakes on zero value on boundary of the geometric domain,ψ0vanishes on the boundary together with itsfirst normal derivative.Therefore,regardless of the chosen functionΦ,functionψ0satisfies thehomogeneous boundary conditions exactly.In many practical situations functionΦcannot be determined exactly,which is why it is usually represented by a linear combination of basis functions{χi}N i=1:Φ=Ni=1C iχi.(14)The basis functions{χi}N i=1have to be smooth enough in order to approximate the differential equation of the problem. Thus,the RFM solution structure(12)may be rewritten as follows:ψ=ψoutletω2outlet+ψinletω2inlet+ψupper wallω2upper wall+ψlower wallω2lower wall1ω2outlet+1ω2inlet+1ω2upper wall+1ω2lower wall+ω2Ni=1C iχi.(15)This solution structure corresponds to the space that contains functions satisfying the prescribed boundary conditions and is sufficiently complete in the sense of being able to approximate the exact solution with an arbitrary degree of accuracy[30].Employment of the RFM solution structures to represent a solution of a physical problem offers several advantages. In particular:an RFM solution structure treats the prescribed boundary conditions exactly;an RFM solution structure contains no information about the differential equation of the problem which means that the same solution structure can be used to represent solutions of different physical problems with similar types of boundary conditions;basis functions in the solution structure can be constructed over a mesh conforming or non-conforming to a geometric model;solution structure can be easily adjusted to a new geometric model—only approximate distancefields have to be reconstructed in order to represent the boundary pieces of new geometric model;an RFM solution structure can be evaluated and differentiated at any point inside the computational domain;finally,an RFM solution structure can be integrated over the geometric model using adaptive numerical procedures[41].2.4Computation of the coefficients in the solution structureSince the RFM solution structure satisfies the given boundary conditions exactly,to solve the problem we need to find the set of the unknown coefficients{C i}N i=1in the RFM solution structure that gives the best approximation to the differential equation of the boundary value problem.Numerical values of these coefficients can be determined via variational or projectional methods.The differential equation(3)for the stream function contains non-linear terms that have to be linearized before the solution method is applied.After substitution of solution structure(12)into differentialequation (3)and application of Newton-Kantorovich linearization scheme we obtain:1Re ∇4ψn +10− ∂ψn +10∂y ∂∇2ψn 0∂x +∂ψn 0∂y ∂∇2ψn +10∂x −∂ψn +10∂x ∂∇2ψn 0∂y −∂ψn 0∂x ∂∇2ψn +10∂y−∂ψn +10∂y ∂∇2ψ1∂x −∂ψ1∂y ∂∇2ψn +10∂x +∂ψn +10∂x ∂∇2ψ1∂y +∂ψ1∂x ∂∇2ψn +10∂y=−1Re ∇4ψ1+∂ψ1∂y ∂∇2ψ1∂x −∂ψ1∂x ∂∇2ψ1∂y −∂ψn 0∂y ∂∇2ψn 0∂x +∂ψn 0∂x ∂∇2ψn 0∂y .(16)This equation is formulated for the function ψ0satisfying the homogeneous boundary conditions ψ0|∂Ω=0,∂ψ0∂n |∂Ω=0.Equation (16)is solved by an iterative algorithm,and the superscripts n +1and n in the equation denote solutions at the current and previous iterations respectively.The iterative process finishes as soon as the difference between twoconsecutive solutions becomes sufficiently small.At each iteration the least squares method is applied to equation (16)minimizing the residual of the equation:F = Ω1Re ∇4ψn +10− ∂ψn +10∂y ∂∇2ψn 0∂x +∂ψn 0∂y ∂∇2ψn +10∂x −∂ψn +10∂x ∂∇2ψn 0∂y −∂ψn 0∂x ∂∇2ψn +10∂y −∂ψn +10∂y ∂∇2ψ1∂x −∂ψ1∂y ∂∇2ψn +10∂x +∂ψn +10∂x ∂∇2ψ1∂y +∂ψ1∂x ∂∇2ψn +10∂y+1Re ∇4ψ1−∂ψ1∂y ∂∇2ψ1∂x +∂ψ1∂x ∂∇2ψ1∂y +∂ψn 0∂y ∂∇2ψn 0∂x −∂ψn 0∂x ∂∇2ψn 0∂y 2d Ω→min.(17)From the necessary condition of the existence of minimum ∂F ∂C i =0,i =1,...,N we obtain a system of linear equations AC =B whose solution gives the numerical values of the unknown coefficients in the solution structure.Elements of the matrix A and vector B are defined as follows:a ij = Ω 1Re ∇4 ω2χi − ∂∂y ω2χi ∂∇2ψn 0∂x +∂ψn 0∂y ∂∂x ∇2 ω2χi −∂∂x ω2χi ∂∇2ψn 0∂y −∂ψn 0∂x ∂∂y ∇2 ω2χi −∂∂y ω2χi ∂∇2ψ1∂x −∂ψ1∂y ∂∂x ∇2 ω2χi +∂∂x ω2χi ∂∇2ψ1∂y +∂ψ1∂x ∂∂y ∇2(ωχi ) 1Re ∇4 ω2χj − ∂∂y ω2χj ∂∇2ψn 0∂x +∂ψn 0∂y ∂∂x ∇2 ω2χj −∂∂x ω2χj ∂∇2ψn 0∂y −∂ψn 0∂x ∂∂y ∇2 ω2χj −∂∂y ω2χj ∂∇2ψ1∂x −∂ψ1∂y ∂∂x ∇2 ω2χj +∂∂x ω2χj ∂∇2ψ1∂y +∂ψ1∂x ∂∂y∇2 ω2χj d Ω;(18)b i = Ω 1Re ∇4 ω2χi − ∂∂y ω2χi ∂∇2ψn 0∂x +∂ψn 0∂y ∂∂x ∇2 ω2χi −∂∂x ω2χi ∂∇2ψn 0∂y −∂ψn 0∂x ∂∂y ∇2 ω2χi −∂∂y ω2χi ∂∇2ψ1∂x −∂ψ1∂y ∂∂x ∇2 ω2χi +∂∂x ω2χi ∂∇2ψ1∂y +∂ψ1∂x ∂∂y ∇2 ω2χi −1Re ∇4ψ1+∂ψ1∂y ∂∇2ψ1∂x −∂ψ1∂x ∂∇2ψ1∂y −∂ψn 0∂y ∂∇2ψn 0∂x +∂ψn 0∂x ∂∇2ψn 0∂yd Ω(19)Integrals (18)and (19)are computed using adaptive integration algorithm based on the Gauss-Legendre quadrature rule in conjunction with hierarchical space decomposition technique [41].。

专利名称：METHOD FOR COMPARING THE SIMILARITY BETWEEN TWO OBJECTS发明人：DEL BIANCO, Alessandro,KURZMANN,Andreas申请号：AT2008000001申请日：20080103公开号：WO08/080185P1公开日：20080710专利内容由知识产权出版社提供摘要：The method is that the infusion is stopped with temporary values of intracranial pressure while in the increasing phase, and when the intracranial pressure values are higher than those determined by the inflexion point tp, and also after the inflexion point value is exceeded, and in addition the inflexion point tp on a curve created by temporary values of the intracranial pressure in the increasing phase, is described by the formula: where: tp is the time at which the inflexion point occurs on the intracranial pressure curve in the increasing phase; whereas the intracranial pressure after the infusion has been stopped in the descending phase, is expressed by the following formula: where: ICP02 is a reference pressure of the cerebrospinal fluid in the descending phase, E2 is an intracranial elastance coefficient in the descending phase, RCSF2 is a resistance of the cerebrospinal fluid reabsorption in the descending phase, and then, from temporary values of the intracranial fluid pressure ICP(t) in the descending phase, the reference pressure of the cerebrospinal fluid ICP02 in the descending phase is determined, as is also E2/RCSF2 the quotient of the intracranial elastance coefficient E2 divided by the resistance of the cerebrospinal fluid reabsorption RCSF2 in the descending phase. The system consists ofthe intracranial pressure sensor (3) through a pressure transducer (4) that converts pressure into an analog signal, and through an analog-to- digital converter (5), is connected to a digital signal recording and processing system (6), and the infusion pump (2) is also connected to the digital signal recording and processing system (6).申请人：DEL BIANCO, Alessandro,KURZMANN, Andreas地址：AT,AT,AT国籍：AT,AT,AT代理机构：WILDHACK, Helmut更多信息请下载全文后查看。

专利名称：METROLOGY METHOD 发明人：RAKOS, Karl申请号：GB2013/052455申请日：20130919公开号：WO2014/045038A1公开日：20140327专利内容由知识产权出版社提供专利附图：摘要：An interferometric method for profiling the topography of a sample surface, said method comprising the steps of: (i)a first interferometric profiling step in which a sample surface is analysed by single-frame interferometry (SFI) at a relatively low first magnification M1 to produce a map comprising pixels with planar (X,Y)-coordinatescorresponding to the area of the sample surface, (ii)interrogating the pixel data obtained from the first profiling step by identifying pixel(s) which meet or exceed a Cut-Off Threshold, and which also meet or exceed a parameter N NAP which is the number of adjacent pixels all of which meet or exceed the Cut-Off Threshold; (iii)interrogating the pixel data obtained from the first profiling step by identifying pixel(s) for which no z-coordinate has been recorded; (iv)generating a Low Magnification Frame File (LMFF) which comprises the (X,Y) coordinates of the pixels derived from steps (ii) and (iii); (v)a second interferometric profiling step in which the sample surface is analysed at a relatively high second magnification M2, wherein M2>M1, wherein only selected regions of the sample surface are analysed at said second magnification M2, wherein said selected regions comprise the features associated with the (X,Y)-coordinates of the pixels in the Low Magnification Frame File; and further comprising a step selected from: (vi)the step of analysis of the output of the second interferometric profiling step to differentiate between an intrinsic defect and an extrinsic defect; (vii)the step of assessing whether said sample surface meets one or more quality control standard(s) and/or one or more target property or properties; and (viii)the step of assessing whether said sample surface is suitable as a surface for subsequent coating.申请人：DUPONT TEIJIN FILMS U.S. LIMITED PARTNERSHIP,RAKOS, Karl地址：3600 Discovery Drive Chester, Virginia 23836 US,2 Warren Close Darlington Durham DL1 4RA GB国籍：US,GB代理人：COCKERTON, Bruce Roger et al.更多信息请下载全文后查看。

(2+1)维extended Kadomtsev-Petviashvili方程的混合型精确解WANG Meineng【摘要】Kadomtsev-Petviashvili方程的类型非常多,描述很多数学物理现象,研究Kadomtsev-Petviashvili方程的精确解是非常有必要的.本文主要讨论(2+1)维extended Kadomtsev-Petviashvili方程.基于Hirota双线性形式和符号计算软件Mathematica,考虑指数函数,三角函数和双曲函数的混合,我们获得了(2+1)维extended Kadomt-sev-Petviashvili方程一些新的混合型精确解,并利用一些三维图形展示了这些解的物理结构和特点.【期刊名称】《南昌大学学报（理科版）》【年(卷),期】2019(043)002【总页数】5页(P123-126,150)【关键词】指数函数;extended Kadomtsev-Petviashvili方程;混合型精确解;三角函数;双曲函数【作者】WANG Meineng【作者单位】【正文语种】中文描述非线性现象的非线性偏微分方程在等离子体物理学、流体力学、孤立波理论、流体力学和湍流理论、光纤、水波、混沌理论和化学物理中有着广泛的应用[1-5]。

对这些方程的研究，可积的或不可积的，近几十年获得了蓬勃发展，旨在揭示这些方程的特点和性质。

孤子脉冲，非线性和色散效应之间的完美平衡，孤子之间的相互作用，递归算子，对称性，可积性，以及其他概念中的许多特征吸引了许多研究者的工作，许多研究非线性发展方程的有效方法不被提出来[6-15]。

Kadomtsev-Petviashvili (KP)型方程已经被许多研究者研究过，使用了很多不同的方法，如用多种方法，例如Wronskian和Grammian方法、多重EXP函数方法、简化Hirota法、TANH法和一类新的有理双曲函数方法[16]。