Existence Theorem on Quasiconformal Mappings

Seddik Gmira Quasiconformal mappings are, nowadays, recognized as a useful, important, and fundamental tool, applied not only in the theory of Teichmüller spaces, but also in various …elds of complex analysis of one variable such as the theories of Riemann surfaces, of Kleinian groups, of univalent functions. In this paper we prove the existence theorem of the solution of the Beltrami di¤erential equation, and we give a fundamental variational formula for quasiconformal mappings, due to L.Ahlfors and L.Bers.

Quasiconformal Mapping

We consider an orientation-preserving homeomorphism f , which is at least partially di¤erentiable almost every where on a domain D in C, satisfying the Beltrami equation

f z = f z
As a natural generalization of the notion of conformal mappings we consider the following

Analytic De…nition

De…nition 1 Let f be an orientation-preserving homeomorphism of a domain D into C. f is quasiconformal (qc) on D if 1. f is absolutly continuous on lines (ACL) 2. There exists a constant k, 0 k < 1 such that jf z j k jf z j almost every where on D:

Setting K = (1 + k) = (1 k),
we say that f is K-quasiconformal mapping on D. We call the in…mum of K > 1 such that f is qc, the maximal dilatation of f , and denote it by K f :

Example 1 An a¢ ne mapping f (z) = az+bz+c; (a; b; c 2 C; jbj < jaj) is qc: k = jbj = jaj and hence K f = (jaj + jbj) =(jaj jbj)

Example 2 Set f (z) = z 1 jzj 2 , z 2 . This f is an orientationpreserving of the unit disk but not quasiconformal.

In general, existence of the partial derivatives f z and f z is not enough to guarantee good properties of f applicable to further investigations. However in the case of homeomorphisms Gehring and Lheto obtained the following remrkable result.

Proposition 1 If a homeomorphism f of a domain D onto C has the partial derevatives f x and f y almost.every where on D, then f is totally di¤erentiable almost.every where on D.

Proposition 2 Let f be a quasiconformal mapping of a domain D. Then the partial derevatives f z and f z are totally square integrable on D. Proof. Let A (E) be the area of f (E) of a measurable Borel subset E of D, and J f (z) the density of the set function A(with respect to the Lebesgue measure dxdy). The Lebesgue theorem implies

Z E J f (z) dxdy A (E)
Then, f is totally di¤erentiable at almost every z 2 D and at each such a point z we have

J f (z) = jf z j 2 jf z j 2 : Since f is quasiconformal, then jf z j 2 jf z j 2 1 1 k 2 J f (z) a.e. on D
Proposition 3 For every quasiconformal mapping f of a domain D, the partial derivatives f z and f z are coincident with those in the sens of distribution. Namely for every element ' 2 C 1 0 (D), the set of all smooth functions on D with compact supports it follows that

Z Z D f z :'dxdy = Z Z D f:' z dxdy Z Z D f z :'dxdy = Z Z D f:' z dxdy
Lemma 1 For a given p > 1, let f be a continuous function on a domain D whose distributional derivatives f z and f z are locally L p on D. Then for every compact subset F of D, there is a sequence

ff n g 1 1 in C 1 0 (D) such that f n converges to f uniformly on F , with lim n !1 Z Z F j(f n ) z f z j p ddxdy = 0 lim n !1 Z Z F j(f n ) z f z j p ddxdy = 0 2 Proof. Fix 2 C 1 0 (D) with = 1 in some neighbourhood of F .
Then f has a compact support. Further ( f ) z and ( f ) z exist and belong to L p (D). Set

' (z) = ( C: exp 1 1 jzj 2 ; z 2 0; z 2 C
Here, is the unit disk. Choose a constant C so that Z Z ' (z) dxdy = 1

Further, for every n > 0, set ' n (z) = n 2 ' (nz) ; z 2 C and

f n (w) = ' n ( f ) (w) = Z Z C ' n (w z) ( f ) (z) dxdy; w 2 C
Then for every n, we have

f z = ' n ( f ) z and f z = ' n ( f ) z . Morever, f n 2 C 1 0 (D)
for every su¢ ciently large n, and f n converges to f uniformly on F as n ! 1. Since ( f ) z = f z and ( f ) z = f z on F and since we can show equalities:

lim n !1 Z Z D j(f n ) z ( f ) z j p dxdy = 0; lim n !1 Z Z D j(f n ) z ( f ) z j p dxdy = 0
Lemma 2 (Weyl) Let f be a continous function on D whose distributional derivative f z is locally integrable on D. If f z = 0 in the sens of distributions on D, then f is holomorphic on D. Proof. For an arbitrarily relatively subdomain D 1 of D, we construct an L 1 -smoothing sequence ff n g 1 1 for f with respect to D 1 as in the proof of Lemma 1, we see that ( f ) z = 0 in some neighbourhood of D 1 , we see also that (f n ) z = 0 on D 1 for every su¢ ciently large n. Since f n converges to f uniformly on D 1 as n ! 1, f is holomorphic on this arbitrary D 1 .

Geometric De…nition

A quadrelateral (Q; q 1 ; q 2 ; q 3 ; q 4 ) is a pair of a Jordan closed domain Q and four points q 1 ; q 2 ; q 3 ; q 4 2 @Q, which are naturally distinct and located in this order with respect to the positive orientation of the boundary @Q.

Proposition 4 For every quadrelateral (Q; q 1 ; q 2 ; q 3 ; q 4 ), there is a homeomorphism h of Q onto some rectangle R = [0; a] [0; b] (a; b > 0) which is conformal in the interior IntQ and satis…es h (q 1 ) = 0; h (q 2 ) = a h (q 3 ) = a + ib; h (q 4 ) = ib Morever a=b is independant of h. The value a=b is called the module of the quadrelateral (Q; q 1 ; q 2 ; q 3 ; q 4 ) and denotet by M (Q). Proof. The Riemann mapping theorem implies the existence of a conformal mapping h 1 : IntQ ! H (upper half-plane). By Carathéodory's Theorem h 1 can be extended to a homeomorphism of Q onto H [ R. Using a suitable Möbius transformation, we may assume that

h 1 (q 1 ) = 1; h 1 (q 2 ) = 1 h 1 (q 3 ) = h 1 (q 4 ) > 1
Set k = 1=h 1 (q 3 ), and

h 2 (w) = Z w 0 dz p (1 z 2 ) (1 k 2 z 2 ) ; z 2 H
Then, h 2 is a conformal mapping of H onto the interior of some rectangle

[ K; K] [0; K 0 ] (K; K 0 > 0). Hence we see that h (z) = h 2 h 1 (z) + K; z 2 Q is a desired mapping. Now, let e h : Q ! e R = [0; e a]

Since

e h (q 1 ) = h (q 1 ) = 0; h (q 1 ) > 0 and e h (q 2 ) > 0 we conclude that c > 0 and d = 0. Hearafter, for every quadrelateral (Q; q 1 ; q 2 ; q 3 ; q 4 ) and a homeomorphism f of Q onto C, we consider f (Q) as a quadrelateral with vertices f (q 1 ) ; f (q 2 ) ; f (q 3 ) ; f (q 4 ). Lemma 3 Every K-qc mapping f of a domain D satis…es

