張永紅,張現(xiàn)峰,付 姣,3,金姍姍,4
1.中國測繪科學研究院,北京100830;2.山東正元地理信息工程有限責任公司 濰坊分公司,山東 濰坊261021;3.遼寧工程技術(shù)大學測繪與地理科學學院,遼寧阜新123000;4.山東科技大學測繪科學與工程學院,山東 青島266510
由于成像不受天氣及時間的影響,基于合成孔徑雷達影像的變化檢測技術(shù)在災害管理(如地震、洪澇等自然災害的災情評估與監(jiān)測)、資源調(diào)查(如農(nóng)作物長勢監(jiān)測)等方面具有獨特的優(yōu)勢。但是,SAR固有的斑點噪聲增加了變化檢測處理的復雜性和不確定性[1]?;诠鈱W影像的變化檢測在很多應用領(lǐng)域?qū)崿F(xiàn)了實用化[2],相比之下,對SAR影像變化檢測的研究還遠未成熟。
早期的SAR變化檢測主要是利用單極化SAR影像。由于斑點噪聲為乘性噪聲,為了減弱斑點噪聲的影響,變化檢測主要是利用前后時相SAR 幅度或強度影像的比值[1,3-6],或者是前后時相SAR影像比值對數(shù)[7-8],生成一個包含變化信息的差異影像,然后對差異影像進行閾值處理提取變化信息。利用小波變化對噪聲不敏感的優(yōu)勢,國內(nèi)學者提出了將比值對數(shù)圖像進行小波變換并通過EM算法確定變化閾值的方法[9],以及將小波變換和區(qū)域級決策融合結(jié)合在一起的SAR變化檢測技術(shù)[10-11]。全極化SAR影像包含更為豐富的地物信息,因此在變化檢測方面比單極化更有優(yōu)勢。國內(nèi)外已經(jīng)發(fā)表的全極化變化監(jiān)測方法主要有3類。第1類是將傳統(tǒng)變化檢測技術(shù)移植到極化SAR影像中的,如文獻[12]將主成份分析的思想應用于多時相極化SAR影像變化檢測。這類方法的不足之處是對極化信息使用不充分。第2類是利用一些極化特征進行變化信息提取,例如文獻[13]通過對目標協(xié)方差矩陣的特征值和特征向量提取廣義標準差、相對比率和橢率角等參數(shù)的分析提取變化區(qū)域,文獻[14]基于復相干因子發(fā)展了一個四維向量,通過對此向量的分析提取變化信息。第3類是基于極化SAR影像統(tǒng)計分布的變化檢測[15-18],其中文獻[17—18]提出的基于極化協(xié)方差矩陣的似然比檢驗的變化檢測方法是其中的代表。當極化通道的后向散射系數(shù)或者通道間的相位差/相關(guān)性這兩類指標中的任何一個發(fā)生變化時,該方法都有較好的效果。但是,這類方法的不足是需要對該極化協(xié)方差矩陣做統(tǒng)計分布假設,高分辨率SAR影像的統(tǒng)計分布比較復雜。同時由于似然比統(tǒng)計量的概率密度函數(shù)無法顯式表達,在利用給定的顯著性水平確定拒絕或接受原假設(前后時相協(xié)方差矩陣相等,即無變化)的臨界值時,使用了大樣本條件下的卡方分布近似替代,因此,該方法理論上只在使用了很高視數(shù)(大于30)進行多視處理的數(shù)據(jù)集上或者變化圖斑很大時有效,對于一般多視的極化SAR數(shù)據(jù)或者小圖斑的變化檢測可靠性無法保證。
針對上述情形,本文從用于極化SAR影像分類的Wishart距離出發(fā),提出一種新的適用于變化檢測的極化距離測度,并發(fā)展了相應的極化變化檢測方法,最后利用北京地區(qū)的多時相全極化RADARSAT-2影像進行試驗分析。該方法無需考慮極化SAR影像的統(tǒng)計分布,計算簡便,而且在像元級進行處理,對小的變化圖斑也能檢測。
全極化SAR影像記錄了地面每個分辨單元在4種基本極化狀態(tài),即hh、hv、vh、vv(hh表示水平發(fā)射/水平接收狀態(tài),其他類推)的散射回波的幅度和相位,形成如下散射矩陣S
通常假定介質(zhì)滿足互易性,此時Svh=Shv。
自然界中的地物絕大多數(shù)為分布式目標,對于這類地物,通常用極化相干矩陣T3或者極化協(xié)方差矩陣C3來描述其散射信息[19]。C3矩陣和T3矩陣均為復對稱矩陣,二者之間酉相似。因此,在下面的敘述中只以極化相干矩陣為例,結(jié)論同樣適用于極化協(xié)方差矩陣。
C3及T3矩陣的跡,被稱為span,表達了地物目標總的散射功率
由于span為4個通道的強度之和,因此比任一通道的強度影像都有更低的斑點噪聲[19]。
多視處理的極化SAR影像的相干矩陣T3服從復 Wishart分布[20],此時,根據(jù)最大似然分類原理,可以推導出如下的樣本相干矩陣與類中心相干矩陣之間的Wishart距離
式中,Vm為第m類中心的相干矩陣,為該類中所有像元相干矩陣的平均值;|·|表示取矩陣的行列式,上標-1表示求矩陣的逆。任一像元可根據(jù)其與所有類別中心相干矩陣的Wishart距離最小值判斷其歸屬于哪一類[20-21]。
在此基礎上,文獻[22]提出了兩個類別中心相干矩陣之間的距離測度
式(4)結(jié)果越小,表明兩個相鄰的類別極化散射特性越相近,則可以合并為同一類別,最終可實現(xiàn)極化SAR影像的非監(jiān)督聚類。
變化檢測與分類之間具有一定的聯(lián)系。極化分類是將極化特性相似的像元歸為同一類別,而變化檢測則是要從前后時相影像上找出極化特性相差較大的像元。由于極化相干矩陣包含了地物極化散射的完整信息,因此,可以通過前后時相的極化相干矩陣之間的距離大小來判斷像元是否發(fā)生變化。
但是,式(4)定義的極化距離測度用于變化檢測還有明顯不足。主要是距離定義中的第1部分d1(Vi,Vj)(即)與像元是否變化無關(guān),因為不同類別的極化協(xié)方差矩陣可能有很接近的行列式值,同時這部分的值與SAR影像是否經(jīng)過輻射定標有關(guān)。但是,距離定義中的第2部分d2(Vi,Vj)(即))對變化敏感。當像元在前后時相完全一樣時,該值為常數(shù)3,若有變化發(fā)生,則必然導致前后時相極化相干矩陣的不同,此時該值將明顯大于3。
為了進一步說明這種情況,從本文所用的試驗影像中(已經(jīng)過輻射定標處理)選擇了一些樣本,這些樣本包括在前后時相上地類沒有發(fā)生變化的像元8對及地類改變的像元8對,針對前后時相的相干矩陣,分別計算了d1(Vi,Vj)及d2(Vi,Vj)的值,列于表1。
表1 變化像元(序號1-8)和非變化像元(序號9-16)的d1(Vi,Vj)和d2(Vi,Vj)Tab.1 d1(Vi,Vj)and d2(Vi,Vj)values for change pixels numbered from 1to 8and unchanged pixels numbered from 9to 16
從表1看出,變化像元與非變化像元的d1(Vi,Vj)值相混合,表明其對變化完全不敏感。但是所有變化像元的d2(Vi,Vj)均顯著高于未變化像元。
因此,本文提出如下的適用于變化檢測的極化距離測度
式中,T3,1、T3,2分別表示前后時相像元的極化相干矩陣;減3是為了保證同質(zhì)像元之間的距離為零,以滿足距離測度的非負性。對此極化距離測度的進一步分析如下。相干矩陣為埃爾米特半正定矩陣,可表示為
式中,λi為特征值;ui為相應的特征向量,也是酉矩陣U的列向量;上標H表示共軛轉(zhuǎn)置。經(jīng)推導,上述極化距離測度可表示為
將上述極化距離作為變化檢測的差異圖像。對差異圖像進行閾值分割從而提取變化圖斑是遙感變化檢測中另一關(guān)鍵問題,國內(nèi)外研究者提出了很多方法,其中馬爾科夫隨機場模型方法較有優(yōu)勢。該方法將圖像的變化提取問題轉(zhuǎn)化為圖像標記問題,通過利用馬爾科夫隨機場與吉布斯隨機場之間的等價關(guān)系,將差異圖像中的灰度信息與空間信息融入尋求全局最優(yōu)的變化提取結(jié)果的過程中[23],因此能可靠、準確地提取變化信息。該方法引入了鄰域信息,對SAR斑點噪聲引起的偽變化有很好的抑制作用[24]。因此本文利用馬爾科夫隨機場模型法完成從極化距離圖像中提取最后的變化圖斑。相應的,基于極化距離測度的變化檢測方法的流程如圖1所示。
試驗數(shù)據(jù)為北京西部石景山地區(qū)2009年9月與2010年7月獲取的兩景RADARSAT-2全極化單視復影像(SLC),均為降軌圖像,波束模式為Fine Quad-Polarization。具體參數(shù)如表2所示。
表2 本文所用RADARSAT-2試驗數(shù)據(jù)的基本信息Tab.2 Information on the RADARSAT-2experimental images
為了保證兩幅影像的輻射一致性,首先對兩幅SLC影像分別進行了如下輻射定標處理
式中,Iσ、Qσ分別為定標后像元的實部和虛部;DNI、DNQ為原始單視復影像像元的實部和虛部;Aj是與斜距有關(guān)的增益系數(shù),由附帶的查找表文件(lutSigma.xml)提供,每一列對應一個增益系數(shù)[25]。
輻射定標后對兩幅SAR影像進行精確配準。由于兩幅SAR影像入射角有近5°的差異,在地形起伏較大區(qū)域無法精確配準,因此將原始影像中相對平坦的子區(qū)(大小為1340像元×2140像元)作為試驗區(qū)。為了消除入射角差異造成像元地面分辨率不一致對精確配準的影響,首先分別將其轉(zhuǎn)換成具有方形像元的地距圖像,且方位分辨率與原圖一致。然后分別提取兩幅地距影像的span圖像,通過span圖像得到兩幅影像間的配準多項式。以2009年的span圖像為參考圖像,沿距離向和方位向均勻選擇了20×25個窗口(窗口大小為64像元×64像元,稱為目標窗口),對每一個目標窗口影像,在待配準影像(2010年span影像)的搜索窗口內(nèi)利用最大互相關(guān)系數(shù)準則尋找同名點[26]。為了確定亞像元級的同名點位置,將目標窗口影像和搜索窗口影像均進行了2×2倍過采樣。對所有匹配到的同名點對,利用二次多項式擬合其對應的幾何位置關(guān)系,對擬合偏差大于1個像元的同名點對進行剔除,迭代求得配準多項式。最后,根據(jù)最終得到的配準多項式將2010年地距影像的每一個極化通道的實部和虛部均變換至參考影像的坐標系統(tǒng)內(nèi),利用雙三次采樣方式賦值,完成了兩幅SAR影像每一個極化通道間的精確配準。
輻射定標和精確配準后的兩幅地距極化散射SAR影像依然用極化散射矩陣表達,記為S2009、S2010,然后進行2×2的多視生成各自的極化相干矩陣T3,2009、T3,2010。多視處理不僅抑制了斑點噪聲,而且使得配準精度再提高一倍。為進一步抑制斑點噪聲,對相干矩陣施行了窗口大小為3×3的增強Lee濾波。至此完成了全部的預處理工作。預處理后的兩幅極化SAR影像顯示如圖2,大小為1077行×1140列。圖中假彩色合成的RGB通道分別為
在預處理后影像的基礎上,根據(jù)式(4)生成了2009年9月23日和2010年7月15日兩幅影像間的極化距離圖像dc(T3)。
作為對比,分別利用極化總功率和各個極化通道的回波強度圖像生成了如下的差異影像
上述基于后向散射強度對數(shù)比值的差異性變量是SAR變化檢測的常用度量。
為了定量比較本文提出的極化距離測度和上述基于后向散射強度的差異性度量的變化檢測能力,從預處理后的影像中選取了4個感興趣區(qū)域為樣本。它們分別代表了研究區(qū)內(nèi)4種典型的地表覆蓋變化:水體變?yōu)椴莸?、建筑物變?yōu)槁愕?、草地變?yōu)榻ㄖ?、草地變?yōu)槁愕亍颖镜倪x擇及地表覆蓋類型的判斷均參考了從Google Earth下載的高分辨率光學影像。在每個感興趣區(qū)域中分別手動勾繪了一個地表覆蓋發(fā)生變化的子區(qū)和一個與其相鄰的未變化子區(qū)。然后統(tǒng)計了這些子區(qū)分別在dc(T3)、dspan、dshh、dshv、dsvv等5個差異圖像上的均值和標準差,結(jié)果分別列于圖3~圖6。在圖3—圖6的圖(a)及圖(b)中,紅色邊框標記的為變化子區(qū),藍色或黃色邊框標記的為非變化子區(qū);可分離性統(tǒng)計圖(c)表示在上述5個圖像中的變化區(qū)域和非變化區(qū)域的均值及標準差,均值以帶顏色的點標示,橫坐標1、3、5、7、9分別表示變化區(qū)域的8、10分別表示非變化區(qū)域的
圖1 變化檢測流程圖Fig.1 Flowchart of the proposed change detection method
圖2 完成輻射定標、斜-地轉(zhuǎn)換及精配準等預處理后的極化SAR影像Fig.2 Two pre-processed RADARSAT-2images
圖3 水體變?yōu)椴莸貥颖炯翱煞蛛x性統(tǒng)計Fig.3 Sample pixels changed from water to grassland and the separability analysis
圖4 建筑物變?yōu)槁愕貥颖炯翱煞蛛x性統(tǒng)計Fig.4 Sample pixels changed from building to bare land and the separability analysis
圖5 草地變?yōu)榻ㄖ飿颖炯翱煞蛛x性統(tǒng)計Fig.5 Sample pixels changed from grassland to building and the separability analysis
圖6 草地變?yōu)槁愕貥颖炯翱煞蛛x性統(tǒng)計Fig.6 Sample pixels changed from grassland to bare land and the separability analysis
由圖3—圖6可以看出,本文提出的極化距離測度對上述全部4種類型的地表覆蓋變化都具有較好的分離能力。dspan對水體變?yōu)椴莸?、建筑物變?yōu)槁愕亍⒉莸刈優(yōu)榻ㄖ锏?種變化類型具有較好的分離能力,但無法區(qū)分草地到裸地的變化。其原因是水體和建筑物有明顯的后向散射特點,水體在每一個極化通道上均只有很弱的后向散射,建筑物由于其幾何特點(墻面與地面構(gòu)成二面角散射),在每一個極化通道上均呈現(xiàn)很強的后向散射,因此地物在前或后時相為水體或建筑物之一時,在dspan上都會呈現(xiàn)較大的值;裸地與草地的后向散射相近,差別主要表現(xiàn)在裸地的后向散射中以奇次散射占優(yōu),而草地的后向散射中漫散射占優(yōu),但在極化總功率上非常接近。dshv對水體變?yōu)椴莸亍⒔ㄖ镒優(yōu)槁愕氐葍煞N變化具有較好區(qū)分能力,Shv與回波中的去極化分量有關(guān),反映了地物的去極化能力[27]。水體的后向散射是比較純粹的鏡面反射(奇次散射),裸地也以奇次散射占優(yōu),因此在水體變?yōu)椴莸?、建筑物變?yōu)槁愕剡@兩種變化中均表現(xiàn)出hv通道回波強度的明顯變化。從圖5中可以發(fā)現(xiàn),草地到建筑物的變化在hv通道上也幾乎可以分離開。圖5同時表明,草地到建筑物的變化在hh通道上可以得到很好的區(qū)分。這主要是因為建筑物后向散射的中的絕大部分為偶次散射(二面角散射),因此其hh通道的回波強度遠大于vv及hv通道的回波強度。
式中σi、σj分別為兩類的方差;μi,μj為其均值。
表3定量地顯示出,與基于后向散射強度生成的各種差異圖像相比,極化距離圖像能更好地區(qū)分變化。以極化距離圖像dc(T3)為基礎(圖7(a)),利用基于Markov隨機場模型方法提取了最終的變化像元,結(jié)果如圖7(b)所示。為了對最終變化檢測結(jié)果進行定量評價,選擇圖中一塊區(qū)域,在光學影像輔助下,手工提取了所有的變化圖斑,作為參考影像。圖8(a)為圖7(b)白框區(qū)域的放大,與其對應的變化參考影像顯示于圖8(b)。對分別由生成的變化圖斑(均通過Markov隨機場模型方法提取變化圖斑)進行了精度評價,結(jié)果列于表4,可見由極化距離測度提取的變化圖具有最低的虛警率和漏檢率,總體變化檢測精度達到94.54%(1-5.46%)。
表3 不同差異圖像上變化與未變化類的分離度Tab.3 Divergence between change and unchanged pixels in 5difference images
表4 由不同差異圖像提取的變化圖的精度Tab.4 Accuracy evaluation of change detection from different algorithms(%)
圖7 極化距離圖像及變化檢測圖斑Fig.7 The polarimetric distance image(a)and finally extracted changes(b)
圖8 圖7(b)中白框區(qū)域處理的放大顯示Fig.8 Zoom-in display of the subset in Fig.7 (b)marked by white box and the manually interpreted change patches
在基于極化協(xié)方差矩陣分類的Wishart距離基礎上,本文提出了適用于多時相極化SAR影像變化檢測的計算簡便的距離測度,并給出了相應的完整的極化變化提取技術(shù)流程。利用北京地區(qū)2009年和2010年獲取的兩景RadarSat-2全極化影像進行了試驗,通過與基于SAR后向散射強度得到的多種差異圖像進行對比,極化距離測度能反映全部不同類型的地物變化,而且在極化距離圖像上變化與未變化區(qū)域的分離度更高。這說明,本文提出的極化距離測度在變化檢測應用中具有很好的潛力。
[1] RIGNOT E J M,VAN ZYL J J.Change Detection Techniques for ERS-1SAR Data[J].IEEE Transactions on Geoscience and Remote Sensing,1993,31(4):896-906.
[2] LU D,MAUSEL P,BRONDIZIO E,et al.Change Detection Techniques[J].International Journal of Remote Sensing,2004,25(12):2365-2407.
[3] CIHLAR J,PULTZ T J,GRAY A L.Change Detection with Synthetic Aperture Radar[J].International Journal of Remote Sensing,1992,13(3):401-414.
[4] VILLASENOR J D,F(xiàn)ATLAND D R,HINZMAN L D.Change Detection on Alaska’s North Slope Using Repeat Pass ERS-1SAR Imagery [J]. IEEE Transactions on Geoscience and Remote Sensing,1993,31(1):227-236.
[5] MOSER G,SERPICO S B.Generalized Minimum Error Thresholding for Unsupervised Change Detection from SAR Amplitude Imagery[J].IEEE Transactions on Geoscience and Remote Sensing,2006,44(10):2972-2982.
[6] KIM D J,MOON W M,KIM Y S.Application of TerraSAR-X Data for Emergent Oil-spill Monitoring [J].IEEE Transactions on Geoscience and Remote Sensing,2010,48(2):852-863.
[7] DEKKER R J.Speckle Filtering in Satellite SAR Change Detection Imagery[J].International Journal of Remote Sensing,1998,19(6):1133-1146.
[8] BOVOLO F,BRUZZONE L.A Detail-preserving Scale-driven Approach to Change Detection in Multitemporal SAR Images[J].IEEE Transactions on Geoscience and Remote Sensing,2005,43(12):2963-2972.
[9] HUANG Shiqi,LIU Daizhi,HU Mingxing,et al.Multitemporal SAR Image Change Detection Technique Based on Wavelet Transform[J].Acta Geodaetica et Cartographica Sinica,2010,39(2):180-186.(黃世奇,劉代志,胡明星,等.基于小波變換的多時相SAR圖像變化檢測技術(shù)[J].測繪學報,2010,39(2):180-186.)
[10] WAN Honglin,JIAO Licheng,XIN Fangfang.Interactive Segmentation Technique and Decision Level Fusion Based Change Detection for SAR Images[J].Acta Geodaetica et Cartographica Sinica,2012,41(1):74-80.(萬紅林,焦李成,辛芳芳.基于交互式分割技術(shù)和決策級融合的SAR圖像變化檢測[J].測繪學報,2012,41(1):74-80.)
[11] WAN Honglin,JIAO Licheng,WANG Guiting,et al.A Regionof-interest Level Method for Change Detection in SAR Imagery[J].Acta Geodaetica et Cartographica Sinica,2012,41(2):239-245.(萬紅林,焦李成,王桂婷,等.在感興趣的區(qū)域?qū)用嫔线M行SAR圖像變化檢測的方法研究[J].測繪學報,2012,41(2):239-245.)
[12] BARONTI S,CARLA R,SIGISMONDI S,et al.Principal Component Analysis for Change Detection on Polarimetric Multitemporal SAR Data[C]∥Proceedings of IEEE International Geoscience and Remote Sensing Symposium.Pasadena:IEEE,1994:2152-2154.
[13] KERSTEN P R,LEE J S,AINSWORTH T L.A Comparison of Change Detection Statistics in POLSAR Images[C]∥Proceedings of 2005IEEE International Geoscience and Remote Sensing Symposium.New York:IEEE,2005:4836-4839.
[14] SABRY R.A New Coherency Formalism for Change Detection and Phenomenology in SAR Imagery:A Field Approach[J].IEEE Geoscience and Remote Sensing Letters,2009,6(3):458-462.
[15] HAO Hongmei,ZHANG Yonghong,SHI Haiyan,et al.Application of Test Statistic Method in Fully Polarimetric SAR Change Detection[J].Journal of Remote Sensing,2012,16(3):526-532.(郝洪美,張永紅,石海燕,等.統(tǒng)計假設檢驗方法在全極化SAR變化檢測中的應用[J].遙感學報,2012,16(3):526-532.)
[16] WU Fan,WANG Chao,ZHANG Hong.Residential Areas Extraction in High Resolution SAR Image Based on Texture Features[J].Remote Sensing Technology and Application,2005,20(1):148-152.(吳樊,王超,張紅.基于紋理特征的高分辨率SAR影像居民區(qū)提取[J].遙感技術(shù)與應用,2005,20(1):148-152.)
[17] CONRADSEN K,NIELSEN A A,SCHOU J,et al.Change Detection in Polarimetric SAR Data and the Complex Wishart Distribution[C]∥Proceedings of IEEE 2001International Geoscience and Remote Sensing Symposium.Sydney:IEEE,2001:2628-2630.
[18] CONRADSEN K,NIELSEN A A,SCHOU J,et al.A Test Statistic in the Complex Wishart Distribution and Its Application to Change Detection in Polarimetric SAR Data[J].IEEE Transactions on Geoscience and Remote Sensing,2003,41(1):4-19.
[19] LEE J S,POTTIER E.Polarimetric Radar Imaging:From Basics to Applications [M ].Boca Raton: CRC Press,2009.
[20] LEE J S,GRUNES M R,KWOK R.Classification of Multilook Polarimetric SAR Imagery Based on the Complex Wishart Distribution[J].International Journal of Remote Sensing,1994,15(11):2299-2311.
[21] FAMIL L F,POTTIER E,LEE J S.Unsupervised Classification of Multifrequency and Fully Polarimetric SAR Images Based on the H/A/Alpha-Wishart Classifier[J].IEEE Transactions on Geoscience and Remote Sensing,2001,39(11):2332-2342.
[22] LEE J S.Unsupervised Terrain Classification Preserving Polarimetric Scattering Characteristics[J].IEEE Transactions on Geoscience and Remote Sensing,2004,42(4):722-731.
[23] MOSER G,SERPICO S,VERNAZZA G.Unsupervised Change Detection from Multichannel SAR Images[J].IEEE Geoscience and Remote Sensing Letters,2007,4(2):278-282.
[24] KASETKASEM T,VARSHNEY P K.An Image Change Detection Algorithm Based on MARKOV Radom Filed Models[J].IEEE Transactions on Geoscience and Remote Sensing,2002,40(8):1815-1823.
[25] MDA.RADARSAT-2Product Format Definition:RN-RP-51-2713:Issue 1/7[R].Richmond:MDA,2008.
[26] LIU Guoxiang,DING Xiaoli,LI Zhilin,et al.Co-registration of Satellite SAR Complex Images[J].Acta Geodaetica et Cartographica Sinica,2001,30(1):60-66.(劉國祥,丁曉利,李志林,等.星載SAR復數(shù)圖像的配準[J].測繪學報,2001,30(1):60-66.)
[27] ZEBKER H A,VAN ZYL J J.Imaging Radar Polarimetry:A Review[J].Proceedings of IEEE,1991,79(11):1583-1606.
[28] SWAIN P H,DAVIS S M.Remote Sensing:The Quantitative Approach [M ].New York:McGraw Hill Book Company,1978.