【正文】 ns of ERS SAR Interferometry,1996,~，Zurich.[3] Galloway .,et of aquifersystem paction and land subsidence using interferometric synthetic aperture radar, AntelopeValley, MojaveDesert, California[J]. WaterResources Research，1998，34(10)：25732585.[4] 王超、張紅、劉智，星載合成孔徑雷達干涉測量[M].北京：科學出版社，2002.[5] 劉智、王超、張紅等、張北，尚義地震同震形變場雷達差分干涉測量[J].科學通報,2000,45(23):2550—2553.[6] 張紅、王超、劉智，獲取張北地震同震形變場的差分干涉測量技術[J].中國圖像圖形學報,2002,5(A6):497500.[7] 王超、張紅、劉智，蘇州地區(qū)地面沉降的星載合成孔徑雷達干涉測量監(jiān)測[J].自然科學進展,2002,12(6):621～624.[8] 沈京湘、王鶴、郎剛，區(qū)域高程內插方法研究[J]，北京測繪，2008,(2)，2326.[9] 王繼業(yè)、劉萬宇，GPS水準多項式曲面擬合模型研究[J]，哈爾濱師范大學自然科學學報，2008,(3)，99103.[10] 魏學文，曲面擬合原理在GPS高程測量中的應用[J]．科技情報開發(fā)與經濟，2006， 175176.[11] 張菊清 、劉平芝,基于抗差趨勢面與正交多面函數(shù)的組合擬合法及其在DEM空間數(shù)據(jù)內插中的應用[J]. 測繪學報，2008(04).[12] 黃立人、陶本藻、趙承坤，多面函數(shù)在地殼垂直運動研究中的應用[J].測繪學報1993，22(1)，2129.[13] 劉經南、施闖、姚宜斌、陶本藻，多面函數(shù)擬合法及其在建立中國地殼平面運動速度場模型中的應用研究[J]．武漢大學學報(信息科學版)，2001，(6)，500504.[14] 楊國華，黃立人．速率面擬合法中多面函數(shù)幾個特性的初步數(shù)值研究[J]．地殼形變與地震，1990，(1)，7083.[15] 封建胡，車剛明，聶玉峰．數(shù)值分析原理[M]．北京：科學出版社2005．[16] 孫家廣,計算機圖形學[M].清華大學出版社，1999.附錄 1 多項式內插MATLAB代碼五次多項式擬合代碼A=load(39。F:\matlab\toolbox\nmm\data\dat\39。)。[x,y]=size(A)。k1=sum(A(:,1))/x。 %求出x坐標平均值矩陣k2=sum(A(:,2))/x。 %求出y坐標平均值矩陣k3=(10*sum(A(:,3)))/x。 %求出z坐標平均值矩陣B=load(39。F:\matlab\toolbox\nmm\data\dat\39。)。[a,b]=size(B)。 P=eye(a)。a1=ones(a,1)。 a1(:,1)=k1。a2=ones(a,1)。 a2(:,1)=k2。a3=ones(a,1)。 a3(:,1)=k3。m=B(:,1)a1。 %中心化后的x坐標矩陣n=B(:,2)a2。 %中心化后的y坐標矩陣l=10.*B(:,3)a3。 %中心化后的z坐標矩陣N=[ones(a,1) m n m.^2 m.*n n.^2 m.^3 (m.^2).*n m.*(n.^2) n.^3 m.^4 (m.^3).*n (m.^2).*(n.^2) m.*(n.^3) n.^4 m.^5 (m.^4).*n (m.^3).*(n.^2) (m.^2).*(n.^3) m.*(n.^4) n.^5]。X=inv(N39。*P*N)*N39。*P*l。V=N*X+a310.*B(:,3)。nihejingdu=sqrt(V39。*P*V/(a1))C=load(39。F:\matlab\toolbox\nmm\data\dat\39。)。[s,t]=size(C)。 a3=ones(s,1)。 a3(:,1)=k1。a4=ones(s,1)。 a4(:,1)=k2。a5=ones(s,1)。 a5(:,1)=k3。z=C(:,1)a3。 %中心化后的x坐標矩陣u=C(:,2)a4。 %中心化后的y坐標矩陣f=10.*C(:,3)a5。 %中心化后的z坐標矩陣L=[z u f]。M=[ones(s,1) z u z.^2 z.*u u.^2 z.^3 (z.^2).*u z.*(u.^2) u.^3 z.^4 (z.^3).*u (z.^2).*(u.^2) z.*(u.^3) u.^4 z.^5 (z.^4).*u (z.^3).*(u.^2) (z.^2).*(u.^3) z.*(u.^4) u.^5]。Q=M*X+a510.*C(:,3)。T=eye(s)。jianhejingdu=sqrt(Q39。*Q/(s1)) %計算外符合精度t=1:a。 t0=1:s。subplot(2,1,2)plot(t0,Q,39。+39。)xlabel(39。點號39。)ylabel(39。殘差(mm)39。)title(39。五次多項式檢核殘差圖39。)count=0。 for i=1:s if abs(Q(i,1))10 count=count+1。 text(i,Q(i,1),39。*39。) end endcount。subplot(2,1,1)plot(t,V39。,39。+39。)xlabel(39。點號39。)ylabel(39。殘差(mm)39。)title(39。五次多項式擬合殘差圖39。)count1=0。for i=1:a if abs(V(i,1))5 count1=count1+1。 text(i,V(i,1),39。*39。) endendcount1附錄 2 多面函數(shù)內插MATLAB代碼基于正交最小二乘多面函數(shù)擬合代碼A=load(39。F:\matlab\toolbox\nmm\data\dat\39。)。[x,y]=size(A)。a1=sum(A(:,1))/x。 a2=sum(A(:,2))/x。 B=load(39。F:\matlab\toolbox\nmm\data\dat\39。)。[a,b]=size(B)。r1=ones(a,1)。 r1(:,1)=a1。 r2=ones(a,1)。 r2(:,1)=a2。 m1=B(:,1)r1。 %中心化后的x坐標矩陣n1=B(:,2)r2。 %中心化后的y坐標矩陣l1=10.*B(:,3)。 k1=[m1 n1 l1]。 P=eye(a)。S0=[]。 for i=1:a。 for j=1:a。 Sij=sqrt((k1(i,1)k1(j,1)).^2+(k1(i,2)k1(j,2)).^2)。 S0=[S0。Sij]。 endend S=reshape(S0,a,a)。D=inv(S39。*P*S)*S39。*P*l1。V=S*Dl1。nihejingdu=sqrt(V39。*P*V/(a1))。K=l139。*S。E=sum(abs(l1))。G=[]。for j=1:a。F=sum(abs(S(:,j)))。G=[G。F]。endH=reshape(G,a,1)。I=K39。./(E.*H)。J=abs(I)。g=[]。 for j=1:aif abs(J(j))==max(abs(I))g=[g。j]。 end endg。T=S(:,g)。S(:,g)=[]。D1=inv(T39。*P*T)*T39。*P*l1。C=load(39。F:\matlab\toolbox\nmm\data\dat\39。)。[s,t]=size(C)。r4=ones(s,1)。 r4(:,1)=a1。 r5=ones(s,1)。 r5(:,1)=a2。 r6=ones(s,1)。 r6(:,1)=a3。 m2=C(:,1)r4。 %中心化后的x坐標矩陣n2=C(:,2)r5。 %中心化后的y坐標矩陣l2=10.*C(:,3)。 k2=[m2 n2 l2]。 S1=[]。 for d=1:s。 for j=1:a。Sdj=sqrt((k2(d,1)k1(j,1)).^2+(k2(d,2)k1(j,2)).^2)。 S1=[S1。Sdj]。 end endAA=reshape(S1,s,a)。BB=AA(:,g)。T=BB*Dl2。 Q=eye(s)。 jianhejingdu=sqrt(T39。*Q*T/(s1))t1=1:a。t0=1:s。subplot(2,1,2)plot(t0,T,39。+39。)xlabel(39。點號39。)ylabel(39。殘差(mm)39。)title(39。多面函數(shù)檢核殘差圖39。)count=0。 for i=1:s if abs(T(i,1))13 count=count+1。 text(i,T(i,1),39。*39。) end endcount。subplot(2,1,1)plot(t1,V39。,39。+39。)xlabel(39。點號39。)ylabel(39。殘差(mm)39。)title(39。多面函數(shù)擬合殘差圖39。)count1=0。for i=1:a if abs(V(i,1))10 count1=count1+1。 text(i,V(i,1),39。*39。) endendcount1。附錄 3 地面沉降等值線圖MATLAB代碼A=load(39。F:\matlab\toolbox\nmm\data\dat\39。)。[x,y]=meshgrid(570000:100:600000,3770000:10:3810000)。 z=griddata(A(:,1),A(:,2),A(:,3),x,y,39。v439。)。 subplot(1,1,1)%figure(2) [c,h] =contour(x,y,z,15)。 clabel(c,h) xlabel(39。x39。)。ylabel(39。y