1 K M (Q) M (f (Q)) KM (Q) Proof. Fix mappings h : Q ! R = [0; a] [0; b] and e h : f (Q) ! e R = [0; e a]
Z a 0 @F @x (x + iy) dx Z a 0 (jF z j + jF z j) dx Since R R R J F dxdy A e R = e
(e ab) 2 Z Z R (jF z j + jF z j) dxdy 2 Z Z R 0 jF z j + jF z j jF z j jF z j dxdy Z Z R 0 jF z j 2 jF z j 2 dxdy Z Z R 0 Kdxdy Z Z R 0 J F dxdy K (ab) e a e b where R 0 = fw 2 R : F z (w) 6 = 0g. Thus we have M (f (Q)) KM (Q).
Next replacing F by i e h f (ih) 1 (or considering (Q; q 2 ; q 3 ; q 4 ; q 1 )), the same argument gives

1 M (f (Q)) K M (Q) De…nition 2 A homeomorphism f of a domain D into C which preserves the orientation is quasiconformal on D, if there is a constant K 1 such that, M (f (Q)) KM (Q) holds for every quadrelateral Q in D Theorem 1 Let f be a K f -quasiconformal mapping of D onto e D, then 1. The inverse f 1 is also K f -quasiconformal 2.
K-quasiconformally is conformally invariant: Namely conformal mappings h and e h of domains D and e D respectively, the composed mapping e h f h is also K f -quasiconformal 3. For every K g -qc mapping g of f (D), the composed mapping g f is

K f K g -quasiconformal.
Proof. 1. Lemma 3 gives

1 K f M f 1 (Q) M (Q)) K f M f 1 (Q ) 2.
By the conformal invariance of the module 3. Clear from the de…nition Lemma 4 Every qc-mapping

f of a domain D satis…es Z Z E jf z j 2 jf z j 2 dxdy = A (E)
for every subset E of D, where

A (E) = R R E dxdy: Proof. Let E be a rectangle contained in D. If f is absolutly continuous on the boundary @E, then in view of Proposition 8, we …nd L 2 -smoothing sequence ff n g 1 1 for f with respect to E (Lemma 1). Set f n = u n + iv n . Green's lemma implies Z Z E n (u m ) x (v n ) y (u m ) y (v n ) x o dxdy = Z @E u m dv n Let m ! 1 and n ! 1, we obtain Z Z E n (u) x (v) y (u) y (v) x o dxdy = Z @E udv
Here we write f = u + iv. The right hand side is interpred as the line integral of udv along the Jordan curve @f (E). By assumption, @f (E) is recti…able. Hence we can show that

Z @f (E) udv = Z Z f (E) dudv = A (E)
Since f is ACL, every rectangle contained in D can be approximated by such rectangles.

Proposition 5 If f is quasiconformal on a domain D, then f z 6 = 0 almost.every where on D. Proof. The set E = fz 2 D : f z = 0g is measurable and f z = 0 a.e. on E. Hence.f (E) has area zero. Since f 1 is also quasiconformal, then E = f 1 (f (E)) has area zero Next, for every quasiconformal mapping f of a domain D, we can consider

f = f z f z
almost every where on D. This f is a bounded measurable function and satis…es ess:

sup z22D f (z) K f 1 K f + 1 < 1
We call f the complex dilatation of f on D. Proposition 6 For every quasiconformal mappings f and g of a domain D,

g f 1 f = f z f z g f 1 f g almost every where on D Proof. g f 1 is quasiconformal on f (D) by Theorem 1. g f 1 is di¤er- entiable on f (D) by Proposition 1. Hence g f 1 is di¤erentiable on f (D)
except for a subset E of measure zero. f 1 (E) is also of measure zero by lemma 4. Hence Theorem 1 implies that, f and g f are di¤erentiable at z and f (z), respectively a.e on D. At such a point z, the chain rule is valid. Writing w = f (z), we get

g z = g f 1 w f:f z + g f 1 w f:f z g z = g f 1 w f:f z + g f 1 w f:f z
By similar argument as before, we can show that f z 6 = 0, g 6 = 0, and g f 1 w f 6 = 0

for almost every z 2 D.

Existence Theorem on Quasiconformal Mappings

We have seen that a quasiconformal mapping f of a domain D induces a bounded measurable function f on D, which satis…es

ess: sup z22D f (z) < 1
Next, we shall prove the converse. Namely, for every measurable function with ess: sup z22D j (z)j < 1, we construct a quasiconformal mapping, whose complex dilatation is equal to . We consider the complex Banach space L 1 (D) of bounded measurable functions on a domain D, with the norm

k k 1 = ess: sup z22D j (z)j Consider the set B (D) 1 = f 2 L 1 (D) ; k k 1 < 1g of Beltrami coe¢ - cients on D.
Proposition 7 Let 2 B (D) 1 . If there exists a quasiconformal mapping f with the complex dilatation f = , then for every conformal mapping h of f (D), the mapping h f has the same complex dilatation . Conversely, for every quasiconformal mapping g with g = , the map g f 1 is a conformal mapping of f (D). Proof. As before, we have h f = f = , and g f 1 = 0 almost.every where on f (D). It follows that g f 1 is 1-quasiconformal and hence is conformal.

To solve the Beltrami di¤erential equation f z = f z ; 2 B (C) 1 , consider …rst, solving the @-problem. If we get a suitable representation f = G (f z ), then we have a relation

f z = G (f z ) z = G ( f z ) z Rewriting this relation in the form f z = F ( ), we obtain a solution f = G ( F ( ))
of the Beltrami equation. On the other hand to reconstruct f from f z we use the classical Cauchy transfotmation. For every p with 1 P < 1, we consider the complex Banach space L P (C) of all measurable functions f on C such that Z Z C jf j p dxdy < 1

The following Pompeiu's formula is essential to solve the @-problem. Proposition 8 Fix p, 2 < p < 1, and let f be a continuous function on C such that, f z and f z 2 L P (C). Then f satis…es

f ( ) = 1 2 i Z @D f (z) dz z 1 Z Z D f z (z) z dxdy, 2 D for every open disk D in C. Proof. Let q < 2 such that, 1 p + 1 q = 1. Since j1= (z &)j q 2 L 1 (D)
, the assertion follows by Hölder's inequality. Take an L p smoothing sequence ff n g 1 1 for f with respect to D and …x a point & 2 D . For every n , Green's formula gives

f n ( ) = 1 2 i Z @D f n (z) dz z 1 2 i Z Z D (f n ) z (z) z dz ^dz Since f n ! f uniformly on D, and (f n ) z ! f z in L p (D) respectively.
Next, we de…ne a linear operator P on L P (C) as follows

P h ( ) = 1 Z Z C h (z) 1 z 1 z dxdy
Lemma 5 For every p with 2 < p < 1 and for every h 2 L p (C), P h is a uniformly Hölder continuous function on C with exponent (1 2p) and satis…es P h (0) = 0. Morever (P h) & = h on C in the sens of distribution. Proof. Hölder equality implies

jP h (&)j 1 khk p (z &) q < 1
Further, if & 6 = 0, by changing the variable, we have

