− z z S z z (17)
where z=(a0 ,a1)' and S is an unbiased estimate of Σ with components
2 11
(1 )
xx
S MSE x n S
= + , 22 ,
xx
S MSE
= S and 12 ( ).
xx
S MSE x
= − S
Kang and Albin (2000) proved that the T2/2 statistic in Equation (17) has the F distribution. Thus, the upper control limit of the T2 chart in Phase I is
2,( 2) ,
2 n k
UCL= F − α. For EWMA and R charts in Phase I, the control limits are modified by substituting σ with MSE . There are other methods to examine the historical data in Phase I. Stover and Brill (1998) proposed a Hotelling’s T2 approach that is similar to the T2 method of Kang and Albin (2000). The distinction between the two is the estimate of the variance-covariance matrix. Another approach is to use a univariate chart based on the first principal component corresponding to the vectors containing the estimators of the intercept and slope. But Kim et al. (2003) advised against this principal component method since it is unable to detect the shifts in the direction perpendicular to the first principal component.
3. Methodologies 3.1 B-splines
In order to extend linear profiles to any smooth functions, a smooth curve fitting technique is needed for profile smoothing from noisy data. A polynomial function is not very flexible for approximating curves with different degrees of smoothness at different locations. One way to overcome this drawback is to use locally a polynomial approximation of low degree. Another way is to allow the derivatives of the approximating function to have discontinuities at certain locations. This can be accomplished by fitting piecewise polynomials or splines. Frequently, a cubic spline, i.e., a piecewise polynomial with continuous first two derivatives, is used for such
approximation. Consider the nonparametric regression model ( ) , 1,..., ,
i i i
y =m x +ε i= n (18)
where m(x) is a regression curve and εi are i.i.d. normally distributed with zero mean and common variance . In this paper, we adopt the B-spline regression method for its popularity and simplicity. See de Boor (1978) for the definition of the B-spline basis. The points where the derivatives of the approximation function could have discontinuities are called knots. Let
2 0
σ >
1,..., b k
t t+ be the knots, where b is the number of bases and each polynomial is of order k. Then each basis of order k is a k-2 continuously differentiable function. Define Bl k, as a B-spline basis that is nonzero only on the interval ( ,t tl l k+ ),l=1, 2,..., .b A B-spline curve of order k can be constructed as
1 ,
( )i b l l k( )i i [ , ], 1,..., , P x =
∑
l= c B x ∀ ∈x u v i= nwhere ’s are the unknown B-spline coefficients and n is the number of set points in the interval [u, v]. So we can modify Equation (18) to
cl
1 , ( ) , 1,..., .
b
i l l l k i i
y =
∑
= c B x +ε i= n (19)The B-spline basis Bl k, of order k can be defined as
1 ,1
1
, , 1
1 1
1 ,
1 , ( )
0 ;
2 , ( ) ( ) ( ).
l l
l
l l
l l k
l k l k l k
l k l l k l
for t t t if k B t
for t t and t t
t t t t
if k B t B t B t
t t t t
+ +
− + +
+ − + +
⎧ ⎧ ≤ <
= = ⎨
⎪ < ≥
⎪ ⎩
⎨ − −
⎪ ≥ = +
⎪ − −
⎩ 1, −1
The following S-PLUS function “spline.des” generates a matrix of B-spline bases evaluated at points provided by users.
spline.des(knots, x, ord, derives),
where “knots” is the vector of knots for the spline, “x” represents x-coordinates at which to evaluate the spline basis functions, “ord” is the order of the spline, and
“derives” is the order of the derivative to evaluate at each of the points. If this vector is given, it must have the same length as x. The default is a vector of zeros of the same length as x. We show an example of B-spline bases in Figure 2.
After constructing B-spline bases, we use these normalized B-spline functions
,
Bl k as the regression basis function. For the given knots, the spline regression method finds the best spline approximation via the following least squares regression:
2 ,
1 1
{ (
n b
i l l k i
i l
min y c B x
= =
∑
−∑
c )} , (20)
where c=( ,...,c1 cb) '. So the least-squares estimator of c is ˆ=( -1 '
c B'B ) B y
wherey=(y1,...,yn)' , cˆ=(c1,...,cb) ', and B is the design matrix with the ( i , l )th element Bl k, ( ),xi l=1,...,b , i=1,...,n . Then c has a multivariate normal distribution with the mean vector c and the variance-covariance matrix
ˆ
2 1
( )
σ −
=
Σ B'B .
In phase II, it is assumed that the reference profile is known. Denote it by . We now establish a B-spline representation of the reference profile. First obtain n pairs of data,
( ) f x
{( , ( )),x f xi i i=1,..., }n . Let be the least squares solution of (20) with replaced by
c
yi f x . We then treat this vector c as the true in-control value. For the ( )i jth sample profile, compute the vector of sample estimators cj =(c1j,...,cb j) ', where
1j,..., b j
c c are the estimated sample B-spline coefficients. Then the T2 statistic of the jth sample is given by . But the upper control limit in phase II needs to be modified as UC
2 1
( j ) ' ( j
Tj = c −c Σ c− −c) L=χb2,α
1
. And in Phase I, we also use historical data containing k sample profiles to construct the control chart. Then the reference curve can be estimated by y B cˆ = with c=( ,...,c cb) ' and σˆ2 =MSE where
1
1 1
1 ,..., ,
k k
j b j
j j
b
c c
c c
k k
= =
=
∑
=∑
and 1
k j j
MSE
MSE k
=
∑
=(21)
with MSEj = yj−B cˆj 2/(n b− ) and ‧ denotes the Euclidean norm. The T2 statistic for the jth sample profile is then modified by
2 1
0 ( ) ' (
1 j j
j
T k k
= − − − ) ,
− c c S c c (22) where is an unbiased estimate of . Thus, the upper control limit of the T
ˆ (2 )
σ −
=
S B'B 1 Σ
2 chart in Phase I is UCL=bFb n b k,( − ) ,α.
The choice of the number of bases, as in the role of the smoothing parameter in any nonparametric regression methods, is an important issue. The boundary effect is another issue. We will address these two issues in Section 4.
3.2 Residual EWMA and R Charts
The second approach we propose is to use the EWMA and the R chart to monitor the residual average and the range, respectively, instead of the combined EWMA/R chart proposed by Kang and Albin (2000). Kang and Albin (2000) found that the value L = 3.1151 yields an in-control ARL of approximately 802 for the EWMA chart and an in-control ARL of approximately 261 for the R-chart. Although ARL0 of their combined procedure is close to 200, the detecting power of the EWMA chart is unequal to that of the R chart. In order to increase the detecting powers of both the EWMA and the R charts, we decide to monitor the process by using the EWMA chart to detect mean shifts and the R chart to detect standard deviation shifts of the residuals.
The regression residual vector for the jth sample is
ˆ { '} .
j = j− j = − -1
e y y I B(B'B) B yj
The details of the EWMA and the R charts are the same as that described in Section 2.
3.3 Metrics
The third approach is to use some metrics to monitor the process. Gardner et al.
(1997) presented a new methodology based on some “metrics” for the equipment fault detection. The equipment fault detection was not only in that it incorporates the use of integrated spatial information in a virtual wafer surface that was fitted by thin-plate splines, but also in that it can be used to detect and classify equipment faults at the same time. Their main focus was to use the differences between the observed and the expected virtual wafer surfaces to construct metrics which can be used to detect and diagnose various types of equipment faults. Thus, the following general and specific metrics were designed:
1 ( )2
M =
∫
R g T dR− (Squared metric) ;2 R
M =
∫
g T dR− (Absolute-value metric) ;2
3 ( ) 1.5
0
o
M R g T dR if g T A otherwise
= − − >
=
∫
(Spec-limit metric);4 ( )2 ( )
( )
R
R
M g T dR if g T
g T dR if g T
= − − >
= − − − ≤
∫
∫
0 0
(Square-above-absolute-value-below metric);
5 ( )2 ( )
( )
R
R
M g T dR if g T
g T dR if g T 0 0
= − −
= − − − >
∫
∫
≤
(Square-below-absolute-value-above metric);
where g is a newly fitted thin-plate spline surface, T is the target surface, and R denotes the wafer surface region. Metrics 1-3 are general metrics used to detect the presence of an equipment fault, and Metrics 4-5 are specific metrics used to detect specific fault patterns. Gardner et al. (1997) proposed an alternate Bayesian (simulation) approach that can be taken to determine the null distribution of the
metrics. According to a procedure given in Green and Silverman (1994), assuming a Gaussian prior distribution, the posterior distribution of the thin-plate spline surface g has the following multivariate normal distribution:
2 ˆ
~ MVN [ ,ˆ ˆσ ( )]λ
g g A ,
where gˆ is the vector of fitted values, σ is calculated as the residual sums of ˆ2 squares about the fitted curve divided by an effective degrees of freedom, and A( )λˆ is the projection matrix which maps the vector of observed values to their predicted values. Simulate M independent sets of n observations from the posterior multivariate normal distribution described above. A spline surface is fitted to each set of observations, and the metrics are calculated. As a result, M independent values are obtained from the null distribution of each metric. The (100×α)th percentile of the M simulated metric values is the critical value for determining if the new curve is far from the target curve. For each newly observed curve, metrics are calculated and compared to their corresponding critical values in order to determine whether or not a fault is detected. In the experimental data, they used the average predicted surface of two wafers as the target surface. 5,000 observations were simulated using parametric bootstrapping from MVN[ ,gˆ ˆσ2A( )]λˆ where is the average of the predicted values using the thin-plate spline fitting for these two wafers, and
gˆ
2 ˆ
ˆ ( )
σ A λ is the average of the covariance matrices from the thin-plate spline fitting for these two wafers. Then these simulated observations were used to obtain 5,000 metric values and the corresponding empirical α=0.01 “critical value” for each metric.
With the same idea, we select two different and popular metrics to monitor the process in all types of shifts. They are Squared metric and Absolute-value metric. Two metrics are defined in the following:
2 1
1 { ( ) ( )}
n
i i
i
M g x T x
=
=
∑
− (Squared metric);1
2 | ( ) ( ) |
n i i
M g x T xi
=
=
∑
− (Absolute-value metric);where g is a newly fitted B-spline curve and T is the reference profile. To construct the control charts, we need to find the critical value of the null distribution of each metric.
Since the distributions of these metrics are difficult to obtain, we use historical data (k sample profiles) to compute the reference profile in Phase I. Fit a B-spline to each of the k sample profiles. Denote the average of k regression estimators for the lth B-spline coefficient by cˆl . So the estimated reference profile , and
1
ˆ ( )
b l l k l
c B x
∑
=ˆ2 MSE
σ = are available by Equation (21). Simulate M sets of n observations from the following model:
1ˆ , ( ) , 1,..., ,
b
i l l l k i i
y =
∑
= c B x +ε i= n )(23) where
. . .
ˆ2
~ (0 ,
i i d
i N
ε σ . For each metric, compute the metric value for each of the simulated profiles. Let the (100×α)th percentile of these M values of the metric be the critical value. The process is claimed out-of-control when the metric of the newly observed profile is greater than the critical value.
3.4 EWMSD Chart
The last approach is using exponentially weighted moving standard deviation (EWMSD) to monitor the process. It is very sensitive in detecting shifts in process variation, particularly when shifts are relatively small. Define the sample standard deviation of the residuals for the jth sample as
ˆ 2
, 1, 2,....
j j
sj j
n b
= − =
− y y
The EWMSD statistic is given by
(1 ) 1
j j
v =θs + −θ vj− ,
where θ ( 0< ≤ ) is a smoothing constant and θ 1 v0 =σ0. In Phase II, assume the in-control value of σ is known. The control limits for the EWMSD chart are 2
(2 )
LCL L
n σ σ θ
= − θ
− and
(2 )
UCL L
n σ σ θ
= + θ
− ,
where L is also a constant chosen to give a specified in-control ARL. When the EWMSD statistic is beyond the control limits, the process is claimed out-of-control.
In Phase I, the control limits are modified by substituting σ with MSE .