PCA基本步骤：
对数据进行归一化处理（直接减去均值）
计算归一化后的数据集的协方差矩阵
计算协方差矩阵的特征值和特征向量
保留最重要的k个特征（通常k要小于n），也可以自己制定，也可以选择一个阈值，然后通过前k个特征值之和减去后面n-k个特征值之和大于这个阈值，则选择这个k
找出k个特征值对应的特征向量
将m * n的数据集乘以k个n维的特征向量的特征向量（n * k）,得到最后降维的数据。

1、进行归一化处理、计算协方差等
def meanX(dataX):
return np.mean(dataX, axis=0) # axis=0表示按照列来求均值，如果输入list,则axis=1
def variance(X):
m, n = np.shape(X)
mu = meanX(X)
muAll = np.tile(mu, (m, 1))
X1 = X - muAll
variance = 1. / m * np.diag(X1.T * X1)
return variance
def normalize(X):
m, n = np.shape(X)
mu = meanX(X)
muAll = np.tile(mu, (m, 1))
X1 = X - muAll
X2 = np.tile(np.diag(X.T * X), (m, 1))
XNorm = X1 / X2
return XNorm
2、定义PCA函数，实现主成分分析
def pca(XMat, k):
average = meanX(XMat) #计算均值,为进行归一化做准备
m, n = np.shape(XMat) #行，列（7）
data_adjust = []
avgs = np.tile(average, (m, 1))
covX = np.cov(data_adjust.T) # 计算协方差矩阵
featValue, featVec = np.linalg.eig(covX) # 求解协方差矩阵的特征值和特征向量
index = np.argsort(-featValue) # 按照featValue进行从大到小排序
#确定主成分的个数
sum=0
max=0
for i in range(len(index)):
sum=sum+featValue[index[i]]
if max/sum<0.99:
max=max+featValue[index[i]]
k=k+1
finalData = []
print ("k = ",k)
if k > n:
print (“k must lower than feature number”)
return
else:
# 注意特征向量是列向量，而numpy的二维矩阵(数组)a[m]
中，a[1]表示第1行值
selectVec = np.matrix(featVec[index[:k]]) # 找出k个特征值对应的特征向量
reconData = (data_adjust * selectVec.T) + average[:,6]
#将m * n的数据集乘以k个n维的特征向量的特征向量（m * k）,得到最后降维的数据。
print( np.shape(reconData))
array = []
array = np.array(array) # 列表转数组
for j in range(k):
array = np.append(array, reconData[:,j]) #将reconData从二维（3列）数组转化从一维数组（1维）
print (np.shape(array))
array1 = array.reshape(k, im_height, im_width) #将 一维数组 转成 3维矩阵（k,高（行），宽（列））
out = ga.SaveArray(array1, os.path.join(path, “after.img”), format=“GTiff”, prototype=img)
return reconData

3、打开img或者tif格式的图像进行PCA
path = “… 主成分分析”
input = os.listdir(path)
data = []
for i in input:
if i == ‘before.img’:
img = gdal.Open(os.path.join(path, i))
im_width = img.RasterXSize # 栅格矩阵的列数
im_height = img.RasterYSize # 栅格矩阵的行数
im_bands = img.RasterCount # 波段数
im_data = img.ReadAsArray(0, 0, im_width, im_height) # 获取数据 将数据写成数组，对应栅格矩阵
#im_data:(波段数，行，列)
for j in range(im_bands):
picture = im_data[j,:,:].flatten() # 变成一维数组
data.append(picture) #将各个波段合并在一个数组里，一行为一个波段（n个波段n行）
print (“data维数”, np.shape(data))
data = np.mat(data) #数组转换为矩阵
data1 = data.T #进行转置
print (“data1维数”, np.shape(data1))
XMat = data1

-得到XMat：是一个numpy的矩阵格式，行表示样本数，列表示特征（每一列表示一个波段的信息，共n个波段）

k = 0 # k：表示取前k个特征值对应的特征向量
pca(XMat,k)