Z Z C 1 z (z &) q dxdy = j&j 2 2q Z Z C 1 z (z 1) q dxdy
Hence there is a costant K p depending only on p such that jP h (&)j K p khk p j&j 1 2=p , & 2 Cr f0g Since P h (0) = 0 by de…nition, this inequality is valid even when & = 0: Now set h 1 (z) = h (z + & 1 ). Then we have

P h 1 (& 2 & 1 ) = P h (& 2 ) P h (& 1 )
We conclude that,

jP h (& 2 ) P h (& 1 )j K p khk p j& 2 & 1 j 1 2p
Thus P h is a uniformly Hölder continuous with exponent 1 2p: For the second assertion we take a sequence fh n g 1

1 in C 1 0 (C) such that kh h n k p ! 0 as n ! 1. Then for every h n (P h n ) & (&) = 1 @ @& Z Z C h n (z + &) z dxdy = 1 Z Z C (h n ) z (z) z & dxdy
Hence Green's formula implies

(P h n ) & (&) = lim " !0 Z fjz &j="g h n (z) z & dz = h n (&)
In particular, for every ' 2 C 1 0 (C), we get

Z Z C h n 'dxdy = Z Z C P h n ' z dxdy, ' 2 C 1 0 (C)
Since kh h n k p ! 0 as n ! 1, P h n converges to P h locally on (uniformaly on any compact subset) C by jP h (&)j K p khk p j&j 1 2=p . Hence, when we let n ! 1 the above equality gives

Z Z C h:'dxdy = Z Z C P h:' z dxdy, ' 2 C 1 0 (C)
Next, we need a suitable integral representation for (P h) z . For this purpose, let h 2 C 1 0 (C). And Green's formula gives

(P h) & (&) = 1 2 i Z Z C h z (z) z & dz ^dz = lim " !0 ( 1 2 i Z fjz &j="g h (z) z & dz + 1 2 i Z Z fjz &j>"g h (z) (z &) 2 dz ^dz )
Since the …rst term in the right side converges to zero as " ! 0, the second term is essential. Let T be the linear operator de…ned by

T h (&) = lim " !0 ( 1 Z Z fjz &j>"g h (z) (z &) 2 dxdy ) , h 2 C 1 0 (C) Lemma 6 Every h 2 C 1 0 (C) satis…es (P h) z = T h, on C and kT hk 2 = khk 2
Proof. We have already seen (P h) z = T h on C for every h 2 C 1 0 (C) and that

kT hk 2 2 = 1 2i Z Z C (P h) z P h z dz ^dz = 1 2i Z Z C (P h) P h zz dz ^dz = 1 2i Z Z C (P h) h z dz ^dz = 1 2i Z Z C h (P h) z dz ^dz = khk 2 2
We see that, the operator T is extended to a bounded linear operator on L 2 (C) into itself with norm 1. Since P is an operator on L p (C) with p > 2,

T is also an operator on L 2 (C). Then, we see by the following classical Calderon-Zygmund's theorem that T gives a bounded linear operator on

L P (C) (p > 2) into itself.
Theorem 2 (Calderon-Zygmund) For every p with 2 p < 1

C p = sup h2C 1 0 (C), khk p =1
kT hk p is …nite. Hence the operator T is extended to a bounded operator of L p (C). Morever, C p is continuous with respect to p. In particular C p satis…es

lim p !2 C p = 1
In fact "Calderon-Zygmund's Theorem" gives the following Proposition 9 For every number p > 2 and every h 2 L p (C)

(P h) z = T h on C in the sens of distribution. Proof. Let fh n g 1 1 be a sequence in C 1 0 (C) approximating h in L P (C). Then Z Z C T h n :' dxdy = Z Z C P h n :' z dxdy, ' 2 C 1 0 (C)
Here P h n ! P h locally uniformally on C and T h n ! T h in L p (C) respectively, as n ! 1. Hence we have the assertion.

Existence of the Normal Solutions

Theorem 3 Fix k such that 0 k < 1 arbitrarily and take p > 2 with kC p < 1. Then for every 2 B (C) 1 with a compact support; k k 1 k, there exists a unique continuous function f such that

f (0) = 0; f z 1 2 L p (C)
and satisfying f z = f z on C in the sens of distribution. This f is called the normal solution of the Beltrami equation for :

Proof. First, we drive a condition which the partial derivative f z of the normal solution should satisfy. Since f z = f z has a compact support, and since f z 1 2 L p (C), f z 2 L p (C) also. Thus we can consider P (f z ).

Set

F (z) = f (z) P (f z ) (z) , z 2 C
Then F (z) is continuous and F (0) = 0. Moreover F z = 0 in the sens of distribution. Hence F (z) is holomorphic on C by Weyl's lemma. On the other hand, since f z 1 and T (f z ) (z) belong to L p (C), so does F 0 1. Thus we conclude that

F 0 (z) = 1, i.e., F (z) = z + a, a 2 C. Since f (0) = 0, then we have f (z) = P (f z ) (z) + z, z 2 C
Then we obtain

f z = T ( f z ) + 1
Suppose that, there is another normal solution g. Then

g z = T ( g z ) + 1
By Calderon-Zygmund's Theorem, we obtain

kf z g z k p = kT ( f z ) T ( g z )k p kC p kf z g z k p
Since kC p < 1 by assumption, we get f z = g z a.e on C. Then the Beltrami equation also gives f z = g z a.e on C. Hence again by Weyl's lemma f g and f g are holomorphic on C, which in turn implies that f g should be constant. Since f (0) = g (0) = 0 we conclude that f = g.

Finally, the existence of the normal solution follows from f z = T ( f z )+1. In fact repeat substituting the whole right hand side for f z on the tight hand side. Then we …nd

f z 1 = T ( f z ) = T ( (1 + T ( f z ))) = T + T ( T ( f z )) = T + T ( T ) + T ( T ( T )) + :::::: This series converges in L p (C), since the operator norm of h ! T ( h) 2 L p (C) is not greater than kC p < 1. Set h = T + T ( T ) + T ( T ( T )) + :::::; then h 2 L p (C)
We shall show that, the desired solution for the Beltrami equation is

f (z) = P ( (h + 1)) (z) + z
In fact (h + 1) 2 L p (C), for has compact support. Hence Lemma 5 implies that f is continuous, f (0) = 0 and f z = (h + 1). Morever, by Proposition 9 we have:

f z = T ( (h + 1)) + 1 = h + 1
Hence, f satis…es the Beltrami equation f z = f z , and f z 1 2 L p (C) :

Basic Properties of Normal Solutions

Corollary 1 Under the conditions of Theorem 3, the following inequalities hold:

kf z k p 1 1 kC p k k p , and jf (& 1 ) f (& 2 )j K p 1 kC p k k p j& 1 & 2 j 1 2p + j& 1 & 2 j for every & 1 ; & 2 2 C.
Proof. Since h = T ( h) + T , then we have:

khk p kC p khk p + k k P
The …rst inequakity holds for f z = (h + 1). For the second we have f (z) = P (f z ) (z) + 1 as before,

jf (& 1 ) f (& 2 )j jP (f z ) (& 1 ) P (f z ) (& 2 )j + j& 1 & 2 j
From this, the second inequality holds easly. Furthermore, the normal solutions depend on the Beltrami coe¢ cients as follows.

Corollary 2 For 0 k < 1 and p > 2. Let f n g 1 1 be a sequence in B (C) 1 with the following conditions:

1. k n k 1 < k for every n 2. every n has a support contained in fz 2 C : jzj < M g with a suitable constant independent of n 3. n converges to some 2 B (C) 1 a.e. on C as n ! 1

Let f n be the normal solution for , and f be the normal solution for . Then f n ! f uniformaly on C as n ! 1, and

lim n !1 k(f n ) z f z k p = 0: Proof. First, since f z = T ( f z ) + 1, then we have k(f n ) z f z k p kT ( n (f n ) z f z )k p + kT ( n ) f z k p kC p k(f n ) z f z k p + C P k( n ) f z k p Hence we have k(f n ) z f z k p C p 1 kC p k( n ) f z k p
Since the support of all n are uniformaly bounded, and since n converges to a.e.on C as n ! 1, then we have

lim n !1 k(f n ) z f z k p = 0 Next, we get jf (&) f n (&)j = jP (f z (f n ) z ) (&)j K p n k( n ) f z k p + k k(f n ) z f z k p o j&j 1 2 p for every & 2 C. Thus f n ! f uniformaly on C as n ! 1. Since f n ! f is
holomorphic in a …xed neighbourhood of 1 for every n, we conclude that f n converges to f uniformaly on C.

Existence Theorem

In fact the existence of a quasiconformal mapping is also valid for a general complex dilatation 2 B (C) 1

Theorem 4 For every Beltrami coe¢ cient 2 B (C) 1 , there exists a homeomorphim f of C onto itself which is a quasiconformal mapping of C with complex dilatation . Morever f is uniquely determined by the following normalization conditions f (0) = 0, f (1) = 1 and f (1) = 1

We call this f the canonical -quasiconformal mapping of C.and denote it by f :

Proof. The uniqueness of f is given by Proposition 7 and the normalization conditions. For the existence, we suppose that has compact support. Let F be the normal solution for . Then Theorem 3 implies that F (z) =F (1) is the desired one. Now suppose that = 0 almost every where in some neighbourhood of the origin. Pulling back by the Möbus transformation

(z) = 1 z If we set e (z) = 1 z z 2 z 2 ; z 2 C
then e 2 B (C) 1 , and has a compact support. Hence as before, there exists the canonical e -quasiconformal mapping f e of C. At every such point 1 z , the quasiconformal mapping

f (z) = 1 f e (1=z)
is also totally di¤erentiable, so using the chain rule we have

f (z) = z 2 z 2 e 1 z = (z) , a.e. on C
Clearly, f sats…es the normalization conditions. Hence f is the desired function. Finally, suppose that is a general Beltrami coe¢ cient. Now, suppose that is a general Beltrami coe¢ cient. We set

1 (z) = (z) ; z 2 C 0; z 2
where is the unit disk. Then f 1 exists as before. Finally Set

2 = 1 1 1 (f 1 ) z (f 1 ) z ! (f 1 ) 1
Then, f 2 exists, because 2 has a compact support. Mreover g = f 2 f 1 is quasiconformal and we can see that g = a.e. Clearly g is the desired function.

Hereafter, we state several applications of the existenc theorem Proposition 10 Every quasiconformal mapping of the disk onto a Jordan domain D is extended to a homeomorphism on onto D Proof. Fix such a quasiconformal mapping f : ! D, and set = f . By setting = 0 on C , we can consider 2 B (C) 1 . Hence by Theorem 3, there exists the canonical -quasiconformal f of C. Set g = f f 1 . Then g is a 1-quasiconformal mapping of D. Hence g is a conformal mapping of D. Since f ( ) is a Jordan domain, Carathéodory's theorem gives the extension g to a homeomorphism of D onto f ( ). Since f = g 1 f , we obtain the assertion.

Proposition 11 There exists no quasiconformal mapping of the disk onto C Proof. Suppose that f : ! C is a quasiconformal conformal mapping Then f 1 is also quasiconformal. Set = f 1 ; then there exists the canonical -quasiconformal mapping of . On the other hand, since g 1 (C) = , Liouville's Theorem implies that g 1 should be a constant.

Proposition 12 Let be an element of B (H) 1 . Then there exists a quasiconformal mapping w : H ! H with complex dilatation : Moreover, such a mapping w(which can be extended to a homeomorphism of H = H [ R onto itself)is uniquely determined by the following normalization conditions: w (0) = 0, w (1) = 1, and w (1) = 1

We call this unique w, the canonical -quasiconformal mapping of H, and denote it by w . Proof. The uniqueness is given by the normalization conditions as before. For the existence, set In particular we see that f e R = R. Since f e preserves orientation, f e (H) = H. Hence the restriction of f e onto H is the desiret qc mapping.

Dependence on Beltrami coe¢ cients

Some of the most important useful facts on the canonical quasiconformal mapping of C concern dependence of the canonical quasiconformal mapping on the Beltrami coe¢ cient.

Theorem 5 Let f (t)g be a family of Beltrami coe¢ cients depending, on t 2 R or C. Suppose that k (t)k 1 ! 0 as t ! 0, and that (t) is di¤erentiable at t = 0 : (t) (z) = t (z) + t (t) (z) ; z 2 C with suitable 2 L 1 (C) and (t) 2 L 1 (C) such that k (t)k 1 ! 0 as t ! 0. Then Proof. Set f t = f (t) (f ) 1 . Then the complex dilatation (t) of f t is given by:

f [ ] (&) = lim
(t) = ft = (t) 1 : (t) (f ) z (f ) z ! (f ) 1
Hence, (t) is written as (t) = t + (jtj) in L 1 (C), where

= 1 j j 2 (f ) z (f ) z ! (f ) 1
Apply Theorem 5 to this family ff t g. Then we can conclude that

f t (&) & t converges to f (&) = 1 Z Z C (z) & (& 1) z (z 1) (z &) dxdy
locally uniformally on C. Hence, changing the variable z in this integrale to (f ) 1 (z) and noting that

f (t) f t = (f t f 0 ) t f
we get the assertion.

  the same conditions as the last proposition. Using the Schwarz re ‡exion principale the map e h h 1 can be extended to an element of Aut (C). Hence with suitable complex numbers c and d, we have e h (z) = ch (z) + d; (z 2 IntQ)

  Then F = e h f h 1 is a qc mapping of the IntR onto the Int e R which maps 0; a, ib and a + ib to 0; e a, i e b and e a + i e b, respectively. In particular, F (z) is ALC on R, and hence for almost every y 2 [0; b], we have e a jF (a + iy) F (iy)j =

  a e b; intgrating both sides of th above inequality over [0; b]

  By the uniqueness theorem, the canonical e -qc mapping f e of C satis…es f e (z) = f e (z)

  for every & 2 C, and the convergence is locally uniform on C.Moreover : f [ ] has the integral representation f [ z (z 1) (z &) dxdy Proof. [IT ] for example. Corollary 3 Let f (t)g be a family of Beltrami coe¢ cients depending on t 2 R or C. Suppose that (t) is di¤erentiable at t = 0 : (t) (z) = (z) + t (z) + t (t) (z) ; z 2 C with suitable 2 B (C) 1 , 2 L 1 (C) such that k (t)k 1 ! 0 as t ! 0: Then f (t) (&) = f (&) + t f [ ] (&) + (jtj) ; & 2 C locally uniformaly on C as t ! 0, where f [ f (&) (f (&) 1) ((f ) z (z)) 2 f (z) (f (z) 1) (f (z) f (&)) dxdy


Yiping Jiang
email: yiping.jiang@student.unsw.edu.au

Jinling Wang
email: bjinling.wang@unsw.edu.au

A New Approach to Calculate the Horizontal Protection Level

Keywords:

A method to compute the minimum Horizontal Protection Level (HPL) using the test statistic of normal distribution, which will exploit advances in computational power to meet the requirement of Time to Alert (TTA), is proposed to improve service availability. To obtain the minimum solution, two approximations used in traditional algorithms need exact solutions: the distribution of the horizontal position error and the determination of the worst case to ensure that the resulted HPL is able to accommodate all possible bias. It is validated with results that the optimal solution is achieved with a pre-defined accuracy and sufficient computational efficiency. Also, the new HPL is used to determine if current approximated methods are conservative, where one of the methods does not meet the integrity requirement with given test statistic, error model and integrity risk definition.

Introduction

Global Positioning System (GPS) is being widely used in safety critical air traffic services, e.g. civil aviation. To guarantee the safety of the navigation solution for these applications, the hazardous analysis is necessary, where the required performances regarding accuracy, integrity, continuity and availability are to be assessed (RTCA DO-229D 2006). The integrity for the non-precision approach (NPA) is of concern here, for which the Receiver Autonomous Integrity Monitoring (RAIM) is used to generate HPL for users to be aware of the level of safety.

There are two major issues when calculating the exact HPL with given integrity risk. First, the distribution of the horizontal position error is too complex to calculate the probability of position error for real time applications. Second, the determination of the worst case is also not straightforward to accommodate all possible bias values. Various methods were proposed to simplify the calculation, resulting in approximated HPLs. Two approximated distributions of the horizontal position error are used in the conventional method including the normal approximation [START_REF] Lee | New techniques relating fault detection and exclusion Performance to GPS Primary Means Integrity Requirements, The 8th international technical meeting of the ION[END_REF] and the chi-squared approximation [START_REF] Ober | Ways to improve RAIM/AAIM availability using position domain performance computations[END_REF]. The former has the disadvantage of underestimating the probability in some cases while the latter one is an overestimate [START_REF] Ober | Ways to improve RAIM/AAIM availability using position domain performance computations[END_REF]. Approximations are also made on the second issue to speed up the process. Without a search for the worst case, specific formulas were given in the conventional methods [START_REF] Brown | GPS RAIM: calculation of thresholds and protection radius using chisquared methods-a geometric approach[END_REF][START_REF] Walter | Weighted RAIM for precision approach[END_REF][START_REF] Brenner | Integrated GPS/Inertial Fault Detection Availability[END_REF]. These approximations can speed up the HPL computation, but they are designed to be conservative with higher protection level than the exact one. The exact HPL value can achieve higher service availability, and therefore this is pursued by solving the two major issues.

A Monte Carlo simulation is commonly used to calculate the probability of position error, which is very timeconsuming. Two ways to calculate the exact position error probability have been identified to be equivalent [START_REF] Johnson | Continuous Univariate Distributions[END_REF][START_REF] Ober | Ways to improve RAIM/AAIM availability using position domain performance computations[END_REF][START_REF] Ober | RAIM Performance: How Algorithms Differ[END_REF]): 1) the integration of the distribution for 2D position error in a quadratic form; 2) the integration of the bivariate normal distribution over a circle. Several implementations are studied and compared in the first approach (Duchesne and Lafaye de Micheaux 2010) and the one in [START_REF] Imhof | Computing the Distribution of Quadratic Forms in Normal Variables[END_REF] showing reliable results is used here. A method was designed in the latter approach [START_REF] Didonato | Integration of the general bivariate Gaussian distribution over an offset circle[END_REF] to speed up the calculation, which was first used in navigation integrity monitoring by [START_REF] Milner | A fast and efficient integrity computation for non-precision approach performance assessment[END_REF]. These two implementations to calculate the probability can achieve a pre-defined accuracy, and they are both referred as the "exact distribution" in this paper.

For the second issue, an iterative search was designed in [START_REF] Milner | Weighted RAIM for APV: The Ideal Protection Level[END_REF] to search the worst case bias with the maximum integrity risk. But it is observed that the iterative method is not able to control the accuracy of the results with uncertainty caused by the number of steps, and computationally heavy with two search loops. Therefore, a new approach is proposed here to improve these two performance criteria based on the method in [START_REF] Jiang | A new approach to calculate the Vertical Protection Level in A-RAIM[END_REF] for the Vertical Protection Level (VPL) computation. A new iterative method with one search loop is designed with results shown that it has higher computational efficiency, but the accuracy of the results is still uncertain. An optimization method is then used to generate results within a pre-defined accuracy, and also the computation is considerably faster than the other methods.

The rest of the paper is organized as follows. The basic model and parameters that are going to be used for calculation of HPL are defined in Section 2. All current methods to calculate HPL are introduced in Section 3. The two issues in calculation of the exact HPL are studied in Section 5 and Section 6 with numerical results in Section 7 and conclusions in Section 8.

Navigational System Model and Integrity Risk

The linearised model to relate GPS observations to the position is,

;

( 1

)
where ∈ is the observation vector which is assumed to be Gaussian; ∈ is the unknown position vector including the three dimensional position error and clock error; ∈ is the design matrix with rank , which is derived by the elevation angle as the transformation from the observation domain to the position domain;

is the covariance matrix of the observations, which is modeled to account for clock and ephemeris errors, troposphere delay, multipath and receiver noise.

The least squares estimate of is,

( 2 ) where , , are the east, north and up component estimation, and is the clock estimation; is the solution estimation matrix with each row as .

RAIM utilizes redundancy in observations to do consistency check. An alarm is sent to the cockpit if an unacceptable large position error is detected to warn the pilot that GPS should not be used for navigation. In the conventional RAIM [START_REF] Brown | GPS RAIM: calculation of thresholds and protection radius using chisquared methods-a geometric approach[END_REF], a chi-squared test is used to detect the failure. A multiple hypothesis structure is used to define the test statistic, which is of normal distribution,

( 3 ) where 1 ⋯ ; ∈ is a zero vector with the i th element as one;

. Under hypothesis :

with as the non-centrality parameter. The probability of false alarm and the probability of missed detection are defined as follows with as the threshold,

| | | ( 4 ) | | | ( 5 )
Protection level is a statistical error bound to guarantee that the probability of the position error exceeding a given bound is within the required integrity risk. If the computed protection level is larger than the given alert limit, an alert within required TTA should be generated by the system. This is an adaptation to the definition in the standards (Annex 10 2006) (DO-229 2006). The position error is the difference between the estimated value and the real one. With HPL of concern, the horizontal position error is a vector with the east position error and the north position error. Also it has been shown that the test statistic and the position error are independent [START_REF] Ober | Ways to improve RAIM/AAIM availability using position domain performance computations[END_REF]. The allocated integrity risk under is defined as [START_REF] Geas | GNSS Evolutionary Architecture Study, GEAS Phase I -Panel Report[END_REF],

‖ ‖ | | | | ( 6 )
where integrity risk is a function of the probability of the horizontal position error ‖ ‖ | , the prior probability of hypothesis and . It should be noted that the integrity risk for each hypothesis is defined to accommodate all possible bias size. In other words, eq. ( 6) is defined to bound the worst case bias by the given integrity risk. The worst case bias is the one that generates the biggest integrity risk. With each computed to protect the worst case bias under each hypothesis, the final HPL is the maximum one among all hypotheses with 0,1 ⋯ .

The risk for each hypothesis is set as [START_REF] Brenner | Integrated GPS/Inertial Fault Detection Availability[END_REF][START_REF] Kelly | The linear model, RNP, and the near optimum fault detection and exclusion algorithm[END_REF]) 1):

=3.33×10 -7 /m, =10 -7 , =10 -4 , and therefore the given value of is 10 -3 . is the fault free alarm probability from the continuity risk as derived in [START_REF] Kelly | The linear model, RNP, and the near optimum fault detection and exclusion algorithm[END_REF]. With the given , the non-centrality parameter in is decided by with eq. ( 4) and (5).

Therefore, is a function of HPL and , and is also a function of HPL and .

The study is aimed at providing a generalized algorithm to calculate HPL for all kinds of navigation service types, since the derivation of the worst case bias is regardless of the service types. The requirement of the NPA is adopted, where HPL is critical. Also, the multiple hypothesis structure is used here, which is also adopted in the design for the next generation RAIM structure ARAIM [START_REF] Pervan | A Multiple Hypothesis Approach to Satellite Navigation Integrity[END_REF]ARAIM report 2012). But The HPL presented here assumes a fixed , as opposed to ARAIM, which can accommodate different prior probabilities. However, it should be noted that the proposed method can be used in a more generalized RAIM scenario.

Current Methods to Calculate HPL

Current methods to calculate HPL are based on approximations due to the complex calculation process to gain the exact probability and the worst case bias. The external reliability method [START_REF] Brown | GPS RAIM: calculation of thresholds and protection radius using chisquared methods-a geometric approach[END_REF] was expanded by a random part to bound the random error [START_REF] Angus | RAIM with Multiple Faults[END_REF]. In this method, two approximated distributions are used to calculate P PE . With the normal approximation, ‖ ‖ is approximated as a one degree Taylor series (Lee 1995), ( 7 ) where and are the non-centrality parameters of and . The approximated variance is of a normal distribution , with which is then projected onto the direction of the faulty observation i with and as the covariance matrix of .

With the normal approximation, HPL BC under is [START_REF] Angus | RAIM with Multiple Faults[END_REF],

1 1 (8)
where is derived by the given and 10 -3 ; K() is the inverse of the cumulative distribution function of a Gaussian random variable with zero mean and unit variance. The slope parameter is defined as the project matrix from the i th observation domain to the horizontal position error domain [START_REF] Brown | GPS RAIM: calculation of thresholds and protection radius using chisquared methods-a geometric approach[END_REF] with , 1

The second approximation is the chi-squared approximation, which is derived by the following inequality [START_REF] Ober | Ways to improve RAIM/AAIM availability using position domain performance computations[END_REF],

‖ ‖ ( 1 0 )
where is of chi-squared distribution with 2 degree of freedom; is the minimum eigenvalue of . With eq. ( 10), it can be concluded that the chi-squared approximation is always safe [START_REF] Ober | Ways to improve RAIM/AAIM availability using position domain performance computations[END_REF].

In a similar way, HPL BC under with the chi-squared approximation is derived as,

2 2,1 (11)
where is derived by the given and 10 -3 ; 2,1 represents the value at probability 1 with the central chi-squared inverse cumulative distribution function and 2 degrees of freedom; The slope factor for each hypothesis is,

2 ( 1 2 )
Another popular method is the weighted RAIM method [START_REF] Walter | Weighted RAIM for precision approach[END_REF] designed for precision approach. The HPL WE is also tested here for NPA,

1 1 ( 1 3 )
where with and as the standard deviation of and respectively, and is derived with the given as 1 .

The solution separation method was proposed by [START_REF] Brenner | Integrated GPS/Inertial Fault Detection Availability[END_REF] and applied in a multiple hypothesis structure [START_REF] Pervan | A Multiple Hypothesis Approach to Satellite Navigation Integrity[END_REF], which is referred as the multiple hypothesis solution separation (MHSS) method,

(

where the east HPL is defined as,

1 ( 1 5 )
where is the standard deviation of the east solution separation, which is equivalent to the east slope parameter The method presented here to calculate HPL can be generalized for multiple faults with the relevant reliability measures for multiple faults described in [START_REF] Knight | Generalised measures of reliability for multiple outliers[END_REF]. With current approximated HPL computation methods introduced as above, the exact HPL is pursued in the following sections, where the distribution of ‖ ‖ is analyzed in Section 4 and the search for the worst case is studied in Section 5.

The Probability of the Horizontal Position Error

Two methods to calculate the exact probability P PE are studied: the Imhof approach [START_REF] Imhof | Computing the Distribution of Quadratic Forms in Normal Variables[END_REF]) and the DJ approach (DiDonato and Jarnagin 1961; Milner and Ochieng 2010). They are chosen in consideration of the ability to control error tolerance for high accuracy, which is set at 10 -9 for both implementations.

The Imhof approach is introduced as follows. With the correlation between the east and north error, ‖ ‖ is not of a chi-squared distribution, but it is a sum of independent chi-squared variables [START_REF] Ober | Ways to improve RAIM/AAIM availability using position domain performance computations[END_REF],

‖ ‖ ∑ ( 1 6 )
where is the observation error vector, which is the difference between the real value and the estimated value; The eigen-decomposition of the symmetric matrix is ⋀ with ⋀ as the diagonal matrix, as the non-zero diagonal elements, as the corresponding eigen-vector, which is the columns of the orthogonal matrix . Therefore, ‖ ‖ is a linear composition of independent chi-squared variables with 1 degree of freedom and non-centrality parameter . The probability of P PE can be computed with the Imhof method,

‖ ‖ | 0.5 ( 1 7 )
where and are,

0.5 ∑ tan 1 0.5 (18) ∏ 1 ( 1 9 )
With the given error tolerance , the upper bound of the upper limit of integral U can be derived by,

∏ ( 2 0 )
With the east and north being de-correlated [START_REF] Milner | A fast and efficient integrity computation for non-precision approach performance assessment[END_REF], the DJ approach is applied, where a rectangular bound is defined to replace the cross area between an ellipse and a circular to reduce the computational load (DiDonato and Jarnagin 1961). The centre position of the circular is the absolute value of the bias position after rotation of the correlation system. The given error tolerance is used to adjust the rectangular bound. The probability obtained by the integration of the cross area between the ellipse with the bivariate normal distribution and a circle is used as a reference which is referred as the Bivariate Normal. Both methods and the two approximated ones (chi-squared, normal) are compared with the reference in Figure 1 with a single epoch example of 7 visible satellites and a 12m outlier on the first observation. As shown in Figure 1, the probabilities based on the Imhof and DJ approaches with 10 -9 error tolerance are almost the same as the Bivariate Normal one. The chi-squared approximation is conservative with larger P PE , while the normal approximation is not. Similar conclusion can be found in [START_REF] Ober | Ways to improve RAIM/AAIM availability using position domain performance computations[END_REF]. Using a standard Intel Core 2 Duo Processor E8400, the computational time per SV to calculate the P PE with different approaches is shown in Table 1. Therefore, the DJ approach should be used with the higher accuracy and computational efficiency. It should also be noted that the Imhof approach is more general and can be used for higher dimensional cases.

The Approaches to Calculate the HPL

The aim is to obtain the exact HPL solution to accommodate all possible bias, which is the HPL solution corresponding to the worst case bias with eq. ( 6). Besides the computation of the exact probability of P PE , the second issue arises when searching for the worst case in a boundary. The maximum HPL within a fixed interval is the corresponding solution. There is a one-to-one relationship between the non-centrality parameter and . The boundary was designed to be in the bias domain from 0 to the Minimal Detectable Bias (MDB) in Milner and Ochieng) 2010, here the domain is used with the boundary as ~1. With the lower boundary given as 10 by eq. ( 6), the new criterion is expressed as, max , subject to 10 1 ( 2 1 )

where stands for the HPL as function of , which is expressed by the non-linear equation ‖ ‖ | 10 with given integrity risk. To obtain the best solution, an iterative method in Milner and Ochieng ( 2010) is introduced first, followed by an improved iterative method, and finally an optimization method is proposed, which are all compared together in terms of accuracy and efficiency.

With a closed boundary, it is guaranteed that there is always an extreme value. The mechanism is shown in Figure 2, with the decrease of within the range, the non-centrality parameter is increased from to . The corresponding HPL increases from HPLa to HPLb and then decreases to HPLc. Therefore, a search is needed to obtain the worst case. To search the worst bias among the range, it is necessary to define the worst case first. As shown in Figure 3, the worst integrity risk happens when is 15.01%, which is the same for the worst HPL in Figure 4. The worst case happens at the same point because the input HPL is 5.197m in Figure 3 corresponding to the worst HPL in Figure 4, and same versa. Still, it is illustrated that there are two ways to define the worst case: the maximum integrity risk with given HPL; the maximum HPL with given integrity risk.

The Original Iterative Method

A procedure is defined in [START_REF] Milner | Weighted RAIM for APV: The Ideal Protection Level[END_REF] with the maximum integrity risk as the worst case to calculate the ideal VPL, which can be applied to HPL as well, and referred to as the iterative A method. There are two iterative search loops over the inner loops over the bias range and the outer one as the input HPL range.

For each inner loop, there is a worst bias point that has the maximum integrity risk with the given input HPL. For all the input HPL values in the outer loop, the one that generates the maximum integrity risk closest to the given integrity risk is the desired value. The process is described as follows: in the outer loop, an upper value U is chosen to decide the range 0 with HPL k as the value at step k. The range 10 1 is used in the inner loop with as the value at step j. In each step j, the integrity risk ‖ ‖ | is calculated and the maximum one among the inner loop is chosen. The HPL k that generates the maximum integrity risk which is the closest to 10 is the final result. Therefore, the accuracy of the results and computational time is dependent on the number of steps in two loops.

The New Iterative Method

The iterative method is designed to derive each HPL value corresponding to each with the non-linear equation and the maximum one is the final result with the worst case bias. The procedure is as follows:

1) the range is divided by the number of steps N with as the value at step j.

2) the corresponding to each are derived by equation

‖ ‖ | 10 ,
where an optimization tool to solve non-linear equations is used. 3) the maximum is therefore the desired result. In this way, the outer loop is removed, and both the accuracy of the results and computational time is dependent on the number of steps in single search loop. This method is referred as the iterative B method.

The Optimization Method

To gain a HPL value with higher accuracy and computational efficiency, an optimization algorithm is used to find the maximum value within a fixed interval. The MATLAB implementation fminbnd is used to find the single local minima with a continuous function of one variable, which is based on golden section search and parabolic interpolation [START_REF] Forsythe | Computer Methods for Mathematical Computations[END_REF]). The criterion is then adapted to be used in this optimization tool, where is treated as the unknown parameter that needs to be solved and HPL is expressed by the unknown with ‖ ‖ | 10 . The maximization problem is adapted as a minimization problem by subtracting the HPL function from a large number (e.g. 100-HPL). The optimization tool is used to find the value that has the maximum HPL value. The termination tolerance, which is the norm of the difference of two consecutive solution vectors during the iterative optimization process, is set as 10 -8 to guarantee the accuracy of the final results. The initial value is set according to experiment data to speed up the solving process. The procedure of this method is described as below:

1) The function to compute is implemented with the DJ approach as . The input parameters include , HPL and the non-central parameter in and as * , * . The output is .

2) The function to compute HPL is implemented as follows:

a. Define an initial value and a search range ∈ 10 , 1 . b. Define other input parameters: , / , and the eigenvalues of . c. A function is defined to compute the equation * 0 with fsolve, where is the input variable and HPL is the output variable. The initial value and the termination tolerance are defined beforehand. d. 100

as a function of variable from is then used in fminbnd as the function to be minimized within a fixed interval of ∈ 10 , 1 . e. The solution from fminbnd is the worst case , and the corresponding HPL is the final solution, which is derived with the equation in .

Although in theory there is always be a solution, the convergence of the optimization algorithm depends on the proper choice of the initial value and the termination tolerance. Besides, to avoid discontinuity, the new iterative method is used in practice whenever there is no extreme value found in optimization algorithm.

Numerical Results

Comparing the new iterative method with the original one, the result would be close to the theoretical value if the steps were close to infinite. However, the number of steps cannot afford to be too big with the requirement of TTA in RAIM. Otherwise the integrity information would not be transmitted to users in time. Normally, TTA should be within a few seconds, for example the TTA requirement for the LPV-200 service is 6.2 seconds. With a smaller number of steps, it is observed that the accuracy of the results is not necessarily better when using more steps. The reason lies in the uncertainty of the location of the worst case and the step size. An example is given in Figure 5 with various steps: 5, 10, 20 and 30. The worst case under each steps are marked with a rectangular box. It can be seen that the HPL using 20 steps is closer to the correct HPL than using 30 steps. With all the uncertainties, it is concluded that the iterative methods are not only too time consuming, but also lack sufficient accuracy. 2, where the HPL is ranged between 0~50m with 50,000 steps to achieve the accuracy of 0.001m in the iterative A method. The result from the optimization method (7.3384m) is used as the reference point. 2407 7.3368 7.3377 7.3319 7.3382 7.3383 7.3384 In Table 2 with the iterative B method, the result with 40 steps is closer to the reference value than with 50 steps. Therefore, it is concluded that the accuracy of the results is not necessarily getting better with larger steps, which is consistent with the analysis in Figure 5. But with extremely large steps of 10000, the value is the same as the reference value. Also, with the number of steps larger than 100, the result stays relatively stable at the accuracy of 10 -4 . With an additional search loop in the iterative A method, the HPL precision is constrained by the outer the number of steps with 0.001m accuracy. Yet, this accuracy is not reliable since it is still susceptible to the uncertainty caused by the search in the range with the example of the results using 40 and 50 steps in Table 2. Therefore, it is concluded that introducing another search loop in the iterative A method does not enable the control of the accuracy of the results.

The iterative B method with 150 steps and 10 steps are compared with the optimization method in Figure 6 using data of the GPS observations captured on UNSW campus within 24 hours. It is assumed that the error is of standard normal distribution without correlation. It is shown in Figure 6 that the results from the iterative method with 150 steps are very close to the results from the optimization method. However, the results from the iterative method with 10 steps are smaller than the optimization method at the scale of 0.1m. Therefore the accuracy of the results from the iterative method is lower and the conservativeness is not guaranteed, depending on the number of search steps, while the optimization method is more reliable.

Using an Intel Core 2 Duo Processor E8400, the average computational time per epoch is shown in Table 3. For both iterative methods, the total search step is 100 in the range. HPL is ranged between 0~50m with 500 steps with 0.1m accuracy. With 100 search steps, the results from the iterative B and optimization method are similar, but it is shown in Table 3 that the optimization method is faster with less than half of the time used per SV than the iterative B method. The iterative A method consumes around 20s per epoch, while the optimization method is more promising to be used in real time applications. Also, an approach to limit the search region was designed in [START_REF] Milner | Determination of the effects of GPS failures on aviation applications[END_REF], as well as [START_REF] Milner | A fast and efficient integrity computation for non-precision approach performance assessment[END_REF] to improve the computational efficiency. But since the two loops impose a considerable computational burden, such improvement is less obvious than the Iterative B method in this test.

All existing approximated algorithms to calculate HPL are examined by comparing with the new optimization value in Figure 7. As shown in Figure 7, the two HPL BC s are conservative, while HPL WE has values higher than the new HPL.

Although the normal approximation distribution is not safe compared with the exact distribution as shown in Figure 1, the HPL BC2 based on it is still conservative.

To quantify the performance improvement with the new HPL, the simulation for worldwide availability distribution is set up as follows. The almanac data for the standard 24 GPS constellation and the and 27 Galileo constellation was used respectively to determine the geometry at each location with a 5×5 degree grid on the world map at 50m altitude. The error model is adopted from GEAS (2008) for ARAIM, and a bias term was added to gain more conservative results. The values for the nominal bias and the maximum bias were 0.1m and 0.75m separately. The User Range Error and User Range Accuracy were 0.25m and 0.5m. The risk definition is the multiple hypothesis requirements as defined in Section 2. Results of HPL at each location are computed every 10min among the 24 hour duration. The mask angle of GPS is set to be 5°. The availability was determined by comparing each HPL value with the alert limit (40m) for each location. The percentage of having over 99% availability over this time is shown worldwide. The simulation software is based on the MATLAB Algorithm Availability Simulation Tool provided by the Stanford University. To examine more obvious improvements, the HAL is set to be arbitrary 35m, It is shown in the above figures and tables the performance is improved with the 99% availability worldwide, especially when the alert limit is tighter in Table 5. Although HPL PB showed lower value than HPL BC1 in Fig. 7, the service availability is higher than HPL BC1 with 24 GPS and lower with 27 Galileo. This is an example to conclude that there is no definite conclusion comparing the size of HPL PB and HPL BC1 , which is dependent on the geometry between satellites and users.

Concluding Remarks

It has been shown that the new approach to calculate HPL has high accuracy and computational efficiency, and therefore service availability can be improved in RAIM and potentially in other GNSS aviation operational applications. Also, the computation of the new iterative method is faster than the old method. With the exact HPL value validated, the conservativeness of current HPLs are analyzed. It is proved that the chi-squared distribution is a safe choice to approximate the exact distribution of horizontal position error, while the normal approximation is not. The results have shown that HPL BC and HPL PB are conservative, while HPL WE does not meet the theoretical integrity requirement. With the optimization method given in the MATLAB toolbox, further efforts can be made to customize an optimization method for this specific problem to maximize the computational efficiency. Although there is no outage of discontinuity with the optimization method in this paper, the proposed method should be further investigated to make sure that 100% convergence can be guaranteed. For the same purpose, the idea proposed by [START_REF] Milner | Determination of the effects of GPS failures on aviation applications[END_REF] to limit the search region may be incorporated in the current proposed method for further study.

(

  [START_REF] Blanch | RAIM with Optimal Integrity and Continuity Allocations under Multiple Failures, Aerospace and Electronic Systems[END_REF]; is the standard deviation of the subset position error , which is derived with the i th observation removed in the estimator; The HPL for the north can be derived in a similar way with is the standard deviation of the east solution separation also equivalent to the north slope parameter | | .

Figure 1 .

 1 Figure 1. P PE as a function of HPL with different distributions

Figure 2 .

 2 Figure 2. Illustration of the search method

  Figure 3.

Figure 5 .

 5 Figure 5. The uncertainty in iterative methods with the number of steps Results from two iterative methods with different steps for range are compared in Table2, where the HPL is ranged between 0~50m with 50,000 steps to achieve the accuracy of 0.001m in the iterative A method. The result from the optimization method (7.3384m) is used as the reference point.

Figure 6 .

 6 Figure 6. The iterative methods and the optimization approach for the new HPL

Figure 7 .

 7 Figure 7. The new HPL and other approximated HPLs

Figure 8 .Figure 10 .

 810 Figure 8. 99% Availability with HPLBC1, 24 GPS Figure 9. 99% Availability with HPLBC2, 24 GPS

  ARAIM report (2012) Interim Report of the EU/US Cooperation on Satellite navigation released, WorkingAvailability with HAL = 40m, Coverage(99%) = 91

Table 1

 1 The average computational time of P PE per SV

		Exact Distributions		Approximated Distributions
	per SV	Bivariate Normal Imhof DJ	Chi-squared	Normal
	Time (ms)	100	85	0.5	0.79	0.35

Table 2

 2 Different HPL values in one epoch with various steps

	Steps	10	30	40	50	100	200	10000
	Iterative A (m) 7.241	7.337	7.338	7.332	7.338	7.338	7.338
	Iterative B(m) 7.

Table 3

 3 The average computational time of HPL in one epoch

		per SV 7SV 8SV 9SV 10SV
	Iterative A(s)	2.60	19.2 20.4 23.2 25.4
	Iterative B (s)	1.30	4.86 9.44 13.8 17.8
	Optimization (s)	0.56	2.31 3.93 5.98 7.66

Table 4 .

 4 99% 40m HAL Availability with Different Algorithms

	RAIM Algorithms HPLBC1 HPLBC2 HPLPB New HPL
	24 GPS	91.52% 89.17% 91.26%	100%
	27 Galileo	100%	100%	100%	100%

Table 5 .

 5 99% 35m HAL Availability with Different Algorithms

	RAIM Algorithms HPLBC1 HPLBC2 HPLPB New HPL
	24 GPS	74.99% 56.88% 83.90%	92.19%
	27 Galileo	89.78% 79.06% 89.06%	95.56%

Availability with HAL = 40m, Coverage(99%) = 91.52% < 50% > 50% > 75% > 85% > 90% > 95% > 99% >99.5% >99.9%